Python rasterio.merge遥感影像融合

最近在尝试用Python做遥感影像融合，用到了rasterio.merge.merge（简称merge）功能，可能还有gdal和arcpy，但是还没有用过。

问题：

merge虽然是非常好用的图像融合工具，但是它只提供四种内置融合方式，first,last,min和max，即

MERGE_METHODS = {
    'first': copy_first,
    'last': copy_last,
    'min': copy_min,
    'max': copy_max
}

因为我想用均值融合，所以不能满足需求。

解决方法：

以max最大值融合为例，

def copy_max(old_data, new_data, old_nodata, new_nodata, **kwargs):
	'''找到old_data和new_data重叠的部分'''
    mask = np.empty_like(old_data, dtype='bool')
    np.logical_or(old_nodata, new_nodata, out=mask)
    np.logical_not(mask, out=mask)
	'''对重叠部分计算最大值'''
    np.maximum(old_data, new_data, out=old_data, where=mask)
	'''将new_data有数据而old_data没有数据的像元，把new_data对应值复制到old_data'''
    np.logical_not(new_nodata, out=mask)
    np.logical_and(old_data, mask, out=mask)
    np.copyto(old_data, new_data, where=mask)

这里用到了np.maximum函数，使用方法可以参考https://numpy.org/doc/stable/reference/generated/numpy.maximum.html
仿照上面代码，改进的平均值融合为：

def copy_mean(merged_data, new_data, merged_mask, new_mask, **kwargs):
    mask = np.empty_like(merged_mask, dtype="bool")
    
    '''因为我的数据中同时存在0值和nodata，其实0也是nodata，因此这里把0也放到mask中，用不到可以删除'''
    mask_merged_zero = merged_data==0.0
    mask_new_zero = new_data==0.0
    np.logical_or(mask_new_zero, new_mask, out=new_mask) 
    np.logical_or(mask_merged_zero, merged_mask, out=merged_mask) 
    
    '''提取重叠区域'''
    np.logical_or(merged_mask, new_mask, out=mask)
    np.logical_not(mask, out=mask)
    
    '''求均值'''
    np.add(merged_data,new_data,out=merged_data,where=mask)
    np.true_divide(merged_data,2,out=merged_data,where=mask)
    
    np.logical_not(new_mask, out=mask)
    np.logical_and(merged_mask, mask, out=mask)
    np.copyto(merged_data, new_data, where=mask, casting="unsafe")

np.add和np.true_divide和np.maximum的用法类似，所以没有用np.mean或者np.average

rasterio.merge的用法可以参考
https://automating-gis-processes.github.io/CSC18/lessons/L6/raster-mosaic.html

如果有更好的方法欢迎分享~！

融合效果：