GeoPandas Geospatial Data Analysis

Python library for geospatial vector data – extends pandas with spatial operations.

When to Use

• Working with geographic/spatial data (shapefiles, GeoJSON, GeoPackage)
• Spatial analysis (buffer, intersection, spatial joins)
• Coordinate transformations and projections
• Creating choropleth maps
• Processing geographic boundaries, points, lines, polygons

Quick Start

import geopandas as gpd

# Read spatial data
gdf = gpd.read_file("data.geojson")

# Basic exploration
print(gdf.head())
print(gdf.crs)  # Coordinate Reference System
print(gdf.geometry.geom_type)

# Simple plot
gdf.plot()

# Reproject to different CRS
gdf_projected = gdf.to_crs("EPSG:3857")

# Calculate area (use projected CRS)
gdf_projected['area'] = gdf_projected.geometry.area

# Save to file
gdf.to_file("output.gpkg")

Reading/Writing Data

# Read various formats
gdf = gpd.read_file("data.shp")       # Shapefile
gdf = gpd.read_file("data.geojson")   # GeoJSON
gdf = gpd.read_file("data.gpkg")      # GeoPackage

# Read with spatial filter (faster for large files)
gdf = gpd.read_file("data.gpkg", bbox=(xmin, ymin, xmax, ymax))

# Write to file
gdf.to_file("output.gpkg")
gdf.to_file("output.geojson", driver="GeoJSON")

# PostGIS database
from sqlalchemy import create_engine
engine = create_engine("postgresql://user:pass@localhost/db")
gdf = gpd.read_postgis("SELECT * FROM table", con=engine, geom_col='geom')

Coordinate Reference Systems

# Check CRS
print(gdf.crs)

# Set CRS (when metadata missing)
gdf = gdf.set_crs("EPSG:4326")

# Reproject (transforms coordinates)
gdf_projected = gdf.to_crs("EPSG:3857")  # Web Mercator
gdf_projected = gdf.to_crs("EPSG:32633")  # UTM zone 33N

# Common CRS codes:
# EPSG:4326 - WGS84 (lat/lon)
# EPSG:3857 - Web Mercator
# EPSG:326XX - UTM zones

Geometric Operations

# Buffer (expand/shrink geometries)
buffered = gdf.geometry.buffer(100)  # 100 units buffer

# Centroid
centroids = gdf.geometry.centroid

# Simplify (reduce vertices)
simplified = gdf.geometry.simplify(tolerance=5, preserve_topology=True)

# Convex hull
hull = gdf.geometry.convex_hull

# Boundary
boundary = gdf.geometry.boundary

# Area and length (use projected CRS!)
gdf['area'] = gdf.geometry.area
gdf['length'] = gdf.geometry.length

Spatial Analysis

Spatial Joins

# Join based on spatial relationship
joined = gpd.sjoin(gdf1, gdf2, predicate='intersects')
joined = gpd.sjoin(gdf1, gdf2, predicate='within')
joined = gpd.sjoin(gdf1, gdf2, predicate='contains')

# Nearest neighbor join
nearest = gpd.sjoin_nearest(gdf1, gdf2, max_distance=1000)

Overlay Operations

# Intersection
intersection = gpd.overlay(gdf1, gdf2, how='intersection')

# Union
union = gpd.overlay(gdf1, gdf2, how='union')

# Difference
difference = gpd.overlay(gdf1, gdf2, how='difference')

Dissolve (Aggregate by Attribute)

# Merge geometries by attribute
dissolved = gdf.dissolve(by='region', aggfunc='sum')

Clip

# Clip data to boundary
clipped = gpd.clip(gdf, boundary_gdf)

Visualization

import matplotlib.pyplot as plt

# Basic plot
gdf.plot()

# Choropleth map
gdf.plot(column='population', cmap='YlOrRd', legend=True)

# Multi-layer map
fig, ax = plt.subplots(figsize=(10, 10))
gdf1.plot(ax=ax, color='blue', alpha=0.5)
gdf2.plot(ax=ax, color='red', alpha=0.5)
plt.savefig('map.png', dpi=300, bbox_inches='tight')

# Interactive map (requires folium)
gdf.explore(column='population', legend=True)

Common Workflows

Spatial Join and Aggregate

# Join points to polygons
points_in_polygons = gpd.sjoin(points_gdf, polygons_gdf, predicate='within')

# Aggregate by polygon
aggregated = points_in_polygons.groupby('index_right').agg({
    'value': 'sum',
    'count': 'size'
})

# Merge back to polygons
result = polygons_gdf.merge(aggregated, left_index=True, right_index=True)

Buffer Analysis

# Create buffers around points
gdf_projected = points_gdf.to_crs("EPSG:3857")  # Project first!
gdf_projected['buffer'] = gdf_projected.geometry.buffer(1000)  # 1km buffer
gdf_projected = gdf_projected.set_geometry('buffer')

# Find features within buffer
within_buffer = gpd.sjoin(other_gdf, gdf_projected, predicate='within')

Best Practices

1. Always check CRS before spatial operations
2. Use projected CRS for area/distance calculations
3. Match CRS before spatial joins or overlays
4. Validate geometries with .is_valid before operations
5. Use GeoPackage format over Shapefile (modern, better)
6. Use .copy() when modifying geometry to avoid side effects
7. Filter during read with bbox for large files

vs Alternatives

Resources

• Docs: https://geopandas.org/
• User Guide: https://geopandas.org/en/stable/docs/user_guide.html
• Gallery: https://geopandas.org/en/stable/gallery/index.html