由于Sentinel影像的存儲(chǔ)格式特殊，通常需要特殊的軟件（如SNAP）來(lái)處理，有局限性。于是使用GDAL對(duì)其進(jìn)行讀取，并重新寫(xiě)為.tif格式，方便進(jìn)一步處理。
影像數(shù)據(jù)用的是大氣校正后的L2A級(jí)數(shù)據(jù)，當(dāng)然L1C數(shù)據(jù)也可以。打開(kāi)SAFE后綴的文件夾，可以看到里面的內(nèi)容，jp2格式的圖片數(shù)據(jù)是存儲(chǔ)在GRANULE文件夾中，這里不需要知道它的存儲(chǔ)規(guī)則，因?yàn)镚DAL可以直接讀取最下面的.xml文件（注意GDAL的版本，2.4以下的版本可能不支持）。

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

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

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

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

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

# -*- coding: utf-8 -*- # @File : S2_read.py # @Author: Freezinghot # @Date : 2021/3/9 # @Desc : 讀取Sentinel2并轉(zhuǎn)為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):# 打開(kāi)柵格數(shù)據(jù)集print(filename)root_ds = gdal.Open(filename)print(type(root_ds))# 返回結(jié)果是一個(gè)list，list中的每個(gè)元素是一個(gè)tuple，每個(gè)tuple中包含了對(duì)數(shù)據(jù)集的路徑，元數(shù)據(jù)等的描述信息# tuple中的第一個(gè)元素描述的是數(shù)據(jù)子集的全路徑ds_list = root_ds.GetSubDatasets() # 獲取子數(shù)據(jù)集。該數(shù)據(jù)以數(shù)據(jù)集形式存儲(chǔ)且以子數(shù)據(jù)集形式組織visual_ds = gdal.Open(ds_list[0][0]) # 打開(kāi)第1個(gè)數(shù)據(jù)子集的路徑。ds_list有4個(gè)子集，內(nèi)部前段是路徑，后段是數(shù)據(jù)信息visual_arr = visual_ds.ReadAsArray() # 將數(shù)據(jù)集中的數(shù)據(jù)讀取為ndarray# 創(chuàng)建.tif文件band_count = visual_ds.RasterCount # 波段數(shù)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()) # 設(shè)置投影坐標(biāo)out_tif.SetGeoTransform(visual_ds.GetGeoTransform())for index, band in enumerate(visual_arr):band = np.array([band])for i in range(len(band[:])):# 數(shù)據(jù)寫(xiě)出out_tif.GetRasterBand(index + 1).WriteArray(band[i]) # 將每個(gè)波段的數(shù)據(jù)寫(xiě)入內(nèi)存，此時(shí)沒(méi)有寫(xiě)入硬盤(pán)out_tif.FlushCache() # 最終將數(shù)據(jù)寫(xiě)入硬盤(pán)out_tif = None # 注意必須關(guān)閉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 + "-----轉(zhuǎn)tif成功")print("----轉(zhuǎn)換結(jié)束----")

轉(zhuǎn)換結(jié)束，在目標(biāo)路徑下生成同名的tif圖像。如果有多幅圖像，也可以實(shí)現(xiàn)批量處理。

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

10m分辨率真彩色：

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

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

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

歡迎在評(píng)論區(qū)交流，如果對(duì)您有幫助，請(qǐng)點(diǎn)個(gè)贊哦~

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。