Python 形状正确，错误地显示多边形内部/外部的点

Question

我正在从 ESRI Shapefile 读取多边形，我想测试某些点（坐标）是否在这些多边形内。我使用 Pyshp 读取 shapefile，并使用 Shapely 检查点是否在内部。 当每个多边形只有一个部分时它工作得很好，但如果多边形有多个部分或者它们有孔，它会给出奇怪的结果。

我正在测试的多边形是不规则的，每个多边形都有数十或数百个点。为了简单起见，我展示的示例仅涉及正方形。

下面的蓝色多边形有两部分。粉红色的多边形有一个洞。

代码首先创建一个名为“test.shp”的简单形状文件。然后我从 shapefile 中读取并使用 shapely 来测试点是否在多边形内部。

from shapely import Point, Polygon
import shapefile


w = shapefile.Writer('test', shapeType=5)
w.field('ID', 'C', size=4)
w.field('Polygon', 'C', size=15)
w.field('Descript', 'C', size=40)

# Add Polygon A, the blue one, with 2 parts
w.poly([
       [[0, 0], [0, 5], [5, 5], [5, 0], [0, 0]],            # poly 1, clockwise
       [[10, 10], [10, 20], [20, 20], [20, 10], [10, 10]]   # poly 2, clockwise
       ])
w.record(w.shpNum, 'Polygon A', 'Two separated squares (Aa and Ab)')


# Add Polygon B, the pink one, with one hole
w.poly([
       [[40, 0], [40, 20], [60, 20], [60, 0], [40, 0]],     # poly 1, clockwise
       [[50, 10], [55, 10], [55, 15], [50, 15], [50, 10]]   # hole 1, counterclockwise
       ])
w.record(w.shpNum, 'Polygon B', 'Square with squared hole')

w.close()


# List of 4 points, the green ones, that I want to see if are inside Polygon A
test_points_1 = (
    Point(3, 3),    # Inside Polygon A, part Aa
    Point(6, 6),    # Outside Polygon A
    Point(15, 15),  # Inside Polygon A, part Ab
    Point(25, 25)   # Outside Polygon A
    )

# List of 4 points, the black ones, that I want to see if are inside Polygon B
test_points_2 = (
    Point(45, 5),   # Inside Polygon B
    Point(52, 12),  # Outside Polygon B (inside hole)
    Point(58, 14),  # Inside Polygon B
    Point(65, 8)    # Outside Polygon B
    )

#
# Open the shapefile, reads the points create a Polygon object with shapely and check
#

sfr = shapefile.Reader("test", encoding = "iso8859-1")

# Test if points from test_points_1 are inside Polygon A
sr_1 = sfr.shapeRecords()[0]
points = sr_1.shape.points
poly = Polygon(points)
for pt in test_points_1:
    if pt.within(poly):
        print("Point ({0},{1}) is inside Polygon A".format(pt.x, pt.y))
    else:
        print("Point ({0},{1}) is outside Polygon A".format(pt.x, pt.y))


# Result:
# Point (3.0,3.0) is outside Polygon A     <-- Incorrect
# Point (6.0,6.0) is outside Polygon A     <-- Correct
# Point (15.0,15.0) is inside Polygon A    <-- Correct
# Point (25.0,25.0) is outside Polygon A   <-- Correct



# Test if points from test_points_2 are inside Polygon B
sr_2 = sfr.shapeRecords()[1]
points = sr_2.shape.points
poly = Polygon(points)
for pt in test_points_2:
    if pt.within(poly):
        print("Point ({0},{1}) is inside Polygon B".format(pt.x, pt.y))
    else:
        print("Point ({0},{1}) is outside Polygon B".format(pt.x, pt.y))

# Result
# Point (45.0,5.0) is outside Polygon B    <-- Incorrect
# Point (52.0,12.0) is outside Polygon B   <-- Correct
# Point (58.0,14.0) is inside Polygon B    <-- Correct
# Point (65.0,8.0) is outside Polygon B    <-- Correct




# However, if I don't check against the whole Polygon B, but I test against each of its parts:
sr_2 = sfr.shapeRecords()[1]
points = sr_2.shape.points
parts = sr_2.shape.parts

# Create one Polygon object with a subset of the points
sub_poly_1 = Polygon(points[0:parts[1]])
for pt in test_points_2:
    if pt.within(sub_poly_1):
        print("Point ({0},{1}) is inside Polygon B, part 0".format(pt.x, pt.y))
    else:
        print("Point ({0},{1}) is outside Polygon B, part 0".format(pt.x, pt.y))

# Point (45.0,5.0) is inside Polygon B, part 0    <-- Correct
# Point (52.0,12.0) is inside Polygon B, part 0   <-- Correct
# Point (58.0,14.0) is inside Polygon B, part 0   <-- Correct
# Point (65.0,8.0) is outside Polygon B, part 0   <-- Correct


# Create another Polygon object with the rest of the points
sub_poly_2 = Polygon(points[parts[1]:])
for pt in test_points_2:
    if pt.within(sub_poly_2):
        print("Point ({0},{1}) is inside Polygon B, part 1".format(pt.x, pt.y))
    else:
        print("Point ({0},{1}) is outside Polygon B, part 1".format(pt.x, pt.y))

# Point (45.0,5.0) is outside Polygon B, part 1    <-- Correct
# Point (52.0,12.0) is inside Polygon B, part 1    <-- Correct
# Point (58.0,14.0) is outside Polygon B, part 1   <-- Correct
# Point (65.0,8.0) is outside Polygon B, part 1    <-- Correct

Answer 1

将读取的形状转换为形状良好的对象的方式是不正确的：

sr_1.shape.points

仅返回所有点的平面列表，而不是考虑多边形的单独环的嵌套列表。

以下代码确实在您的测试中给出了正确的结果。

from shapely import Point
from shapely.geometry import shape
import shapefile

w = shapefile.Writer("test", shapeType=5)
w.field("ID", "C", size=4)
w.field("Polygon", "C", size=15)
w.field("Descript", "C", size=40)

# Add Polygon A, the blue one, with 2 parts
w.poly(
    [
        [[0, 0], [0, 5], [5, 5], [5, 0], [0, 0]],  # poly 1, clockwise
        [[10, 10], [10, 20], [20, 20], [20, 10], [10, 10]],  # poly 2, clockwise
    ]
)
w.record(w.shpNum, "Polygon A", "Two separated squares (Aa and Ab)")


# Add Polygon B, the pink one, with one hole
w.poly(
    [
        [[40, 0], [40, 20], [60, 20], [60, 0], [40, 0]],  # poly 1, clockwise
        [[50, 10], [55, 10], [55, 15], [50, 15], [50, 10]],  # hole 1, counterclockwise
    ]
)
w.record(w.shpNum, "Polygon B", "Square with squared hole")

w.close()


# List of 4 points, the green ones, that I want to see if are inside Polygon A
test_points_1 = (
    Point(3, 3),  # Inside Polygon A, part Aa
    Point(6, 6),  # Outside Polygon A
    Point(15, 15),  # Inside Polygon A, part Ab
    Point(25, 25),  # Outside Polygon A
)

# List of 4 points, the black ones, that I want to see if are inside Polygon B
test_points_2 = (
    Point(45, 5),  # Inside Polygon B
    Point(52, 12),  # Outside Polygon B (inside hole)
    Point(58, 14),  # Inside Polygon B
    Point(65, 8),  # Outside Polygon B
)

#
# Open the shapefile, reads the points create a Polygon object with shapely and check
#

sfr = shapefile.Reader("test", encoding="iso8859-1")

# Test if points from test_points_1 are inside Polygon A
sr_1 = sfr.shapeRecords()[0]
poly = shape(sr_1.shape)
for pt in test_points_1:
    if pt.within(poly):
        print("Point ({0},{1}) is inside Polygon A".format(pt.x, pt.y))
    else:
        print("Point ({0},{1}) is outside Polygon A".format(pt.x, pt.y))

# Test if points from test_points_2 are inside Polygon B
sr_2 = sfr.shapeRecords()[1]
poly = shape(sr_2.shape)
for pt in test_points_2:
    if pt.within(poly):
        print("Point ({0},{1}) is inside Polygon B".format(pt.x, pt.y))
    else:
        print("Point ({0},{1}) is outside Polygon B".format(pt.x, pt.y))