Python shapely, incorrectly showing point inside/outside polygon

Question

I am reading polygons from an ESRI Shapefile and I want to test if certain points (coordinates) are inside those polygons. I am using Pyshp to read the shapefiles and Shapely to check if the points are inside. It works fine when each polygon has only one part, but it gives strange results if the polygons have several parts or it they have holes.

The polygons I am testing are irregular and they have tens or hundreds of points each. For simplicity's sake, the examples I am showing involve only squares.

The blue polygon below has two parts. The pink polygon has one hole.

The code first creates a simple shapefile named 'test.shp'. Then I read from the shapefile and use shapely to test if points are inside the polygons.

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

They way you convert the shapes read to shapely objects isn't correct: sr_1.shape.points just returns a flat list of all points instead of a nested list that takes in account the seperate rings of the polygons.

The following code does give the correct results in your tests.

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))