前言
Sentinel-2数据的免费获取无疑为广大遥感、地信等领域的科技工作者提供了强大的数据支撑，然而往往一些项目/工程需要对某一地区进行长时序变化分析，自然也就需要更多景遥感影像支持。为了数据处理的方便，本篇博客主要介绍了利用Python来实现Sentinel-2 L1C级影像批量大气校正并另存为Geotiff格式，而且还可以生成影像的边界范围Shp，以便为后续操作提供便利。
一、Sentinel-1数据介绍
关于Sentinel-2数据各个波段的介绍在网上已有较多的资料，不再过多赘述。此处仅对sentinel-2文件的命名规则进行简单介绍。
sentinel-2文件的命名形式如下：
MMM _ MSIL1C (or L2A) _ YYYYMMDDHHMMSS _ Nxxyy _ ROOO _ Txxxxx _ < Product Discriminator >.SAFE
| 形式 | 含义 |
|---|---|
| MMM | 任务ID（S2A / S2B） |
| MSIL1C (or L2A) | Level-1C (Level-2A) 产品等级 |
| YYYYMMDDHHMMSS | 数据库感应开始时间 |
| Nxxyy | 处理基线编号 |
| ROOO | 相对轨道编号（R001-R143） |
| Txxxxx | 拼接域编号 |
| < Product Discriminator > | 数据开始生产时间 |
| .SAFE | 产品格式（欧洲标准存档格式） |
例如S2B_MSIL1C_20220715T051659_N0400_R062_T44TQR_20220715T071617.SAFE表示标识Sentinel-2B于2022年07月15日凌晨05:16:59（欧洲时间）获得的Level-1C产品。它是在相对轨道062期间通过Tile 44TQR 获得的，并且用PDGS处理基线04.00处理。(The Payload Data Ground Segment (PDGS) :有效载荷数据地面站 )
再补充一点关于Sentinel-2产品级别的信息：
| 产品级别 | 含义 |
|---|---|
| Level-0 | 原始数据 |
| Level-1A | 包含元信息的几何粗校正产品 |
| Level-1B | 包含元信息的几何粗校正产品 |
| Level-1C | 经正射校正和亚像元级几何精校正后的大气表观反射率产品 |
| Level-2A | 主要包含经过过大气校正的大气底层反射率数据 |
目前，欧空局（ESA）仅发布了哨兵2号（S2）的L1C级和L2A级多光谱数据（MSI）。如上表所示，Sentinel-2 L1C是经过正射校正和几何精校正的大气表观反射率产品，并没有进行大气校正。L2A级数据主要包含经过大气校正的大气底层反射率数据（Bottom-of-Atmosphere corrected reflectance），由插件Sen2cor生产得到。
作为用户，我们在使用中可以图省事下载L2A级数据，然而欧空局并不会将云量比较大（Sentinel-2 成像自带检测云量的算法）的L1C处理成L2A级产品，尽管多数情况下，我们的确不需要云量过量的产品。但是这也难免会影响基于影像的对某一地区进行长时序的变化分析。
因此，接下来主要介绍一下利用Python实现Sentinel-2 L1C级影像批量大气校正并另存为Geotiff格式，从而生成L2A级产品的生产流程。
二、Sentinel-2 L2A级批量化生产
1. 环境配置
① 安装Sen2Cor插件
ESA提供了一个将L1C级数据转换为L2A级数据的离线软件Sen2Cor，软件下载地址为: http://step.esa.int/main/snap-supported-plugins/sen2cor/
可以看到，网站提供多个版本的Sen2Cor下载，本文选择下载Sen2Cor_v2.10。安装步骤如下：
- 将下载的文件放到想要安装的位置并解压；
- 解压后双击L2A_Process.bat文件完成软件安装；
- 按Win+R打开CMD命令行窗口，输入命令“cd + 软件安装绝对路径”将路径设置到程序存放目录下，输入下列命令查看程序使用说明，程序不报错，即表示安装成功。
“L2A_Process --help”
② 安装GDAL库
本文主要使用GDAL库读取Sen2Cor所生成的L2A.SAFE文件，并重新写为.tif格式，方便后续处理。（注意GDAL的版本，2.4以下的版本可能不支持）。
这里贴一下GDAL库whl文件的下载链接: https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal
2. 代码运行
以下为 利用Python实现Sentinel-2 L1C级影像批量大气校正并另存为Geotiff格式的完整代码：
# -*- coding: utf-8 -*-
"""
首先读取Sentinel-2 L1C级产.zip文件，将其解压之后利用Sen2Cor进行大气校正；
之后读取得到的Sentinel-2 L2A级产品，删选波段进行波段合成并转为Geo_tiff文件；
最后利用得到的Geo_tiff文件提取其轮廓范围，生成shapefile文件。
"""
import subprocess
import zipfile
import os
from osgeo import gdal
import os
import numpy as np
from osgeo import gdal, osr, ogr
import glob
import fnmatch
import os
# 第一步：解压下载的Sentinel L1C级产品ZIP文件
# 定义一个解压函数
def unzip_file(zip_file_name, mode='rb'):
"""
Return a unzipped file name which comes from the zipped file of A Sentinel-2 L1C data
返回Sentinel-2 L1C级.zip文件解压后的文件名
zip_file_name - the zipped file name of A Sentinel-2 L1C data
mode - Optional parameter, the mode of the zipped file reader
"""
# 打开.zip文件
zip_file = open(zip_file_name, mode)
# 利用zipfile包解压文件
zip_fn = zipfile.ZipFile(zip_file)
# 获取.zip文件的所有子目录名和文件名
namelist = zip_fn.namelist()
for item in namelist:
# 提取.zip文件的子目录及文件，解压在当前文件夹（'.'表示当前文件夹）
zip_fn.extract(item, '.')
# 关闭.zip文件
zip_fn.close()
zip_file.close()
# 打印解压完成
print("Congratulations! Unzipping finished!")
# 返回解压后的文件名（字符串，带.SAFE后缀）
return namelist[0]
# 第二步 利用Sen2cor.bat插件实现Sentinel L1C级产品大气校正
# Sen2cor.bat文件所在路径
sen2cor_path = r"C:\Users\XJG\.snap\auxdata\Sen2Cor-02.10.01-win64\L2A_Process.bat"
# Sentinel-2 L1C原始数据（.zip格式文件，无需解压）所在目录
origin_dir = r"D:\111111"
# 文件模式
pattern = ".zip"
# 输出目录
output_dir = r"D:\111111"
for in_file in os.listdir(origin_dir):
# 判断是否是.zip文件
if pattern in in_file:
# 获取.zip文件的完整路径
zip_file_path = os.path.join(origin_dir, in_file)
# 解压.zip文件，产生.safe格式文件
safe_in_file_path = unzip_file(zip_file_path)
# 设置Sen2cor命令行参数，按照原始分辨率处理
cmd_args = [sen2cor_path, safe_in_file_path, \
'--output_dir', output_dir]
# 打印处理开始的消息
print("{} processing begin!".format(safe_in_file_path))
# 传入命令行参数并调用命令行(cmd)执行命令
subprocess.call(cmd_args)
# 打印处理完成的消息
print("{} processing finished!\n".format(safe_in_file_path))
# 打印所有文件处理完成的消息
print("Congratulations! All zipped file finished!")
# 第三步 提取不同分辨率的波段，并进行相同分辨率波段合成，输出为Geo_TIFF文件
def S2tif(filename):
# 打开栅格数据集
print(filename)
root_ds = gdal.Open(filename)
print(type(root_ds))
# 返回结果是一个list，list中的每个元素是一个tuple，每个tuple中包含了对数据集的路径，元数据等的描述信息
# tuple中的第一个元素描述的是数据子集的全路径
ds_list = root_ds.GetSubDatasets() # 获取子数据集。该数据以数据集形式存储且以子数据集形式组织
visual_ds = gdal.Open(ds_list[0][0]) # 打开第1个数据子集的路径。ds_list有4个子集，内部前段是路径，后段是数据信息。
visual_arr = visual_ds.ReadAsArray() # 将数据集中的数据读取为ndarray
# 创建.tif文件
band_count = visual_ds.RasterCount # 波段数
xsize = visual_ds.RasterXSize
ysize = visual_ds.RasterYSize
out_tif_name = filename.split(".SAFE")[0] + ".tif"
driver = gdal.GetDriverByName("GTiff")
out_tif = driver.Create(out_tif_name, xsize, ysize, band_count, gdal.GDT_Float32)
out_tif.SetProjection(visual_ds.GetProjection()) # 设置投影坐标
out_tif.SetGeoTransform(visual_ds.GetGeoTransform())
for index, band in enumerate(visual_arr):
band = np.array([band])
for i in range(len(band[:])):
# 数据写出
out_tif.GetRasterBand(index + 1).WriteArray(band[i]) # 将每个波段的数据写入内存，此时没有写入硬盘
out_tif.FlushCache() # 最终将数据写入硬盘
out_tif = None # 注意必须关闭tif文件
if __name__ == "__main__":
from osgeo import gdal
SAFE_Path = (r'D:\111111')
data_list = glob.glob(SAFE_Path + "\\*.SAFE")
for i in range(len(data_list)):
data_path = data_list[i]
filename = data_path + "\\MTD_MSIL2A.xml"
S2tif(filename)
print(data_path + "-----转tif成功")
print("Congratulations! ----转换结束----")
# 第四步 提取.tif影像的范围轮廓，生成shapefile文件
ogr.RegisterAll()
shape_path = r"D:\111111\shp" # shape输出位置
if not os.path.exists(shape_path):
os.mkdir(shape_path)
img_list = fnmatch.filter(os.listdir(origin_dir), '*.tif')
for img in img_list:
p_img=origin_dir+'/'+img
outfilename = shape_path+'/'+img[:-4]+".shp"
dataset = gdal.Open(p_img)
oDriver = ogr.GetDriverByName('ESRI Shapefile')
oDS = oDriver.CreateDataSource(outfilename)
srs = osr.SpatialReference(wkt=dataset.GetProjection())
geocd = dataset.GetGeoTransform()
oLayer = oDS.CreateLayer("polygon", srs, ogr.wkbPolygon)
oDefn = oLayer.GetLayerDefn()
row = dataset.RasterXSize
line = dataset.RasterYSize
geoxmin = geocd[0]
geoymin = geocd[3]
geoxmax = geocd[0] + (row) * geocd[1] + (line) * geocd[2]
geoymax = geocd[3] + (row) * geocd[4] + (line) * geocd[5]
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(geoxmin, geoymin)
ring.AddPoint(geoxmax, geoymin)
ring.AddPoint(geoxmax, geoymax)
ring.AddPoint(geoxmin, geoymax)
ring.CloseRings()
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
outfeat = ogr.Feature(oDefn)
outfeat.SetGeometry(poly)
oLayer.CreateFeature(outfeat)
outfeat = None
oDS.Destroy()
# 打印所有文件处理完成的消息
print("Congratulations! All file finished!")
具体说明详见代码中的注释。其中需要注意的是，文件的打开方式与其他格式的遥感影像相同，用gdal.Open(filename)即可读入。返回结果是一个list，list中的每个元素是一个tuple，每个tuple中包含了对数据集的路径，元数据等的描述信息。获取其中的子数据集，使用data.GetSubDatasets()，这里面就包含了子数据集的路径，元数据等的描述信息。然后用子数据集的路径打开，读取为ndarry格式。这里第1个子集存储的是10m分辨率的B2, B3, B4, B8这4个波段，如果需要其他波段，就改ds_list[0][0]换别的子集，如ds_list[1][0]是第2个子集的路径，ds_list[1][1]是第2个子集的描述信息。
ds_list = root_ds.GetSubDatasets() # 获取子数据集信息
visual_ds = gdal.Open(ds_list[0][0]) # 打开第1个数据子集的路径。
visual_arr = visual_ds.ReadAsArray() # 将数据集中的数据读取为ndarray
因此，假如要获得20m分辨率的波段集合，只需要修改这一句：
visual_ds = gdal.Open(ds_list[1][0])
==========================================================================
运行成功后就会生成如下文件：
其中，“.zip”文件为所下载的Sentinel-2 L1C级数据，“.SAFE”文件为利用Sen2Cor生成的大气校正后的Sentinel-2 L2A级数据，“.tif”文件为利用GDAL筛选的分辨率为10m的波段合成影像，而“shp”文件夹内即为所生成的影像边界范围矢量文件。可以在envi中打开影像及shp文件，如下：
参考资料
本篇博客在撰写过程中主要参考了以下几位大神的资料，在此拜谢！
- 哨兵2号（Sentinel-2）介绍、下载、预处理及批处理
- 哨兵2数据‘文件名’详解（吐血整体，仅供参考）
- 02-Sentinel-2 L1C级数据bat和Python脚本批量大气校正
- Python脚本批量读取哨兵2号（Sentinel2）影像并另存为Geotiff格式
噢耶(o)/！
2532
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
