Spatial Join Highest Intersection

find intersection of largest area joined

Import

Generate Test Data

Load a sample admin boundary file using geopandas

admin_bounds_gdf = gpd.read_file("../data/geoboundary.geojson")

admin_bounds_gdf.head(3)

	shapeName	shapeISO	shapeID	shapeGroup	shapeType	geometry
0	Abra	None	PHL-ADM2-3_0_0-B1	PHL	ADM2	MULTIPOLYGON (((120.96795 17.95706, 120.97803 ...
1	Agusan del Norte	None	PHL-ADM2-3_0_0-B2	PHL	ADM2	MULTIPOLYGON (((125.57724 9.45679, 125.59687 9...
2	Agusan del Sur	None	PHL-ADM2-3_0_0-B3	PHL	ADM2	MULTIPOLYGON (((125.91087 8.85625, 125.91461 8...

admin_bounds_gdf.dtypes

shapeName       object
shapeISO        object
shapeID         object
shapeGroup      object
shapeType       object
geometry      geometry
dtype: object

# admin_bounds_gdf.explore()

Generate a sample grid

grid_generator5k = grids.SquareGridGenerator(50_000)  # 5 km x 5 km square cells

grid_gdf = grid_generator5k.generate_grid(admin_bounds_gdf)

CPU times: user 2.38 s, sys: 112 ms, total: 2.49 s
Wall time: 2.65 s

grid_gdf.plot();

ax = admin_bounds_gdf.plot(facecolor="grey", edgecolor="grey", alpha=0.2)
ax = grid_gdf.plot(ax=ax, facecolor="none", edgecolor="green")

grid_gdf.describe()

	x	y
count	323.000000	323.000000
mean	11.773994	15.018576
std	4.861594	8.791210
min	0.000000	0.000000
25%	8.000000	8.000000
50%	11.000000	14.000000
75%	16.000000	21.000000
max	21.000000	37.000000

grid_gdf.head(3)

	x	y	geometry
0	3	3	POLYGON ((118.27581 5.92871, 118.72497 5.92871...
1	2	3	POLYGON ((117.82665 5.92871, 118.27581 5.92871...
2	3	4	POLYGON ((118.27581 6.37528, 118.72497 6.37528...

grid_gdf.dtypes

x              int64
y              int64
geometry    geometry
dtype: object

Spatial join with highest intersection

---

get_highest_intersection

def get_highest_intersection(
    gdf1:GeoDataFrame, # gdf1 will be the basis of output geometry
    gdf2:GeoDataFrame, # gdf2 data will all be included during intersection
    proj_crs:str, # metric CRS (e.g., Philippines uses EPSG:32651)
)->GeoDataFrame:

Gets the intersection based on the largest area joined

get_highest_intersection(grid_gdf, admin_bounds_gdf, "EPSG:32651")

	geometry	x	y	shapeName	shapeISO	shapeID	shapeGroup	shapeType
0	POLYGON ((118.27581 5.92871, 118.72497 5.92871...	3	3	Tawi-Tawi	None	PHL-ADM2-3_0_0-B77	PHL	ADM2
1	POLYGON ((117.82665 5.92871, 118.27581 5.92871...	2	3	Tawi-Tawi	None	PHL-ADM2-3_0_0-B77	PHL	ADM2
2	POLYGON ((118.27581 6.37528, 118.72497 6.37528...	3	4	Tawi-Tawi	None	PHL-ADM2-3_0_0-B77	PHL	ADM2
3	POLYGON ((118.27581 6.82146, 118.72497 6.82146...	3	5	Tawi-Tawi	None	PHL-ADM2-3_0_0-B77	PHL	ADM2
4	POLYGON ((116.92834 7.26723, 117.37749 7.26723...	0	6	Palawan	None	PHL-ADM2-3_0_0-B59	PHL	ADM2
...	...	...	...	...	...	...	...	...
318	POLYGON ((125.46233 5.03451, 125.91149 5.03451...	19	1	Davao del Sur	None	PHL-ADM2-3_0_0-B28	PHL	ADM2
319	POLYGON ((125.46233 9.93163, 125.91149 9.93163...	19	12	Dinagat Islands	None	PHL-ADM2-3_0_0-B30	PHL	ADM2
320	POLYGON ((125.46233 10.37375, 125.91149 10.373...	19	13	Eastern Samar	None	PHL-ADM2-3_0_0-B31	PHL	ADM2
321	POLYGON ((125.91149 9.93163, 126.36065 9.93163...	20	12	Surigao del Norte	None	PHL-ADM2-3_0_0-B74	PHL	ADM2
322	POLYGON ((125.91149 10.37375, 126.36065 10.373...	20	13	Eastern Samar	None	PHL-ADM2-3_0_0-B31	PHL	ADM2

323 rows × 8 columns

Check plot of output

output = get_highest_intersection(grid_gdf, admin_bounds_gdf, "EPSG:32651")

output.plot()

# output.explore()