实现不同尺寸和不同波段数的图像叠加(layerstaking)——按地理位置对准GDAL,Numpy

遥感影像叠加，栅格图像叠加，不同尺寸和不同波段数
背景
高分四号卫星搭载两个不同的传感器，获取的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

完成后的波段组合显示效果，位置精确，达到想要的效果