Wednesday 28 February 2018

python - Does shapely within function identify inner holes?


I have a series of shapefiles that include both polygons and multipolygons. I am using ogr to extract them in wkb, and them I am using shapely to check if a point fall inside a polygon (or a multipolygon). However, examining some of the polygons in wkb I realized that many of them have inner rings (i.e. I have a set of coordinates describing the first, outer polygon, and a set of coordinates describing the inner polygon/s). Does the shapely within function take care of this, I mean, does it return True only if a point falls in the area delimited by outer and inner polygons? And are the inner polygons always supposed to be holes?


I am sorry, perhaps these are really dumb questions, but I am completely new to this kind of analyses and, despite my efforts, I could not find some clear information about these doubts. Thanks in advance for your help.



Answer



'Yes' is the short answer to your question. Interior rings in a single polygon are always holes. The exterior ring comes first followed by a list of interior rings (holes). A polygon must also be non-self-intersecting. Multipolygons are merely a list of polygons bundled together, but the structure of each consituent polygon in the list follows the same structure as for a single 'stand-alone' polygon. So, The interior rings of the constituent polygons of a multipolygon are also holes. Shapely honours this, so if your point falls in the hole, the tests .within and .contains will both return False.


You can test this very easily as follows:


from shapely.geometry import Polygon
from shapely.geometry import Point


# construct a valid polygon (i.e. not self-intersecting) with a hole
ext = [(0, 0), (0, 2), (2, 2), (2, 0), (0, 0)]
int = [(1, 1), (1, 1.5), (1.5, 1.5), (1.5, 1)]
myPoly = Polygon(ext,[int])

# construct a point that is within the bounds of the exterior ring but inside the hole
myPoint = Point(1.25, 1.25)
print(myPoint.within(myPoly)) # returns 'False'

# construct a point that is inside the polygon (avoiding the hole)

myOtherPoint = Point(0.5, 0.5)
print(myOtherPoint.within(myPoly)) # returns 'True'

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...