使用GDAL读取Sentinel数据

深碍√TFBOYSˉ_ 2022-05-15 12:52 218阅读 0赞

# 使用GDAL读取Sentinel数据 #

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

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

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

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

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

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

![使用GDAL命令行读取Sentinel数据][GDAL_Sentinel]

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

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

![使用GDAL命令行读取Sentinel数据的元数据信息][GDAL_Sentinel 1]

## 使用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

[HDF_NetCDF]: https://theonegis.github.io/geos/%E5%8F%96HDF%E6%88%96%E8%80%85NetCDF%E6%A0%BC%E5%BC%8F%E7%9A%84%E6%A0%85%E6%A0%BC%E6%95%B0%E6%8D%AE/
[GDAL_Sentinel]: /images/20220515/8d3afca88a08447c9f19c4588bf83085.png
[GDAL_Sentinel 1]: /images/20220515/bba6ead098cf460aa0e4321c1a7f06ee.png