P i p
Point in polygon. There are two basic methods with numerous variants. One or the other in some form forms the foundation of spatial queries.
I have opted to present the winding number approach since it is the one I understood the least. Pictures of polygon boundaries flipping around the point wasn't cutting it for me, nor was the fact that most implementations checked one point at a time to see if it fit into one polygon.
So much for the fantasy of the points automagically finding their polygon home in an instance. I started of with cursors and point geometry and polygon objects and soon realized that there has to be an alternate way. Dr Google was useless. The same references kept coming up. A useful history, but the academic literature provided nothing more than variants on a theme.
I was hoping that you might be able to check a whole bunch of points at once to see if they were in a polygon. Nadda, Zilch. Except for one glimmer of hope found buried on GitHub. I had my inspiration and the wheel would once again be reinvented.
Some sample shapes appear to the right. Two polygons were created in ArcGIS Pro. The upper is their array representation as an svg graphic and the lower two are arcpy polygons similarly represented.
The winding number algorithm in its simplest form is show below.
It is pretty simple. Parse the points into those that are obviously beyond the extent of the bounds of the polygon, then get fussy and check every side of a polygon edge to see if it is enclosed in the bounds of the polygon.
The sided-ness of a point relative to a polygon edge is handled by _is_right_side.
Practically, the y value of a point is inspected to see if it is between two points on a segment (line 18-19). If it is, then the sided-ness check is made and if the point's y is on the right side, then the w counter is incremented. Failing that test, it is decrease. If neither condition applies, nothing happens. In the end any nonzero locations are within the polygon, the rest outside. I should note here, that some implementations provide nonzero initial values, so that the algorithm looks for zero values instead. The nice thing about the winding number vs the crossing number approach is that the former only involves only side checks and addition and multiplication, the latter has a division requiring a denominator check for 0.
Of course there are the perennial questions of "what if it is on the boundary?"... "what about holes?"... "what if I walk counterclockwise?" Interesting questions, but none really of interest to me. I couldn't get past the fact that apparently the inclusion test was supposed to be done one point at a time.
So it is heartwarming to see that all the points are returning the same points in those 2 highly irregular polygons. My preference is for the winding number since no segment intersection is required, hence no division. Also, winding number is much easier to vectorize (see np_wn below) since the last thing I want to do is test each point against each line segment.
|Winding number for point in polygon inclusion tests.|
The numpy incarnation of winding number simply compares all the points to all sides rather than one point to one side a time.
For those that are interested in time testing... read on. To the rest of you... until next time.
Creating arcpy Points and Polygons from x,y coordinates
Making the points and polygons from numpy array data is fairly simple. Arcpy's object model forces you to create point objects, then throw them into an Array object, from which you can form a Polygon object. Sadly, there is no shortcut to this even if you don't ever need the properties and methods associated with the objects used in the construction path from numbers to polygon.
Simple tests for inclusion
Not bad... 100 points against 73 and 22 edges respectively.
A valiant effort, hardly a sip of coffee. I wonder how long it will take for a much larger dataset.
Where my interests lie. Can this whole thing be vectorized? A lot faster it appears.
No slouch, especially when you have to consider the extra baggage a multipoint carries.