由于Sentinel影像的存儲格式特殊，通常需要特殊的軟件（如SNAP）來處理，有局限性。于是使用GDAL對其進行讀取，并重新寫為.tif格式，方便進一步處理。
影像數據用的是大氣校正后的L2A級數據，當然L1C數據也可以。打開SAFE后綴的文件夾，可以看到里面的內容，jp2格式的圖片數據是存儲在GRANULE文件夾中，這里不需要知道它的存儲規則，因為GDAL可以直接讀取最下面的.xml文件（注意GDAL的版本，2.4以下的版本可能不支持）。

文件的打開方式與其他格式的遙感影像相同，用gdal.Open(filename)即可讀入。返回結果是一個list，list中的每個元素是一個tuple，每個tuple中包含了對數據集的路徑，元數據等的描述信息。獲取其中的子數據集，使用data.GetSubDatasets()，這里面就包含了子數據集的路徑，元數據等的描述信息。

ds_list = root_ds.GetSubDatasets() # 獲取子數據集信息

然后用子數據集的路徑打開，讀取為ndarry格式。這里第1個子集存儲的是10m分辨率的B2, B3, B4, B8這4個波段，如果需要其他波段，就改ds_list[0][0]換別的子集，如ds_list[1][0]是第2個子集的路徑，ds_list[1][1]是第2個子集的描述信息。

visual_ds = gdal.Open(ds_list[0][0]) # 打開第1個數據子集的路徑。 visual_arr = visual_ds.ReadAsArray() # 將數據集中的數據讀取為ndarray

然后要寫新的tif圖像，包括圖像數據和投影等，具體可以參考以下代碼。

# -*- coding: utf-8 -*- # @File : S2_read.py # @Author: Freezinghot # @Date : 2021/3/9 # @Desc : 讀取Sentinel2并轉為Geo_tiff from osgeo import gdal import os import numpy as np from osgeo import gdal, osr, ogr import glob # os.environ['CPL_ZIP_ENCODING'] = 'UTF-8'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.RasterXSizeysize = visual_ds.RasterYSizeout_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 gdalSAFE_Path = (r'E:\RSDATA\Sentinel2\L2A')data_list = glob.glob(SAFE_Path + "\\*.SAFE")#filename = ('E:\\RSDATA\\Sentinel2\\L2A\\S2A_MSIL2A_20210220T024731_N9999_R132_T51STA_20210306T024402.SAFE\\MTD_MSIL2A.xml')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("----轉換結束----")

轉換結束，在目標路徑下生成同名的tif圖像。如果有多幅圖像，也可以實現批量處理。

可以直接拖到ENVI查看波段組合的效果，是保存了4個10m分辨率的波段。

10m分辨率真彩色：

—————————————————————————————————————————————
增加內容：
如果需要其他波段的數據又不知道數據結構的話，可以debug模式看一下文件的結構。下面是我找到的子數據集的結構，位于程序中ds_list中。如圖所示，假如要獲得20m分辨率的波段集合，只需要修改這一句：

visual_ds = gdal.Open(ds_list[1][0])

20m分辨率（B12,B11,B8A）組合：

歡迎在評論區交流，如果對您有幫助，請點個贊哦~

總結

以上是生活随笔為你收集整理的Python脚本批量读取哨兵2号（Sentinel2）影像并另存为Geotiff格式的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。