使用GDAL读取Sentinel数据

GDAL 2.1已经原生支持对于Sentinel数据的读取，我这里使用Sentinel-2光学卫星数据给出使用GDAL工具对其进行读取的方法。

这里我们要大概知道Sentinel数据的组织。下载下来的Sentinel数据是一个ZIP压缩包，里面包含了JPEG2000格式的影像数据以及一些XML格式的元数据文件。

GDAL将Sentinel数据看做一个数据集（概念上类似HDF格式的数据集），里面包含了很多子数据文件。所以，对于Sentinel数据的读取就和对于HDF数据的读取是相同的啦。

对于HDF或者NetCDF格式数据的读取参考我的博文：读取HDF或者NetCDF格式的栅格数据

使用GDAL命令行读取Sentinel数据的元数据信息

直接使用gdalinfo [文件名]可以查看Sentinel文件的元信息，如下图所示：

[外链图片转存中...(img-6pmZXhFt-1694159802233)]

从上面的图中我们可以看到所有的Subdatasets的文件全名，这样我们可以继续使用gdalinfo [子数据集全路径]的方式查看具体的子数据集的元数据信息

下图显示的数据子集中包含四个波段的数据（红，绿，蓝，近红外）

[外链图片转存中...(img-8HQ2DbPN-1694159802234)]

使用GDAL命令行工具将Sentinel数据转为GeoTIFF格式

转换是针对具体的子数据集而言的，所以使用gdal_translate [sentinel subdataset full name] [output filename]命令进行

下面的例子将包含红绿蓝近红外波段的数据子集转为GeoTIFF影像

gdal_translate SENTINEL2_L1C:/vsizip/S2A_MSIL1C_20180504T173911_N0206_R098_T13TGF_20180504T212111.zip/S2A_MSIL1C_20180504T173911_N0206_R098_T13TGF_20180504T212111.SAFE/MTD_MSIL1C.xml:10m:EPSG_32613 B2-3-4-8.tif

使用Python脚本读取Sentinel数据

from osgeo import gdal

import os
os.environ['CPL_ZIP_ENCODING'] = 'UTF-8'

filename = ('/Users/tanzhenyu/Desktop/'
            'S2A_MSIL1C_20180504T173911_N0206_R098_T13TGF_20180504T212111.zip')
# 打开栅格数据集
root_ds = gdal.Open(filename)
# 返回结果是一个list，list中的每个元素是一个tuple，每个tuple中包含了对数据集的路径，元数据等的描述信息
# tuple中的第一个元素描述的是数据子集的全路径
ds_list = root_ds.GetSubDatasets()

visual_ds = gdal.Open(ds_list[0][0])  # 取出第12个数据子集（MODIS反射率产品的第一个波段）
visual_arr = visual_ds.ReadAsArray()  # 将数据集中的数据转为ndarray
del visual_arr

# 获得栅格数据的一些重要信息
print(f'打开数据为：{ds_list[0][1]}')
print(f'投影信息：{visual_ds.GetProjection()}')
print(f'栅格波段数：{visual_ds.RasterCount}')
print(f'栅格列数（宽度）：{visual_ds.RasterXSize}')
print(f'栅格行数（高度）：{visual_ds.RasterYSize}')

程序输出如下：

打开数据为：Bands B2, B3, B4, B8 with 10m resolution, UTM 13N
投影信息：PROJCS["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32613"]]
栅格波段数：4
栅格列数（宽度）：10980
栅格行数（高度）：10980