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 if
x, y
will 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 map