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

遥感图像叠加，网格图像叠加，不同尺寸和波段数
背景
高分4号卫星配备了两个不同的传感器PMS数据包含5个波段，IRS只有一个中红外波段，它们的尺寸和数量不同，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

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