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

Generate a sample grid

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

%%time
grid_gdf = grid_generator5k.generate_grid(admin_bounds_gdf)

CPU times: user 3.67 s, sys: 19.1 ms, total: 3.69 s
Wall time: 3.69 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[source]

get_highest_intersection(gdf1:GeoDataFrame, gdf2:GeoDataFrame, proj_crs:str)

Gets the intersection based on the largest area joined

	Type	Default	Details
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)

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()

<AxesSubplot:>