2

I have a shapefile with a large grid of rectangular polygons (approximately 6M polygons) and I would like to pull an attribute value from a polygon in this set based on an arbitrary point (x,y) somewhere within the extents of these polygons; e.g., for the point (xi,yi) the 'value' would be in some polygon Pn.

I am uncertain where to start - I have another project where I used the rasterio.sample() method to get a raster value at a specific point so, I feel like some similar approach would be applicable to look through the polygon space...

Any pointers?

1

2 Answers 2

3

Create a point df using gpd.points_from_xy and then spatial join. It is very fast:

import geopandas as gpd
from sqlalchemy import create_engine as ce
con = ce("postgresql://postgres@localhost:5432/data")
sql = "select * from test.grid_2e6"

grid = gpd.read_postgis(sql=sql, con=con)
print(grid.shape)
#(1981926, 2)

print(grid.head(2))
#    id                                               geom
# 0   1  POLYGON ((244550.000 6110830.000, 244550.000 6...
# 1   2  POLYGON ((244550.000 6111560.000, 244550.000 6...

y, x = 6922071.8,470881.6
point = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=[x], y=[y]), crs=grid.crs)
print(point)
#                          geometry
# 0  POINT (470881.600 6922071.800)

sj = gpd.sjoin(left_df=point, right_df=grid, how="left", rsuffix="grid_")
print(sj)
#                          geometry  index_grid_      id
# 0  POINT (470881.600 6922071.800)       663891  663892
1
  • 1
    spatial join turned out to be the approach I needed! Thank you. (the Geopackage suggestion below was also very helpful for the initial loading of data)
    – Jimbobauer
    Dec 9, 2023 at 15:42
3

Depending on what you need exactly, geopandas also supports immediate filtering on the file. If you only need the row(s) for one or some points (the mask can also be a GeoDataFrame with multiple points), this will be faster than reading the entire file.

When using a file type with a proper spatial index like e.g. Geopackage this will be even faster (100x). Timings using a file of ~380 MB with ~500.000 (more complex) polygons on a Windows PC.

from datetime import datetime
import geopandas as gpd
from shapely import Point

start = datetime.now()
df = gpd.read_file("C:/Temp/prc2023/prc2023.shp", mask=Point(150000, 185000))
print(f"read .shp, found {len(df)} rows, took {datetime.now() - start}")

start = datetime.now()
df = gpd.read_file("C:/Temp/prc2023/prc2023.gpkg", mask=Point(150000, 185000))
print(f"read .gpkg, found {len(df)} rows, took {datetime.now() - start}")

Output:

read .shp, found 1 rows, took 0:00:04.417004
read .gpkg, found 1 rows, took 0:00:00.057872
1
  • Thank you - this was very helpful
    – Jimbobauer
    Dec 9, 2023 at 15:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.