【Python&RS】基于GDAL哨兵二号波段合成
学遥感的避免不了使用哨兵数据，毕竟10m的分辨率可以满足大部分的定量分析，同时也是最重要的一点，它免费！！！
之前好像ENVI5.3打不开哨兵数据，易智瑞已给出了解决办法。我想说的是大家能下载L2A级数据就去下，省的麻烦。如果需要的数据只有L1C级，那就使用欧空局发布的Sen2Cor去进行辐射定标和大气校正。
这样就出现了问题，因为Sen2Cor是对单波段进行导出，如果我们一个个导入ENVI去波段组合，那岂不是浪费了很多时间，所以我这里使用Python实现Sentinel-2数据的批量波段组合（PS：Sen2Cor也可以批量大气校正，后面可以出一篇教程）
注意：遥感影像的本质就是多维数组，所以原理并不复杂！看懂这篇文章后同样也可以对其他卫星数据进行波段组合，都是一样的！
一、获取数据的基本信息
这一步主要是为了给我们波段组合后的栅格数据赋值，如投影、分辨率等信息。这段代码在我之前的博文中出现过无数次了，就不多说了。同时对GDAL库的安装之前也有教过，可以自行翻阅：。
这里的基本信息可以从任意一个波段获取，因为他们都是一样的。
ds = gdal.Open(image_path) # 打开数据集dataset
ds_width = ds.RasterXSize # 获取数据宽度
ds_height = ds.RasterYSize # 获取数据高度
ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
ds_prj = ds.GetProjection() # 获取投影信息
二、将基本信息写入新的栅格驱动
这里创建驱动时，输入的参数分别是输出文件名，影像宽度、影像高度、输出的波段数。这里波段数可以依据不同卫星自己修改，同时如果你想加入某些卫星的辅助波段（如云掩膜、水汽等）也可以修改这个参数。
driver = gdal.GetDriverByName('ENVI') # 载入数据驱动，用于存储内存中的数组
ds_result = driver.Create(out_path, ds_width, ds_height, bands=11, eType=gdal.GDT_Int32)
# 创建一个数组，宽高为原始尺寸
ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数
ds_result.SetProjection(ds_prj) # 导入投影信息
# 上述代码用于创建空白数组以及获取投影信息三、遍历文件夹写入波段数据
我这里是写死了波段名称，如果你们想对其他卫星数据进行组合的话需要自己修改，我这里只适用于哨兵二号数据。此外，如果你有多个哨兵数据需要组合，可以在外面再添加一个for循环文件夹即可实现。
band_lists = ["B2.img", "B3.img", "B4.img", "B5.img", "B6.img", "B7.img", "B8.img", "B8A.img", "B9.img", "B11.img", "B12.img"]
i = 1
for band_list in band_lists:
# 遍历波段列表
for image_list in image_lists:
# 遍历文件夹中所有的文件
if band_list == image_list:
# 如果波段列表与文件一致，则将该文件打开
print(str(i) + ".正在融合%s" % band_list)
filepath1 = os.path.join(file_path, image_list)
ds = gdal.Open(filepath1) # 打开数据集dataset
band_data = ds.GetRasterBand(1).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float)
# 以数组的形式读取当前波段
ds_result.GetRasterBand(i).WriteArray(band_data)
i += 1
del ds_result
# 删除内存中的结果，否则结果不会写入图像中四、完整代码
# -*- coding: utf-8 -*-
"""
@Time ： 2023/3/31 15:35
@Auth ： RS迷途小书童
@File ：Band Synthesis.py
@IDE ：PyCharm
@Purpose ：哨兵2号波段组合
"""
import os
import numpy as np
from osgeo import gdal
from osgeo.gdalconst import *
def Layerstack(file_path, out_path, image_path):
"""
:param file_path: 输入波段存放的文件夹
:param out_path: 输出文件路径
:param image_path: 输入任意影像数据，用于定义坐标系等信息
:return:
"""
print("-----开始组合波段-----")
ds = gdal.Open(image_path) # 打开数据集dataset
ds_width = ds.RasterXSize # 获取数据宽度
ds_height = ds.RasterYSize # 获取数据高度
ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
ds_prj = ds.GetProjection() # 获取投影信息
driver = gdal.GetDriverByName('ENVI') # 载入数据驱动，用于存储内存中的数组
ds_result = driver.Create(out_path, ds_width, ds_height, bands=11, eType=gdal.GDT_Int32)
# 创建一个数组，宽高为原始尺寸
ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数
ds_result.SetProjection(ds_prj) # 导入投影信息
image_lists = os.listdir(file_path)
del ds
# 上述代码用于创建空白数组以及获取投影信息
band_lists = ["B2.img", "B3.img", "B4.img", "B5.img", "B6.img", "B7.img", "B8.img", "B8A.img",
"B9.img", "B11.img", "B12.img"]
i = 1
for band_list in band_lists:
# 遍历波段列表
for image_list in image_lists:
# 遍历文件夹中所有的文件
if band_list == image_list:
# 如果波段列表与文件一致，则将该文件打开
print(str(i) + ".正在融合%s" % band_list)
filepath1 = os.path.join(file_path, image_list)
ds = gdal.Open(filepath1) # 打开数据集dataset
band_data = ds.GetRasterBand(1).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float)
# 以数组的形式读取当前波段
ds_result.GetRasterBand(i).WriteArray(band_data)
i += 1
del ds_result
# 删除内存中的结果，否则结果不会写入图像中
print("全部波段合并完成")
if __name__ == "__main__":
file_path1 = r"G:/resampled/"
# 存放单波段的文件夹目录
image_path1 = r"G:/resampled/B2.img"
# 以B2波段的基本信息作为新栅格投影、地理变换的标准
out_path1 = r"G:/GDAL_try/3"
# 输出的文件名，这里是输出ENVI格式，所以不用加后缀名
Layerstack(file_path1, out_path1, image_path1)
# 执行函数
本博文代码是我之前在编写使用GDAL计算NDVI时顺便写出来的，自己已经运行过很多次了。效率相比ENVI手动合成肯定提升了不少，但代码感觉还是有些冗余，能用就行还要啥自行车啊。
如果大家在学习Python或者RS时有什么问题，可以随时留言交流！如果大家对批量处理有兴趣同样可以留言给博主，博主会分享相关代码以供学习！
