sentinel2波段合成

　　此代碼是從sentinel2原始文件包中，選出B、G、R、NIR進行波動合成，保存成tif.

python代碼

# -*- coding: utf-8 -*-
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,out,name):
    # 打開柵格數據集
    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 = out+"\\"+name + ".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:\data\實驗結果\哨兵數據下載\重點湖泊矢量_缺2\sentinel2A\6天際湖\S2A_\S2A_MSIL2A_20210622T050651_N0300_R019_T45STT_20210622T081855.SAFE')
    out=r"D:\data\實驗結果\哨兵數據下載"
    name=os.path.basename(SAFE_Path).split('.')[0]
    data_list = glob.glob(SAFE_Path)
    #MSIL2A傳感器
    #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_MSIL1C.xml"    #MTD_MSIL2A   根據傳感器類型選擇

        S2tif(filename,out,name)

        print(data_path + "-----轉tif成功")
    print("----轉換結束----")