Mask R-CNN农田预测多边形芯片拼接问题及优化方案问询

问题背景 #

我们团队近期采用Mask-R-CNN模型，基于卫星影像预测农田空间范围。流程为将卫星影像分割为512×512像素芯片，逐芯片预测农田后拼接多边形，发现模型未对芯片边缘进行预测，导致芯片边缘存在数据缺失：

现有拼接修复方案 #

提出以下拼接修复流程：

• 利用原始未裁剪影像的角点坐标与空间分辨率，生成模拟512×512芯片的fishnet多边形；
• 按每个fishnet对shapefile进行子集提取；
• 检查子集内每个多边形的顶点坐标，找出距离边界2米内且重复出现的东向、北向坐标（南北边界对应北向坐标，东西边界对应东向坐标）——因多边形在2米范围被截断且呈直线，会存在重复坐标；
• 将包含上述坐标的多边形顶点，替换为对应fishnet的边界坐标；
• 对所有fishnet内的多边形完成处理后，使用gpd.unary_union将数据集合并为单个多面多边形，再进行拆分（explode）。

处理后的效果如下：

现有方案存在的问题 #

该方案多数情况有效，但存在部分多边形仅部分合并、少数多边形完全未合并的问题：

当前实现代码 #

import os
import pandas as pd
import geopandas as gpd
import rasterio
from shapely.geometry import Polygon
from collections import Counter
from pathlib import Path


def stitchshp(img, polys, endpath):
    img = rasterio.open(img)
    polys = gpd.read_file(polys)
   #Making Fishnets from the image
    firsteast = img.bounds.left
    firstnorth = img.bounds.top
    easting = [img.bounds.left]
    northing = [img.bounds.top]

    def funeast(image):
        while image < img.bounds.right:
            image = image +(0.5*512)
            easting.append(image)
            return funeast(image)
        
    def funnorth(image):
        while image > img.bounds.bottom:
            image = image -(0.5*512)
            northing.append(image)
            return funnorth(image)


    funeast(firsteast)
    funnorth(firstnorth)

    geom = []

    for i in range(len(northing)):
        for j in range(len(easting)):
            coords = [(easting[j], northing[i]), (easting[j]+(0.5*512), northing[i]), (easting[j]+(0.5*512), northing[i]-(0.5*512)), (easting[j], northing[i]-(0.5*512))]
            geom.append(Polygon(coords))

    table = gpd.GeoDataFrame({"geometry": geom}, crs ={'init' :'epsg:32644'} )


    #using fishnets to subset manipulate, and stitch polygons
    allpolyslist = []
    for fish in range(len(table.geometry)):
        net = table.geometry[fish]
        east = max([i[0] for i in [*net.exterior.coords]])
        west = min([i[0] for i in [*net.exterior.coords]])
        north = max([i[1] for i in [*net.exterior.coords]])
        south = min([i[1] for i in [*net.exterior.coords]])
        subpolys = polys[polys.intersects(net)]
        net = gpd.GeoSeries(net)

        maxeast = []
        mineast = []
        maxnorth = []
        minnorth = []
        for sequ in [[*i.exterior.coords] for i in subpolys.geometry]:
            maxeast.append(max(i[0] for i in sequ))
            mineast.append(min(i[0] for i in sequ))
            maxnorth.append(max(i[1] for i in sequ))
            minnorth.append(min(i[1] for i in sequ))
        
        eastdict = Counter(maxeast)
        eastchange = [key for key, value in eastdict.items() if key >east-2 and key <east+2 and value>1]
        westdict = Counter(mineast)
        westchange = [key for key, value in westdict.items() if key <west+2 and key >west-2 and value>1]
        northdict = Counter(maxnorth)
        northchange = [key for key, value in northdict.items() if key >north-2 and key <north+2 and value>1]
        southdict = Counter(minnorth)
        southchange = [key for key, value in southdict.items() if key <south+2 and key >south-2 and value>1]

        polylists = [[*i.exterior.coords] for i  in subpolys.geometry]

        properpolys = []
        for poly in polylists:
            if len(poly)>6:
                properpolys.append(poly)


        for poly in properpolys:
            for j in range(len(poly)):
                if poly[j][0] in eastchange:
                    poly[j] = (east, poly[j][1])
                if poly[j][0] in westchange:
                    poly[j] = (west, poly[j][1])
                if poly[j][1] in northchange:
                    poly[j] = (poly[j][0], north)
                if poly[j][1] in southchange:
                    poly[j] = (poly[j][0], south)
        

        
        allpolyslist = allpolyslist + properpolys
    print("ok, iterated thru!")
    #gpd.GeoSeries([Polygon(i) for i in allpolyslist]).plot()   

    allpolylist = [Polygon(i) for i in allpolyslist]
    print("converted to polygons!")
    allpolydict = gpd.GeoDataFrame({"geometry": allpolylist}, crs = {'init' :'epsg:32644'})
    merge = allpolydict.unary_union
    merge = gpd.GeoSeries(merge)
    merge = merge.explode()
    print("exploded polys!")
    united = gpd.GeoDataFrame(geometry = merge)
    print("Made final df!")
    united.to_file(endpath, crs = allpolydict.crs)

优化思路探讨 #

针对部分多边形未完全合并或未合并的问题，可从以下方向优化：

1. 修正Fishnet生成逻辑 #

当前递归函数funeast和funnorth存在逻辑错误：递归调用在while循环内直接return，导致仅生成一个步长的坐标，无法遍历整个影像范围。需修改为循环遍历而非递归，同时结合影像分辨率生成正确的网格步长：

# 替换原递归函数，改用循环生成完整坐标
def generate_eastings(start, end, step):
    eastings = [start]
    current = start
    while current + step < end:
        current += step
        eastings.append(current)
    eastings.append(end)  # 确保包含右边界
    return eastings

def generate_northings(start, end, step):
    northings = [start]
    current = start
    while current - step > end:
        current -= step
        northings.append(current)
    northings.append(end)  # 确保包含底边界
    return northings

# 基于影像分辨率计算正确步长（512像素 × 单个像素尺寸）
pixel_size = img.res[0]
step = 512 * pixel_size
easting = generate_eastings(firsteast, img.bounds.right, step)
northing = generate_northings(firstnorth, img.bounds.bottom, step)

2. 改进边界顶点替换逻辑 #

• 移除重复次数判断：直接判断顶点是否在fishnet边界的2米缓冲范围内，无需依赖重复次数，避免遗漏仅出现一次的边界截断顶点；
• 批量处理边：如果某条边的两个端点都靠近同一fishnet边界，直接将整条边替换为fishnet的对应边界线段，提升修复精度。

3. 合并前的预处理 #

• 拓扑修复：对所有多边形执行shapely.make_valid()，修复无效几何结构；
• 碎片过滤：移除面积过小的多边形碎片，避免干扰合并过程；
• 跨网格预合并：处理单个fishnet时，同步提取相邻网格的多边形，提前合并跨网格的部分。

4. 合并后的后处理 #

• 缝隙填补：对合并后的多边形执行微小缓冲操作（如先缓冲0.1米，再缓冲-0.1米），自动填补小缝隙：

# 缓冲修复缝隙
united['geometry'] = united['geometry'].buffer(0.1).buffer(-0.1)

• 二次拓扑检查：对拆分后的多边形再次验证有效性，修复可能存在的重叠或断裂。

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