Post #3: Spatial Joins and Attribute Queries using GeoPandas

 

Title: How to Perform Spatial Joins and Attribute Filters with Python and GeoPandas

---

📍Introduction

Spatial joins are at the heart of GIS analysis. They allow you to combine layers based on their location—like finding which villages fall inside a floodplain, or assigning administrative regions to road segments.

With GeoPandas, spatial joins and attribute filtering are easy, fast, and scriptable. In this post, we’ll walk through a simple example using two shapefiles:
📌 Polygons (e.g., districts)
📌 Points (e.g., schools, hospitals, wells, etc.)

---

🧰 Step 1: Load the Data

python
import geopandas as gpd

# Load polygons and points
districts = gpd.read_file("data/districts.shp")
wells = gpd.read_file("data/water_wells.shp")

Make sure both layers have the same CRS:

python
districts = districts.to_crs("EPSG:4326")
wells = wells.to_crs("EPSG:4326")

---

🔗 Step 2: Perform a Spatial Join

Join wells to districts to find out which well lies in which district:

python
joined = gpd.sjoin(wells, districts, how="left", predicate="within")

Now each well will have the attributes of the district it falls within!

python
print(joined.head())

---

🧠 What Just Happened?

• how="left": Keep all wells, even if they don’t fall in a district.

• predicate="within": Spatial relationship used to join features.

• Other predicates: "intersects", "contains", "touches"

---

🔍 Step 3: Run Attribute Queries

Filter wells that fall in a specific district:

python
riyadh_wells = joined[joined["district_name"] == "Riyadh"]
print(f"Number of wells in Riyadh: {len(riyadh_wells)}")

You can also filter based on depth or any other attribute:

python
deep_wells = riyadh_wells[riyadh_wells["depth"] > 100]

---

💾 Step 4: Save the Result

Export the result to a new shapefile or GeoJSON:

python
riyadh_wells.to_file("output/riyadh_wells.shp")

Or as GeoJSON:

python
riyadh_wells.to_file("output/riyadh_wells.geojson", driver="GeoJSON")

---

📊 Bonus: Count Wells per District

If you want a summary table (non-spatial):

python
counts = joined.groupby("district_name").size().reset_index(name="well_count")
print(counts)

---

🧠 Use Cases for Spatial Joins:

• Link building points to land parcels

• Assign municipalities to GPS points

• Match roads to regions

• Join elevation points to geology zones

• Count features within a boundary

---

🎯 Conclusion

Spatial joins let you combine geospatial layers based on real-world relationships. With GeoPandas, you can script these workflows, run them in batches, and make your GIS work smarter—not harder.

---

📌 Next Up:

➡️ Post 4: Plotting and Styling GIS Data with Python