Reverse geocoding from a shapefile of polygons.
The ReverseGeolocShapefile class improves the shapefile.Reader class by
adding a intersection method to find polygons that intersects with a single
point.
I use a simple spatial index that store the bounding box of each polygons.
I save these informations into a .sidx.npy file next to the .dbf/.shp
files.
Works with Python 3.6+ (I'm using f-string).
The intersection method of the ReverseGeolocShapefile class return a list
of row index.
import intersectshape
shp = intersectshape.ReverseGeolocShapefile('/path/of/shapefile.shp')
# You can use all methods of shapefile.Reader object like...
row = shp.record(5)
geo = shp.shape(2)
# Find polygons that intersects a point
x, y = 200, 500 # must be at the same projection system
ids = shp.intersection(x, y)
print(ids)
# return [] or [0, 5, 9, 7] for example.
This code are under the GNU GPL v3 licence.
Data of French communal polygons (Geofla shapefile located in the data folder) are under Etalab Open License.