查找单栅格在另一栅格各分割区域内的最大值及小对象移除无效问题

该问题由两个常见使用错误导致：

1. remove_small_objects设置in_place=True时，函数会直接修改输入数组，返回值为None，你同时执行了赋值操作，会直接将raster1覆盖为空值
2. 该函数默认对所有非0值统一识别连通域，若你需要按[1,2,3,4]不同取值分别处理小区域，需要先对栅格做标签标注

修正后的代码如下：

import skimage
import numpy as np
from skimage.measure import label

raster1 = raster1.astype(int)
# 按不同取值生成独立连通域标签
labeled_raster = label(raster1, connectivity=2)
# 移除小于20像素的小区域，不需要加in_place参数，直接接收返回值
raster1 = skimage.morphology.remove_small_objects(labeled_raster, min_size=20, connectivity=2)
# 可选：如果需要恢复原有的1/2/3/4分类值，可做映射处理

不需要手写循环遍历，直接调用现成的统计函数即可，效率远高于自定义循环，两种实现方案如下：

方案1：scipy.ndimage 快速统计 #

适合仅需要获取各分类最大值的场景：

from scipy import ndimage

# 直接输入分类栅格和待统计的数值栅格，指定要统计的分类值即可
category_max = ndimage.maximum(raster2, labels=raster1, index=[1,2,3,4])
# 输出的category_max为数组，顺序对应1、2、3、4类的最大值

如果需要生成和原栅格同尺寸、每个位置填充对应分类最大值的结果栅格，可使用以下代码：

max_raster = ndimage.maximum(raster2, labels=raster1, index=raster1)

方案2：skimage.regionprops 全量统计 #

适合需要同时获取区域面积、平均值、最小值等更多属性的场景：

from skimage.measure import regionprops

props = regionprops(label_image=raster1, intensity_image=raster2)
result = {}
for prop in props:
    # prop.label对应原分类栅格的取值，prop.max_intensity就是对应区域的最大值
    result[prop.label] = prop.max_intensity

内容的提问来源于stack exchange，提问作者Ash