遥感影像叠加，栅格图像叠加，不同尺寸和不同波段数
背景
高分四号卫星搭载两个不同的传感器，获取的PMS数据包含5个波段，IRS只有1个中红外波段，他们的尺寸和波段数量是不同的，IRS数据包含在PMS数据内部。叠加效果如下图：
思路
考虑到它们的左上角并不是对齐的，所以想要实现波段叠加，不能直接将它们的数据矩阵叠加，否则必然错位。而是要找到IRS左上角在PMS中对应的位置（x，y），然后再写到同PMS一样大的array中。
代码实现
def layerstaking(PMS,IRS,merge):
"实现不同尺寸和不同波段数的影像叠加——按地理位置对准"
"针对高分4号PMS和IRS数据"
"输入影像1（上层），输入影像2（下层），输出文件名"
# 读取
IRS_ds = gdal.Open(IRS)
IRS_geotrans = IRS_ds.GetGeoTransform()
IRS_lefttop_x,IRS_lefttop_y = IRS_geotrans[0],IRS_geotrans[3]
IRS_cols = IRS_ds.RasterXSize
IRS_rows = IRS_ds.RasterYSize
IRS_array = IRS_ds.GetRasterBand(1).ReadAsArray()
PMS_ds = gdal.Open(PMS)
PMS_cols = PMS_ds.RasterXSize
PMS_rows =PMS_ds.RasterYSize
PMS_geotrans = PMS_ds.GetGeoTransform()
PMS_lefttop_x, PMS_lefttop_y ,PMS_x_res, PMS_y_res = PMS_geotrans[0], PMS_geotrans[3], PMS_geotrans[1], PMS_geotrans[5]
# 计算IRS在PMS中的相对位置，即行列号
x_situ = int(abs((PMS_lefttop_x-IRS_lefttop_x)/PMS_x_res))
y_situ = int(abs((PMS_lefttop_y-IRS_lefttop_y)/PMS_y_res))
# 创建与PMS同样尺寸的图层来存储IRS
new_band = np.full([PMS_rows, PMS_cols], -999.)
# 将IRS数据放入对应位置
new_band[y_situ:(y_situ+IRS_rows), x_situ:(x_situ+IRS_cols)] = IRS_array
# 数据写出
Driver = PMS_ds.GetDriver()
nbands = PMS_ds.RasterCount + 1
out_name = merge
outDataset = Driver.Create(out_name, PMS_cols, PMS_rows, nbands, gdal.GDT_Float32)
for m in range(nbands - 1): # PMS
ReadBand = PMS_ds.GetRasterBand(m + 1)
outband = outDataset.GetRasterBand(m + 1)
outband.SetNoDataValue(-999)
Image = ReadBand.ReadAsArray()
outband.WriteArray(Image)
outDataset.GetRasterBand(6).WriteArray(new_band) # IRS
outDataset.FlushCache()
outDataset.SetGeoTransform(PMS_ds.GetGeoTransform())
outDataset.SetProjection(PMS_ds.GetProjection())
print("end")
del IRS_ds, PMS_ds, outDataset
工具
基于python，GDAL, Numpy
完成后的波段组合显示效果，位置精确，达到想要的效果
