Python中多边形/Shapefile内生成网格点：无法判断点是否在多边形内

我的目标是在矩形（也可为不规则多边形）内生成带计算间距的点坐标。以下是我使用的代码，但它无法判断点是否在多边形内——即使生成的点在多边形内部，输出的点数也为0。请问我哪里出错了？

import numpy as np
import geopandas as gpd
from shapely.geometry import Point


df = gpd.read_file(r'C:\Users\User\Desktop\PyTrials\Rectangle.shp')

def grid(row, pointcount):


    finished = 0
    maxiter = 20 
    iterations = 0
    #prep_polygon =prep(row)
    #pointcount=2
    while finished != 1 or iterations<maxiter: #Until there are n points inside the polygon or maxiter is reached
        iterations+=1
        xmin, ymin, xmax, ymax = row.geometry.bounds #Find polygon bounds


        #Find suitable spacing. The 50 buffer is to prevent points at polygon border.
        spacing=(((row.geometry.buffer(50).area/pointcount)**0.5)*0.99)*(1.3-np.random.random()*0.6)
       
        nspacingx=np.ceil((xmax-xmin)/spacing) #Number of points in x
        nspacingy=np.ceil((ymax-ymin)/spacing) #and y

        randomstart=[xmin+spacing*np.random.random(),ymin+spacing*np.random.random()] #A random start coordinate

        xlist=[randomstart[0]+(x*spacing) for x in range(int(nspacingx)+1)] #Create a list of x coordinates
        ylist=[randomstart[1]+(y*spacing) for y in range(int(nspacingy)+1)] #And y

        #Using meshgrid, list all points combinations [Point(x1,y1), Point(x1,y2), ... and check if they are inside the polygon NOT WORKING
        points = [Point(coords) for coords in np.array(np.meshgrid(xlist, ylist)).T.reshape(-1, 2).tolist() if Point(coords).intersects(row.geometry)] #This step is when it is unable to evaluate if #points are inside the polygon

        numberpoints=(len(points))
        print (numberpoints) #when I print this, it gives me 0 even though #the generated points are inside the polygon

        if len(points)==pointcount:
            #print(sum(points))
            return points
            finished =1

df['gridlist'] = df.apply(grid, args=(10,), axis=1) #Apply the function with 10 m spacing
df = df.explode('gridlist') #Explode so each point in gridlist become a row
df['geometry'] = df['gridlist']
df = df.drop('gridlist', axis=1)

df.to_file(r'C:\Users\User\Desktop\PyTrials\PointsRectangle.shp')

问题分析与修复 #

1. 点与多边形的包含判断逻辑错误 #

你使用Point(coords).intersects(row.geometry)判断点是否在多边形内，这个方法会包含点与边界接触的情况，且在带投影的坐标体系下容易因精度问题返回错误结果。应该改用**row.geometry.contains(Point(coords))或Point(coords).within(row.geometry)**，这两个方法能更准确地判断点是否严格在多边形内部（若需要包含边界，可替换为covers方法）。

2. while循环条件逻辑错误 #

原代码while finished != 1 or iterations<maxiter的逻辑是“未完成或未到最大迭代次数就继续循环”，会导致循环无法正常终止。正确逻辑应为while finished != 1 and iterations < maxiter——“未完成且未到最大迭代次数时才继续循环”。另外finished =1写在return之后，永远不会被执行，可直接删除。

3. Buffer过大导致Spacing计算错误 #

用row.geometry.buffer(50).area计算面积会让多边形面积大幅膨胀，进而算出的spacing远大于实际需求，生成的点几乎全部落在原始多边形外部，最终导致点数为0。若要避免点落在边界上，应在生成点后过滤边界点，而非先给多边形扩边计算面积。

修正后的代码 #

import numpy as np
import geopandas as gpd
from shapely.geometry import Point


df = gpd.read_file(r'C:\Users\User\Desktop\PyTrials\Rectangle.shp')

def grid(row, pointcount):
    maxiter = 20 
    iterations = 0
    polygon = row.geometry
    xmin, ymin, xmax, ymax = polygon.bounds

    while iterations < maxiter:
        iterations += 1
        # 基于原始多边形面积计算间距，移除不必要的buffer
        spacing = (((polygon.area / pointcount)**0.5) * 0.99) * (1.3 - np.random.random()*0.6)
        # 计算x、y方向的点数量
        nspacingx = np.ceil((xmax - xmin) / spacing)
        nspacingy = np.ceil((ymax - ymin) / spacing)
        # 生成随机起始点
        randomstart = [xmin + spacing*np.random.random(), ymin + spacing*np.random.random()]
        # 生成网格坐标列表
        xlist = [randomstart[0] + x*spacing for x in range(int(nspacingx)+1)]
        ylist = [randomstart[1] + y*spacing for y in range(int(nspacingy)+1)]
        # 生成所有网格点并筛选在多边形内部的点
        coords_list = np.array(np.meshgrid(xlist, ylist)).T.reshape(-1, 2).tolist()
        points = [Point(coords) for coords in coords_list if polygon.contains(Point(coords))]
        
        numberpoints = len(points)
        print(f"Iteration {iterations}: {numberpoints} points found")
        
        # 若找到足够数量的点，直接返回（可切片取刚好pointcount个）
        if len(points) >= pointcount:
            return points[:pointcount]
    # 达到最大迭代次数仍未找到，返回当前结果
    return points

# 应用函数生成10个点
df['gridlist'] = df.apply(grid, args=(10,), axis=1)
df = df.explode('gridlist')
df['geometry'] = df['gridlist']
df = df.drop('gridlist', axis=1)

df.to_file(r'C:\Users\User\Desktop\PyTrials\PointsRectangle.shp')

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