Suppose we wish to understand if a given X/Y pair falls within a set of non-intersecting arbitrary polygons. For example, we are trying to see if a given lat,long falls within a Michigan minor civil division. There are 1500+ such divisions each one of which could be a polygon with hundreds or thousands of points. The following is a practical solution that is fast enough to return to a user.
- Break up Michigan into
$N$ box oriented tiles. We can easily calculate ifx, ywill only fall within a tile based on the value(x % WIDTH, y % HEIGHT). Select$N$ such that most (e.g. ~60%?) tiles are fully contained within a polygon, and the vast majority (e.g. 95%) are in no more than 4 polygons. I imagine this would happen easy at$N$ less than 100,000. - Preprocess the polygon shape files using Shapely to determine a mapping from the tile to the intersection of the tile with it's polygons. This will be a one-to-many mapping. Store this into a json blob.
- At run time, you can easily determine which tile a point is in. Using the mapping computed in #2, we can run a fairly a straightforward point in polygon typescript algorithm for the few polygon intersections which map be relevant. Once we determine the polygon intersection, we can then return the name of the polygon.
Steps #1 and #2 are only for performance. You can obviously brute force #3.