Geopandas

GeoPandas Skill

geopandas spatial data library for Python: manipulation, analysis, and visualization of geographic data. Covers GeoDataFrames, spatial joins, CRS/projections, vector operations, raster integration (rasterio, xarray), choropleth mapping, interactive maps (folium), basemap tiles (contextily), spatial autocorrelation, and the PySAL ecosystem. Use when working with geographic data, reading/writing spatial files (Shapefile, GeoPackage, GeoParquet), making maps, or running spatial statistics. For interactive web-based geographic charts without spatial analysis, use plotly.

Comprehensive skill for spatial data analysis with geopandas and the broader Python geospatial stack. Use the decision trees below to find the right guidance, then load detailed references as needed.

Version Notes

This skill targets geopandas 1.x (tested with 1.1.3). Key changes from earlier versions:

• Shapely >= 2.0 required (PyGEOS backend removed, vectorized ops built-in)
• pyogrio is the default I/O engine (replacing fiona, 5-10x faster)
• cascaded_union removed — use union_all() instead
• GeoSeries.unary_union property renamed to GeoSeries.union_all() method

What is GeoPandas?

GeoPandas extends pandas with spatial data types and operations:

• GeoDataFrame: A pandas DataFrame with a geometry column — tabular data meets spatial operations
• Spatial operations: Joins, overlays, dissolve, clip, buffer, and distance calculations on vector geometries
• CRS handling: Coordinate reference system management via pyproj for correct spatial computations
• Visualization: Static maps (matplotlib), interactive maps (folium via .explore()), and GPU-accelerated rendering (lonboard)
• Ecosystem hub: Integrates with PySAL (spatial statistics), rasterio (rasters), contextily (basemaps), and mapclassify (classification schemes)

How to Use This Skill

Reference File Structure

Reading Order

1. New to geopandas? Start with quickstart.md then spatial-operations.md 2. Making maps? Read visualization.md (relies on crs-projections.md for projection choices) 3. Spatial statistics? Read pysal-spatial-stats.md (for methodology context, also load data-scientist skill's geospatial-analysis.md) 4. Having issues? Check gotchas.md first

Related Skills

• data-scientist (geospatial-analysis.md, geospatial-operations.md): Spatial methodology — when/why to use spatial methods, interpretation guidance, MAUP, ecological fallacy. Load alongside this skill for research workflows.
• polars: If spatial data is combined with large tabular datasets, use polars for non-spatial transformations before converting to GeoDataFrame.
• plotnine / plotly: For non-map visualizations of spatial analysis results (coefficient plots, distributions).

Quick Decision Trees

"I need to read or write spatial data"

Loading/saving spatial data?
├─ Read vector file (Shapefile, GeoPackage, GeoJSON) → ./references/data-io.md
├─ Read GeoParquet → ./references/data-io.md
├─ Read from PostGIS / DuckDB Spatial → ./references/data-io.md
├─ Download boundaries (Census, OSM) → ./references/data-io.md
├─ Create GeoDataFrame from lat/lon columns → ./references/quickstart.md
├─ Write to file → ./references/data-io.md
└─ Read raster data (GeoTIFF) → ./references/raster-integration.md

"I need to combine or transform spatial data"

Spatial operations?
├─ Join by location (point-in-polygon, etc.) → ./references/spatial-operations.md
├─ Join by nearest feature → ./references/spatial-operations.md
├─ Overlay (intersection, union, difference) → ./references/spatial-operations.md
├─ Dissolve (merge polygons by attribute) → ./references/spatial-operations.md
├─ Clip to boundary → ./references/spatial-operations.md
├─ Buffer features → ./references/spatial-operations.md
├─ Compute distances → ./references/spatial-operations.md
├─ Compute centroids or areas → ./references/spatial-operations.md
└─ Areal interpolation (mismatched boundaries) → ./references/spatial-operations.md

"I need to fix CRS or projection issues"

CRS/projection issues?
├─ Check current CRS → ./references/crs-projections.md
├─ Reproject to different CRS → ./references/crs-projections.md
├─ Choose a projection for analysis → ./references/crs-projections.md
├─ Data has no CRS (set it) → ./references/crs-projections.md
├─ CRS mismatch error → ./references/gotchas.md
└─ Area/distance calculations wrong → ./references/crs-projections.md

"I need to make a map"

Making maps?
├─ Quick static choropleth → ./references/visualization.md
├─ Classification schemes (quantiles, Fisher-Jenks) → ./references/visualization.md
├─ Add basemap tiles → ./references/visualization.md
├─ Interactive map (pan/zoom/hover) → ./references/visualization.md
├─ Large dataset (millions of features) → ./references/visualization.md
├─ Multi-panel / faceted maps → ./references/visualization.md
├─ LISA cluster map → ./references/pysal-spatial-stats.md
└─ Export to PNG/SVG/HTML → ./references/visualization.md

"I need spatial statistics"

Spatial statistics?
├─ Build spatial weights matrix → ./references/pysal-spatial-stats.md
├─ Test for spatial autocorrelation (Moran's I) → ./references/pysal-spatial-stats.md
├─ Find hot spots / cold spots (LISA) → ./references/pysal-spatial-stats.md
├─ Spatial regression (lag, error, Durbin) → ./references/pysal-spatial-stats.md
├─ Point pattern analysis → ./references/pysal-spatial-stats.md
└─ Methodology guidance (interpretation, MAUP) → data-scientist skill: geospatial-analysis.md

"I need to work with rasters"

Raster operations?
├─ Read GeoTIFF → ./references/raster-integration.md
├─ Zonal statistics (summarize raster by polygons) → ./references/raster-integration.md
├─ Clip/mask raster by polygon → ./references/raster-integration.md
├─ Extract raster values at points → ./references/raster-integration.md
├─ Multidimensional raster (xarray) → ./references/raster-integration.md
└─ Rasterize vectors / vectorize rasters → ./references/raster-integration.md

"Something isn't working"

Having issues?
├─ CRS mismatch errors → ./references/gotchas.md
├─ Invalid geometry errors → ./references/gotchas.md
├─ Spatial join produced wrong row count → ./references/gotchas.md
├─ Memory issues with large files → ./references/gotchas.md
├─ Shapely version confusion → ./references/gotchas.md
├─ Coordinate order (lon/lat vs lat/lon) → ./references/gotchas.md
└─ General troubleshooting → ./references/gotchas.md

File-First Execution in Research Workflows

Important: In data research pipelines (see CLAUDE.md), spatial operations are executed through script files, not interactively. This ensures auditability and reproducibility.

The pattern: 1. Write spatial analysis code to scripts/stage{N}_{type}/{step}_{task-name}.py 2. Execute via Bash with automatic output capture wrapper script 3. Validation results get automatically embedded in scripts as comments 4. If failed, create versioned copy for fixes

Closely read agent_reference/SCRIPT_EXECUTION_REFERENCE.md for the mandatory file-first execution protocol covering complete code file writing, output capture, and file versioning rules.

See:

• agent_reference/SCRIPT_EXECUTION_REFERENCE.md — Script execution protocol and format with validation

The examples in reference files show geopandas syntax. In research workflows, wrap them in scripts following the file-first pattern.

---

Quick Reference

Essential Imports

import geopandas as gpd
from shapely.geometry import Point, Polygon, LineString, MultiPolygon

Core Operations

Common CRS Codes

Topic Index

Citation

When this library is used as a primary analytical tool, include in the report's Software & Tools references:

Jordahl, K. et al. geopandas: Python tools for geographic data [Computer software]. https://geopandas.org/

Cite when: geopandas is used for spatial operations, spatial joins, or map visualization central to the analysis. Do not cite when: Only used to read a shapefile for a simple reference lookup.

If PySAL spatial analysis functions are also used (spatial weights, Moran's I, etc.), additionally cite:

Rey, S.J. et al. (2022). "The PySAL Ecosystem of Open-Source Python Packages for the Analysis of Spatial Data." Geographical Analysis, 54(3), 467-487.

For method-specific citations (e.g., spatial statistics techniques), consult the reference files in this skill and agent_reference/CITATION_REFERENCE.md.

CRS and Projections

Coordinate Reference System handling in geopandas — checking, setting, transforming, and choosing projections. Getting the CRS right is a prerequisite for correct spatial operations; getting it wrong produces silently wrong results.

---

CRS Fundamentals

Geographic vs Projected CRS

The critical rule: Buffer, area, distance, and centroid operations on a geographic CRS (longitude/latitude) produce wrong or misleading results because one degree of longitude varies from ~111 km at the equator to ~0 km at the poles. Always reproject to an appropriate projected CRS before computing.

---

Checking CRS

# Check the current CRS
print(gdf.crs)
# Output: EPSG:4326 (or the full WKT2 string)

# Check EPSG code
print(gdf.crs.to_epsg())
# Output: 4326

# Check if geographic (lon/lat) or projected
print(gdf.crs.is_geographic)
# Output: True

print(gdf.crs.is_projected)
# Output: False

# Check the axis units
print(gdf.crs.axis_info)
# Shows units (degree vs metre)

# Detailed CRS information
print(gdf.crs.to_wkt(pretty=True))

---

Setting CRS (Declaring, Not Transforming)

Setting a CRS tells geopandas what system the coordinates are already in. No coordinate values change. Use this when data arrives without CRS metadata (common with CSV files containing lat/lon columns).

# Set CRS on creation
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lon, df.lat), crs="EPSG:4326")

# Set CRS on existing GeoDataFrame (no CRS currently set)
gdf = gdf.set_crs("EPSG:4326")

# Override an incorrect CRS (use cautiously — you're saying "the existing CRS label is wrong")
gdf = gdf.set_crs("EPSG:4326", allow_override=True)

When to use `set_crs`: Only when the GeoDataFrame has no CRS (gdf.crs is None) or when you know the existing CRS label is incorrect. If the data already has a CRS and you want to change the projection, use to_crs instead.

---

Reprojecting (Transforming Coordinates)

Reprojecting recomputes all coordinate values from one CRS to another. Use this when you need a different projection for analysis or visualization.

# Reproject to NAD83 Conus Albers (equal-area, good for US maps)
gdf_albers = gdf.to_crs(epsg=5070)

# Reproject using CRS string
gdf_utm = gdf.to_crs("EPSG:32617")

# Reproject using PROJ string (less common)
gdf_custom = gdf.to_crs("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96")

# Reproject to match another GeoDataFrame's CRS
gdf_matched = gdf.to_crs(other_gdf.crs)

Setting vs Transforming: The Difference

# WRONG: Using set_crs to "reproject" — coordinates don't change, CRS label changes
# This makes the data plot in the wrong location
gdf_wrong = gdf.set_crs("EPSG:5070")  # DON'T DO THIS if gdf is in EPSG:4326

# RIGHT: Using to_crs to reproject — coordinates are recomputed
gdf_right = gdf.to_crs(epsg=5070)     # Coordinates change, CRS changes

---

Choosing a Projection

Decision Guide

What does your analysis need?
├─ Area calculations or thematic maps
│   └─ Equal-area projection
│       ├─ Continental US → EPSG:5070 (NAD83 Conus Albers)
│       ├─ Single state → State Plane (equal-area variant) or LAEA
│       ├─ Global → Mollweide or Equal Earth
│       └─ Custom study area → LAEA centered on study centroid
├─ Distance or local-scale analysis (<500 km)
│   └─ UTM zone for study area
│       Use pyproj to find the right zone:
│       from pyproj import CRS
│       CRS.from_authority("EPSG", pyproj.database.query_utm_crs_info(
│           datum_name="WGS 84", area_of_interest=pyproj.aoi.AreaOfInterest(
│               west_lon_degree=-77, south_lat_degree=38,
│               east_lon_degree=-76, north_lat_degree=39))[0].code)
├─ Web map display
│   └─ Web Mercator (EPSG:3857) — display only, never for analysis
├─ Data storage or exchange
│   └─ WGS84 (EPSG:4326)
└─ Unsure
    └─ Start with EPSG:5070 for US, UTM for local

Common US Projections

State Plane Coordinate Systems

Each US state has one or more State Plane zones optimized for local accuracy. Find the right one:

import pyproj

# List available State Plane CRS for a state
results = pyproj.database.query_crs_info(
    auth_name="EPSG",
    area_of_interest=pyproj.aoi.AreaOfInterest(
        west_lon_degree=-78.0, south_lat_degree=38.0,
        east_lon_degree=-75.0, north_lat_degree=40.0
    ),
    crs_types=["PROJECTED_CRS"]
)
# Filter results for "State Plane" in the name

Custom LAEA (Lambert Azimuthal Equal-Area)

For study areas not well-served by standard projections, create a custom equal-area projection centered on your data:

# Center the projection on your study area
centroid = gdf.to_crs(epsg=4326).dissolve().centroid.iloc[0]
custom_crs = f"+proj=laea +lat_0={centroid.y} +lon_0={centroid.x} +datum=WGS84 +units=m"
gdf_custom = gdf.to_crs(custom_crs)

---

CRS Matching Before Spatial Operations

All spatial operations (join, overlay, distance, etc.) require both inputs to be in the same CRS. geopandas raises a CRSMismatchError if they differ.

# Check if two GeoDataFrames share the same CRS
if gdf1.crs == gdf2.crs:
    result = gpd.sjoin(gdf1, gdf2)
else:
    # Reproject one to match the other
    gdf2_reprojected = gdf2.to_crs(gdf1.crs)
    result = gpd.sjoin(gdf1, gdf2_reprojected)

Best Practice: Reproject Early

# Standard workflow: reproject all inputs to a common CRS at the start
TARGET_CRS = "EPSG:5070"

counties = gpd.read_file("counties.gpkg").to_crs(TARGET_CRS)
schools = gpd.read_file("schools.gpkg").to_crs(TARGET_CRS)
tracts = gpd.read_file("tracts.gpkg").to_crs(TARGET_CRS)

# Now all spatial operations work without CRS concerns

---

pyproj CRS Objects

geopandas uses pyproj for CRS handling. The pyproj.CRS object provides rich CRS information:

from pyproj import CRS

# Create from EPSG
crs = CRS.from_epsg(5070)

# Create from PROJ string
crs = CRS.from_proj4("+proj=aea +lat_1=29.5 +lat_2=45.5 +lon_0=-96 +datum=NAD83")

# Create from WKT
crs = CRS.from_wkt(wkt_string)

# Inspect properties
print(crs.name)             # NAD83 / Conus Albers
print(crs.is_geographic)    # False
print(crs.is_projected)     # True
print(crs.axis_info)        # Axis names, directions, units
print(crs.area_of_use)      # Valid geographic extent
print(crs.to_epsg())        # 5070

# Get UTM zone for a point
from pyproj import database
utm_results = database.query_utm_crs_info(
    datum_name="WGS 84",
    area_of_interest=pyproj.aoi.AreaOfInterest(
        west_lon_degree=-77.1, south_lat_degree=38.8,
        east_lon_degree=-76.9, north_lat_degree=39.0
    )
)
print(f"EPSG:{utm_results[0].code}")  # e.g., EPSG:32618

---

Common CRS Pitfalls

---

References and Further Reading

Dorman, M., Graser, A., Nowosad, J., and Lovelace, R. (2025). Geocomputation with Python, Chs. 1 and 6: "Geographic data" and "Reprojecting geographic data." https://py.geocompx.org/

Tenkanen, H., Heikinheimo, V., and Whipp, D. (2024). Introduction to Python for Geographic Data Analysis, Ch. 5: "Map projections." https://pythongis.org/

pyproj documentation. https://pyproj4.github.io/pyproj/

Battersby, S. (2017). "Map Projections." The Geographic Information Science & Technology Body of Knowledge. https://gistbok.ucgis.org/

Spatial Data I/O

Reading and writing spatial data in geopandas — file formats, engines, web data sources, and spatial databases.

---

File Format Comparison

---

The pyogrio Engine

Since geopandas 1.0, pyogrio is the default I/O engine (replacing fiona). It is 5-10x faster for most operations because it uses GDAL's vectorized column-oriented reading.

# pyogrio is used automatically — no explicit engine needed
gdf = gpd.read_file("data.gpkg")

# Explicitly specify engine if needed
gdf = gpd.read_file("data.gpkg", engine="pyogrio")

# Use Arrow for even faster reads (returns Arrow-backed geometries)
gdf = gpd.read_file("data.gpkg", engine="pyogrio", use_arrow=True)

# Fall back to fiona if needed (must be installed separately)
gdf = gpd.read_file("data.gpkg", engine="fiona")

pyogrio Direct Usage

For maximum control or when not using GeoDataFrames:

import pyogrio

# List layers in a file
pyogrio.list_layers("multi_layer.gpkg")

# Read file info without loading data
info = pyogrio.read_info("data.gpkg")
print(info)  # CRS, geometry type, feature count, bounds, etc.

# Read to Arrow for interop with other tools
table = pyogrio.read_arrow("data.gpkg")

---

Reading Vector Data

Standard Read

gdf = gpd.read_file("counties.gpkg")

Filtered Reads (Essential for Large Files)

# Bounding box filter — only features intersecting the box
gdf = gpd.read_file("us_counties.gpkg", bbox=(-80, 35, -75, 40))

# Geometry mask — only features intersecting a geometry
from shapely.geometry import box
mask = box(-80, 35, -75, 40)
gdf = gpd.read_file("us_counties.gpkg", mask=mask)

# Column filter — read only needed columns (reduces memory)
gdf = gpd.read_file("us_counties.gpkg", columns=["GEOID", "NAME", "POP", "geometry"])

# Row limit — read a sample
gdf = gpd.read_file("us_counties.gpkg", rows=100)

# Specific layer from multi-layer file
gdf = gpd.read_file("census.gpkg", layer="tracts")

# SQL-based filter (with pyogrio engine)
gdf = gpd.read_file("us_counties.gpkg", where="STATE_FIPS = '06'")

Reading GeoParquet

GeoParquet is the best format for analytical workflows — columnar storage enables fast column selection, predicate pushdown, and compact storage.

gdf = gpd.read_parquet("counties.parquet")

# Column selection (only reads selected columns from disk)
gdf = gpd.read_parquet("counties.parquet", columns=["GEOID", "NAME", "geometry"])

# Bounding box filter
gdf = gpd.read_parquet("counties.parquet", bbox=(-80, 35, -75, 40))

# From cloud storage (requires fsspec + storage backend)
gdf = gpd.read_parquet("s3://bucket/counties.parquet")

Reading GeoFeather

gdf = gpd.read_feather("counties.feather")
gdf = gpd.read_feather("counties.feather", columns=["GEOID", "geometry"])

---

Writing Vector Data

To Standard Formats

# GeoPackage (recommended)
gdf.to_file("output.gpkg", driver="GPKG")

# GeoPackage with specific layer name
gdf.to_file("output.gpkg", driver="GPKG", layer="counties")

# Append to existing GeoPackage
gdf.to_file("output.gpkg", driver="GPKG", layer="new_layer", mode="a")

# Shapefile (legacy)
gdf.to_file("output.shp")

# GeoJSON
gdf.to_file("output.geojson", driver="GeoJSON")

# FlatGeobuf
gdf.to_file("output.fgb", driver="FlatGeobuf")

To GeoParquet (Preferred for Pipelines)

gdf.to_parquet("output.parquet")

# With compression
gdf.to_parquet("output.parquet", compression="snappy")  # default
gdf.to_parquet("output.parquet", compression="zstd")     # better compression

To GeoFeather

gdf.to_feather("output.feather")

To Other Formats

# GeoJSON string (for web APIs)
geojson_str = gdf.to_json()

# Well-Known Text (WKT)
gdf["wkt"] = gdf.geometry.to_wkt()

# Well-Known Binary (WKB)
gdf["wkb"] = gdf.geometry.to_wkb()

---

Web Data Sources

OpenStreetMap via osmnx

import osmnx as ox

# Download street network
G = ox.graph_from_place("Washington, DC", network_type="drive")

# Download building footprints
buildings = ox.features_from_place("Washington, DC", tags={"building": True})

# Download specific POIs
schools = ox.features_from_place("Washington, DC", tags={"amenity": "school"})

# Download administrative boundary
dc = ox.geocode_to_gdf("Washington, DC")

Census Boundaries via pygris

import pygris

# Download county boundaries
counties = pygris.counties(year=2022)

# Download census tracts for a state
tracts = pygris.tracts(state="DC", year=2022)

# Download school districts
districts = pygris.school_districts(state="VA", year=2022)

# Download block groups
bgs = pygris.block_groups(state="MD", county="031", year=2022)

Census Boundaries via Direct URL

# TIGER/Line Shapefiles (direct download)
url = "https://www2.census.gov/geo/tiger/TIGER2022/COUNTY/tl_2022_us_county.zip"
counties = gpd.read_file(url)

Natural Earth (Global Boundaries)

# Via URL
url = "https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip"
world = gpd.read_file(url)

WFS (Web Feature Service)

# Example WFS endpoint
wfs_url = "https://example.com/wfs?service=WFS&version=2.0.0&request=GetFeature&typeName=layer_name&outputFormat=application/json"
gdf = gpd.read_file(wfs_url)

---

Spatial Databases

PostGIS

from sqlalchemy import create_engine

engine = create_engine("postgresql://user:pass@host:5432/dbname")

# Read from PostGIS
gdf = gpd.read_postgis(
    "SELECT * FROM counties WHERE state_fips = '06'",
    con=engine,
    geom_col="geom"
)

# Write to PostGIS
gdf.to_postgis("counties_clean", engine, if_exists="replace")

DuckDB Spatial

import duckdb

con = duckdb.connect()
con.install_extension("spatial")
con.load_extension("spatial")

# Read spatial data via DuckDB
result = con.execute("""
    SELECT *, ST_AsWKB(geom) as geometry
    FROM read_parquet('counties.parquet')
    WHERE state_fips = '06'
""").fetchdf()

# Convert to GeoDataFrame
gdf = gpd.GeoDataFrame(
    result.drop(columns=["geometry"]),
    geometry=gpd.GeoSeries.from_wkb(result["geometry"]),
    crs="EPSG:4326"
)

---

Format Selection Guide

What's the use case?
├─ Analytical pipeline (parquet ecosystem) → GeoParquet
├─ GIS interoperability (QGIS, ArcGIS) → GeoPackage
├─ Web display / API response → GeoJSON
├─ Streaming large files over network → FlatGeobuf
├─ Local temporary cache → GeoFeather
├─ Must support legacy GIS tools → Shapefile (reluctantly)
└─ Multi-layer container → GeoPackage

---

References and Further Reading

Jordahl, K. et al. (2024). geopandas I/O documentation. https://geopandas.org/en/stable/docs/user_guide/io.html

pyogrio documentation. https://pyogrio.readthedocs.io/

GeoParquet specification. https://geoparquet.org/

Tenkanen, H., Heikinheimo, V., and Whipp, D. (2024). Introduction to Python for Geographic Data Analysis, Ch. 6: "Reading and writing spatial data." https://pythongis.org/

Boeing, G. (2017). "OSMnx: New methods for acquiring, constructing, analyzing, and visualizing complex street networks." Computers, Environment and Urban Systems, 65, 126-139.

Common Gotchas and Troubleshooting

Frequent issues when working with geopandas and the spatial Python stack — symptoms, causes, and fixes.

---

CRS Mismatch Errors

Symptom

ValueError: GeoDataFrame does not have a CRS set.
# or
CRSMismatchError: CRS mismatch between the CRS of left geometries and the CRS of right geometries.

Cause

Spatial operations (join, overlay, distance) require all inputs in the same CRS. This error fires when CRS is missing or doesn't match.

Fix

# Check CRS of both inputs
print(gdf1.crs)  # e.g., EPSG:4326
print(gdf2.crs)  # e.g., EPSG:5070 or None

# If missing CRS, set it (don't reproject — set declares what system the data is already in)
gdf2 = gdf2.set_crs("EPSG:4326")

# If mismatched, reproject one to match the other
gdf2 = gdf2.to_crs(gdf1.crs)

# Now the spatial operation will work
result = gpd.sjoin(gdf1, gdf2)

Prevention

Reproject all inputs to a common CRS at the start of the script:

TARGET_CRS = "EPSG:5070"
gdf1 = gdf1.to_crs(TARGET_CRS)
gdf2 = gdf2.to_crs(TARGET_CRS)

---

Invalid Geometries

Symptom

• Overlay or spatial join returns empty or incomplete results
• TopologyException: found non-noded intersection error
• Unexpected geometry fragments after overlay

Cause

Invalid geometries violate the Simple Features specification — self-intersections, unclosed rings, duplicate vertices, etc.

Detection

# Check which geometries are invalid
invalid_mask = ~gdf.geometry.is_valid
print(f"Invalid geometries: {invalid_mask.sum()} / {len(gdf)}")

# Inspect a sample
if invalid_mask.any():
    from shapely.validation import explain_validity
    for idx in gdf[invalid_mask].index[:5]:
        print(f"  Row {idx}: {explain_validity(gdf.loc[idx, 'geometry'])}")

Fix

# Option 1: make_valid() — preferred (Shapely 2.0+)
gdf["geometry"] = gdf.geometry.make_valid()

# Option 2: Buffer by zero — classic approach
gdf["geometry"] = gdf.geometry.buffer(0)

# Verify fix
assert gdf.geometry.is_valid.all(), "Some geometries still invalid after repair"

When to Validate

• Before any overlay operation
• Before spatial joins with contains or within predicates
• After reading Shapefiles (the format is prone to validity issues)
• After any geometry transformation (simplification, projection, dissolve)

---

Spatial Join Produces Wrong Row Count

Symptom

After gpd.sjoin(), the result has more rows than the left input (unexpected duplicates) or fewer rows (unexpected drops).

Cause: Many-to-Many Matches

A point on a polygon boundary may match multiple polygons (with intersects). Overlapping polygons cause a single feature to match several targets.

Diagnosis

print(f"Input rows: {len(left_gdf)}")
print(f"Result rows: {len(result)}")
print(f"Duplicated indices: {result.index.duplicated().sum()}")
print(f"Null join columns: {result['right_col'].isna().sum()}")

Fix

# Strategy 1: Use stricter predicate
result = gpd.sjoin(points, polygons, predicate="within")  # Instead of "intersects"

# Strategy 2: Drop duplicates (keep first match)
result = result[~result.index.duplicated(keep="first")]

# Strategy 3: Aggregate after join
result = result.groupby(result.index).agg({
    "county_name": "first",
    "population": "sum"
})

Cause: Unmatched Features (Fewer Rows)

Default how="inner" drops features with no spatial match.

# Use left join to keep all left features
result = gpd.sjoin(schools, counties, predicate="within", how="left")
# Unmatched schools will have NaN in county columns

---

Memory Issues with Large Shapefiles

Symptom

Out of memory when reading a large Shapefile, or very slow performance.

Fixes

# 1. Read only needed columns
gdf = gpd.read_file("huge.shp", columns=["GEOID", "NAME", "geometry"])

# 2. Read only a geographic subset
gdf = gpd.read_file("huge.shp", bbox=(-80, 35, -75, 40))

# 3. Use pyogrio with Arrow backend (lower memory)
gdf = gpd.read_file("huge.shp", engine="pyogrio", use_arrow=True)

# 4. Use GeoParquet instead (much faster reads, column selection)
# Convert once:
gdf_full = gpd.read_file("huge.shp")
gdf_full.to_parquet("huge.parquet")
# Then read efficiently:
gdf = gpd.read_parquet("huge.parquet", columns=["GEOID", "geometry"])

# 5. For visualization of millions of features, use lonboard instead of matplotlib
from lonboard import viz
m = viz(gdf)

---

Shapely 1.x vs 2.x Differences

geopandas 1.x requires Shapely 2.0+, which introduced significant changes.

Key Changes

What Broke

# Shapely 1.x (no longer works)
from shapely.ops import cascaded_union    # Removed
cascaded_union(geom_list)

# Shapely 2.x equivalent
from shapely import union_all
union_all(geom_list)
# Or in geopandas:
gdf.geometry.union_all()

# Shapely 1.x
geom.is_empty  # Was an attribute in 1.x, still works in 2.x

# The PyGEOS backend option is gone — Shapely 2.x IS the vectorized engine
# Delete any old environment variable:
# unset USE_PYGEOS  (no longer needed)

If You See Deprecation Warnings

Some Shapely 1.x functions were removed; others still work but have modern equivalents:

# Removed (will raise ImportError)
from shapely.ops import cascaded_union  # Use shapely.union_all() instead

# Deprecated — use the modern equivalents
from shapely.ops import unary_union     # Deprecated in shapely; use union_all()
from shapely import union_all           # Preferred vectorized equivalent
# In geopandas: gdf.geometry.union_all() (unary_union is deprecated in 1.1.3)

# Repair pattern
geom.buffer(0)                # Classic fix — still works
shapely.make_valid(geom)      # Preferred in 2.x (more predictable)

---

Coordinate Order Confusion (lon/lat vs lat/lon)

The Problem

• Mathematics/GIS convention: x = longitude, y = latitude → (lon, lat) = (-77.036, 38.901)
• Everyday convention: "latitude and longitude" → people say (38.901, -77.036)
• Google Maps URLs: @38.901,-77.036 (lat, lon)
• GeoJSON spec: [longitude, latitude] (lon, lat)

Symptoms of Getting It Wrong

• Points plot in the ocean or in the wrong hemisphere
• Spatial joins return zero matches
• Data appears reflected across the equator or prime meridian

How geopandas Expects Coordinates

# geopandas follows the GIS convention: x = longitude, y = latitude
gpd.points_from_xy(
    x=df["longitude"],    # x-axis = longitude (horizontal)
    y=df["latitude"]      # y-axis = latitude (vertical)
)

# Shapely Point: Point(x, y) = Point(longitude, latitude)
from shapely.geometry import Point
capitol = Point(-77.009, 38.890)  # (lon, lat)

Quick Diagnostic

# If total_bounds look wrong, coordinates may be swapped
print(gdf.total_bounds)
# Expected for US data: [-125, 24, -66, 50] (minx, miny, maxx, maxy)
# If you see [24, -125, 50, -66], lat/lon are swapped

# Fix: swap x and y
from shapely.ops import transform
gdf["geometry"] = gdf.geometry.map(lambda geom: transform(lambda x, y: (y, x), geom))

---

dtype Warnings with Mixed Geometry Types

Symptom

UserWarning: Geometry column does not contain geometry.

Or unexpected behavior when a GeoDataFrame contains mixed geometry types (e.g., both Point and MultiPoint, or Polygon and MultiPolygon).

Fix

# Check geometry types
print(gdf.geom_type.value_counts())

# Explode MultiGeometries to single parts
gdf = gdf.explode(index_parts=False)

# Or force to Multi type for consistency
from shapely.geometry import MultiPolygon

def to_multi(geom):
    if geom.geom_type == "Polygon":
        return MultiPolygon([geom])
    return geom

gdf["geometry"] = gdf.geometry.map(to_multi)

---

GeoDataFrame Loses CRS After Pandas Operations

Symptom

After a pandas operation (concat, merge, groupby), gdf.crs is None.

Cause

Some pandas operations return a plain DataFrame, losing the GeoDataFrame type and CRS.

Fix

# After pd.concat — re-wrap as GeoDataFrame
import pandas as pd
result = pd.concat([gdf1, gdf2])
result = gpd.GeoDataFrame(result, crs=gdf1.crs)

# After merge where geometry might be lost
result = gdf.merge(df, on="key")
if not isinstance(result, gpd.GeoDataFrame):
    result = gpd.GeoDataFrame(result, geometry="geometry", crs=gdf.crs)

# After groupby — dissolve preserves geometry; manual groupby does not
# Use dissolve instead of groupby for spatial data:
result = gdf.dissolve(by="state", aggfunc="sum")  # Preserves geometry + CRS

---

Plotting Issues

Map Appears Stretched or Distorted

# Add aspect ratio correction
ax.set_aspect("equal")

# Or use a projection that minimizes distortion for your study area
gdf.to_crs(epsg=5070).plot()

Legend Overlaps Map

gdf.plot(
    column="value",
    legend=True,
    legend_kwds={
        "loc": "lower right",
        "fontsize": 8,
        "shrink": 0.6,       # For continuous legends
        "pad": 0.02
    }
)

Basemap Not Showing

# Ensure data is in Web Mercator for basemap alignment
import contextily as ctx
ax = gdf.to_crs(epsg=3857).plot(alpha=0.5)
ctx.add_basemap(ax)

# Or pass CRS explicitly
ax = gdf.plot(alpha=0.5)
ctx.add_basemap(ax, crs=gdf.crs)

---

geopolars: Not Production-Ready

geopolars (a potential geopandas-like interface backed by Polars) is experimental and not suitable for production use. Stick with geopandas for spatial operations and convert to/from Polars for non-spatial data manipulation:

import polars as pl
import geopandas as gpd

# Polars → GeoDataFrame
df_polars = pl.read_parquet("data.parquet")
df_pandas = df_polars.to_pandas()
gdf = gpd.GeoDataFrame(df_pandas, geometry=gpd.points_from_xy(df_pandas.lon, df_pandas.lat), crs="EPSG:4326")

# GeoDataFrame → Polars (drop geometry, or convert to WKT/WKB)
df_polars = pl.from_pandas(gdf.drop(columns=["geometry"]))

---

References and Further Reading

geopandas FAQ and migration guide. https://geopandas.org/en/stable/docs/user_guide.html

Shapely 2.0 migration guide. https://shapely.readthedocs.io/en/stable/migration.html

Tenkanen, H., Heikinheimo, V., and Whipp, D. (forthcoming). Introduction to Python for Geographic Data Analysis. https://pythongis.org/

Dorman, M., Graser, A., Nowosad, J., and Lovelace, R. (2025). Geocomputation with Python. https://py.geocompx.org/

PySAL Spatial Statistics

Spatial weights, autocorrelation, LISA cluster analysis, spatial regression, and point pattern analysis using the PySAL ecosystem. For methodology and interpretation guidance, see the data-scientist skill's geospatial-analysis.md and geospatial-operations.md.

---

PySAL Ecosystem Overview

PySAL (Python Spatial Analysis Library) is a federation of packages:

Install all at once:

pip install pysal

---

Spatial Weights

Spatial weights formalize the concept of "neighbor" — they define which observations are connected and how strongly. Every spatial statistic depends on this choice.

Contiguity Weights (Polygons)

from libpysal.weights import Queen, Rook

# Queen: neighbors share an edge or vertex
w = Queen.from_dataframe(gdf)

# Rook: neighbors share an edge only (stricter)
w = Rook.from_dataframe(gdf)

# Inspect
print(w.n)                    # Number of observations
print(w.mean_neighbors)       # Average number of neighbors
print(w.min_neighbors)        # Minimum (watch for islands with 0)
print(w.max_neighbors)        # Maximum
print(w.islands)              # Observations with no neighbors
print(w.histogram)            # Distribution of neighbor counts

Distance-Based Weights (Points or Polygons)

from libpysal.weights import KNN, DistanceBand, Kernel

# K-Nearest Neighbors (every observation gets exactly k neighbors)
w = KNN.from_dataframe(gdf, k=6)

# Distance band (all neighbors within threshold distance)
w = DistanceBand.from_dataframe(gdf, threshold=10000)  # 10 km (projected CRS!)

# Kernel weights (distance-weighted, continuous)
w = Kernel.from_dataframe(gdf, bandwidth=15000)

CRS requirement: Distance-based weights compute distances between observations. The GeoDataFrame must be in a projected CRS (meters) — geographic CRS (degrees) produces meaningless distances.

Other Weight Types

from libpysal.weights import block_weights

# Block weights (same group = neighbor)
w = block_weights(gdf["state_fips"])

# Higher-order contiguity (neighbors of neighbors)
from libpysal.weights import higher_order
w2 = higher_order(w, k=2)  # 2nd-order neighbors

Graph API (Newer, Recommended)

libpysal's Graph class is the modern interface, backed by pandas/sparse matrices:

from libpysal.graph import Graph

# Contiguity
g = Graph.build_contiguity(gdf.geometry, rook=False)  # Queen (rook=False)
g = Graph.build_contiguity(gdf.geometry, rook=True)    # Rook

# KNN
g = Graph.build_knn(gdf.geometry, k=6)

# Distance band
g = Graph.build_distance_band(gdf.geometry, threshold=10000)

# Kernel
g = Graph.build_kernel(gdf.geometry, kernel="gaussian", k=10)

# Inspect
print(g.n)                  # Number of observations
print(g.n_edges)            # Number of neighbor pairs
print(g.cardinalities)      # Neighbor count per observation
print(g.isolates)           # Observations with no neighbors

# Convert to W for use with esda/spreg
w = g.to_W()

# Spatial lag (weighted average of neighbors' values)
spatial_lag = g.lag(gdf["poverty_rate"])

Row Standardization

Most applications use row-standardized weights (each row sums to 1), so the spatial lag is a weighted average:

# W objects
w.transform = "r"  # Row-standardize

# Graph objects
g_std = g.transform("R")

Saving and Loading Weights

# Save as parquet (Graph API)
g.to_parquet("weights.parquet")
g = Graph.read_parquet("weights.parquet")

# Save as GAL/GWT (legacy W format)
from libpysal.io import open as ps_open
gal = ps_open("weights.gal", "w")
gal.write(w)
gal.close()

---

Global Spatial Autocorrelation

Moran's I

Tests whether a variable is spatially clustered (positive I), dispersed (negative I), or random (I ≈ 0).

from esda.moran import Moran

# Compute Moran's I
mi = Moran(gdf["poverty_rate"], w, permutations=999)

# Results
print(f"Moran's I: {mi.I:.4f}")
print(f"Expected I: {mi.EI:.4f}")
print(f"p-value (permutation): {mi.p_sim:.4f}")
print(f"p-value (analytical): {mi.p_norm:.4f}")
print(f"z-score: {mi.z_sim:.4f}")

Geary's C

Complementary to Moran's I — more sensitive to local differences:

from esda.geary import Geary

gc = Geary(gdf["poverty_rate"], w, permutations=999)
print(f"Geary's C: {gc.C:.4f}")  # C < 1: positive autocorrelation, C > 1: negative
print(f"p-value: {gc.p_sim:.4f}")

Getis-Ord G

Tests for clustering of high values (hot spots) vs low values (cold spots):

from esda.getisord import G

go = G(gdf["poverty_rate"], w, permutations=999)
print(f"G: {go.G:.4f}")
print(f"p-value: {go.p_sim:.4f}")

Join Count Statistics

For binary variables (e.g., urban/rural, treatment/control), join counts test whether like values cluster:

from esda.join_counts import Join_Counts

# Binary variable (1 = urban, 0 = rural)
jc = Join_Counts(gdf["urban"].values, w, permutations=999)

print(f"BB joins (both 1): {jc.bb}")     # Black-black joins
print(f"WW joins (both 0): {jc.ww}")     # White-white joins
print(f"BW joins (mixed): {jc.bw}")      # Black-white joins
print(f"p-value (BB): {jc.p_sim_bb:.4f}")  # Clustering of 1s

---

Local Spatial Autocorrelation (LISA)

Local Moran's I

Identifies local clusters and outliers — where spatial autocorrelation is strongest.

from esda.moran import Moran_Local

# Compute LISA
lisa = Moran_Local(gdf["poverty_rate"], w, permutations=999)

# Results per observation
gdf["lisa_I"] = lisa.Is              # Local Moran's I value
gdf["lisa_q"] = lisa.q               # Quadrant (1=HH, 2=LH, 3=LL, 4=HL)
gdf["lisa_p"] = lisa.p_sim           # p-value (permutation)
gdf["lisa_sig"] = lisa.p_sim < 0.05  # Significant at 0.05

# Cluster labels
quadrant_labels = {1: "HH", 2: "LH", 3: "LL", 4: "HL"}
gdf["lisa_cluster"] = gdf["lisa_q"].map(quadrant_labels)
gdf.loc[~gdf["lisa_sig"], "lisa_cluster"] = "Not Significant"

LISA Quadrants

LISA Cluster Map

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# Standard LISA color scheme
lisa_colors = {
    "HH": "#d7191c",          # Red (hot spot)
    "LL": "#2c7bb6",          # Blue (cold spot)
    "HL": "#fdae61",          # Orange (high-low outlier)
    "LH": "#abd9e9",          # Light blue (low-high outlier)
    "Not Significant": "#f0f0f0"  # Light gray
}

fig, ax = plt.subplots(figsize=(12, 8))
gdf.plot(
    color=[lisa_colors[c] for c in gdf["lisa_cluster"]],
    edgecolor="white",
    linewidth=0.3,
    ax=ax
)

# Legend
patches = [mpatches.Patch(color=c, label=l) for l, c in lisa_colors.items() if l != "Not Significant"]
patches.append(mpatches.Patch(color=lisa_colors["Not Significant"], label="Not Significant"))
ax.legend(handles=patches, loc="lower right", fontsize=9)
ax.set_title("LISA Cluster Map: Poverty Rate", fontsize=14)
ax.set_axis_off()
plt.tight_layout()
plt.savefig("lisa_clusters.png", dpi=300, bbox_inches="tight")

Local Getis-Ord Gi*

Identifies statistically significant hot spots and cold spots:

from esda.getisord import G_Local

gi = G_Local(gdf["poverty_rate"], w, permutations=999)

gdf["gi_z"] = gi.Zs                # Z-scores
gdf["gi_p"] = gi.p_sim             # p-values
gdf["hotspot"] = "Not Significant"
gdf.loc[(gi.Zs > 0) & (gi.p_sim < 0.05), "hotspot"] = "Hot Spot"
gdf.loc[(gi.Zs < 0) & (gi.p_sim < 0.05), "hotspot"] = "Cold Spot"

---

Spatial Regression

OLS with Spatial Diagnostics

spreg's OLS includes Lagrange Multiplier tests that help choose the right spatial model:

import numpy as np
from spreg import OLS

# Prepare arrays
y = gdf[["poverty_rate"]].values           # Dependent variable (n, 1)
X = gdf[["median_income", "pct_rural"]].values  # Independent variables (n, k)

# OLS with spatial diagnostics
ols = OLS(
    y, X, w=w,
    name_y="poverty_rate",
    name_x=["median_income", "pct_rural"],
    name_ds="counties",
    spat_diag=True       # Include LM tests for spatial dependence
)

print(ols.summary)
# Look for:
# - Moran's I on residuals (significant = spatial dependence)
# - LM-lag and LM-error tests → guide model choice
# - Robust LM tests → which specification is preferred

LM Test Decision Rule

LM test results from OLS:
├─ LM-lag significant, LM-error not → Spatial Lag Model (SAR)
├─ LM-error significant, LM-lag not → Spatial Error Model (SEM)
├─ Both significant → Check robust versions:
│   ├─ Robust LM-lag significant → SAR
│   ├─ Robust LM-error significant → SEM
│   └─ Both robust significant → Spatial Durbin Model or SAR+SEM combo
└─ Neither significant → OLS is fine, no spatial model needed

Spatial Lag Model (SAR)

The outcome is influenced by neighbors' outcomes (Wy):

from spreg import ML_Lag, GM_Lag

# Maximum Likelihood estimation
sar_ml = ML_Lag(
    y, X, w=w,
    name_y="poverty_rate",
    name_x=["median_income", "pct_rural"],
    name_ds="counties"
)
print(sar_ml.summary)
# Key output: rho (spatial autoregressive coefficient), betas, log-likelihood

# GMM estimation (S2SLS — more robust to misspecification)
sar_gm = GM_Lag(
    y, X, w=w,
    name_y="poverty_rate",
    name_x=["median_income", "pct_rural"]
)

Spatial Error Model (SEM)

Spatial dependence is in the error term (unobserved spatially correlated factors):

from spreg import ML_Error, GM_Error

# Maximum Likelihood
sem_ml = ML_Error(
    y, X, w=w,
    name_y="poverty_rate",
    name_x=["median_income", "pct_rural"],
    name_ds="counties"
)
print(sem_ml.summary)
# Key output: lambda (spatial error coefficient), betas

# GMM (robust to heteroskedasticity)
sem_gm = GM_Error(
    y, X, w=w,
    name_y="poverty_rate",
    name_x=["median_income", "pct_rural"]
)

Spatial Durbin Model

Includes both spatial lag of Y and spatial lags of X (most flexible):

from spreg import ML_Lag

# Spatial Durbin = ML_Lag with slx_lags
sdm = ML_Lag(
    y, X, w=w,
    slx_lags=1,     # Include WX terms
    name_y="poverty_rate",
    name_x=["median_income", "pct_rural"],
    name_ds="counties"
)
print(sdm.summary)

Residual Diagnostics

After fitting a spatial model, verify that residuals no longer exhibit spatial autocorrelation:

from esda.moran import Moran

# Check residuals
mi_resid = Moran(model.u, w, permutations=999)  # model.u = residuals
print(f"Moran's I on residuals: {mi_resid.I:.4f}, p={mi_resid.p_sim:.4f}")
# p > 0.05 indicates spatial dependence has been adequately modeled

---

Point Pattern Analysis

Centrography (Descriptive Statistics)

from pointpats import centrography

points = np.column_stack([gdf.geometry.x, gdf.geometry.y])

# Mean center
mc = centrography.mean_center(points)

# Weighted mean center
wmc = centrography.weighted_mean_center(points, gdf["enrollment"].values)

# Standard distance (spatial spread)
sd = centrography.std_distance(points)

# Standard deviational ellipse
sx, sy, theta = centrography.ellipse(points)

Quadrat Analysis

Test whether points are uniformly distributed:

from pointpats import QStatistic

q = QStatistic(points, nx=10, ny=10)  # 10x10 grid
print(f"Chi-squared: {q.chi2:.2f}, p-value: {q.chi2_pvalue:.4f}")
# Significant p → points are not uniformly distributed

Kernel Density Estimation (KDE)

import numpy as np
from scipy.stats import gaussian_kde

# Compute KDE
xy = np.vstack([gdf.geometry.x, gdf.geometry.y])
kde = gaussian_kde(xy)

# Evaluate on grid
xmin, ymin, xmax, ymax = gdf.total_bounds
xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([xx.ravel(), yy.ravel()])
density = kde(positions).reshape(xx.shape)

# Plot
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10, 8))
ax.contourf(xx, yy, density, cmap="YlOrRd", levels=20)
gdf.plot(ax=ax, color="black", markersize=1, alpha=0.3)
ax.set_title("Point Density (KDE)")

Ripley's Functions

Ripley's functions characterize point patterns at multiple scales — testing whether points are clustered, dispersed, or random at different distances.

from pointpats import PointPattern, Genv, Fenv, Kenv

# Create point pattern from coordinates
pp = PointPattern(np.column_stack([gdf.geometry.x, gdf.geometry.y]))

# Ripley's G (nearest-neighbor distance distribution)
# Tests if points are closer together than expected under randomness
g = Genv(pp, intervals=20, realizations=99)
# g.observed  = observed G function
# g.low, g.high = simulation envelope (under CSR)
# If observed > high → clustering at that distance

# Ripley's F (empty space function)
# Tests if empty spaces are smaller than expected (= clustering)
f = Fenv(pp, intervals=20, realizations=99)

# Ripley's K (cumulative neighbor count by distance)
# Most commonly used; L function = variance-stabilized K
k = Kenv(pp, intervals=20, realizations=99)

If the observed function falls outside the simulation envelope, the point pattern deviates significantly from Complete Spatial Randomness (CSR) at that distance. For interpretation guidance, see the data-scientist skill's geospatial-analysis.md.

---

Bivariate Spatial Autocorrelation

Test whether the spatial pattern of one variable is related to the spatial pattern of another:

from esda.moran import Moran_BV

# Bivariate Moran's I: is poverty spatially correlated with unemployment?
bv = Moran_BV(
    gdf["poverty_rate"].values,
    gdf["unemployment_rate"].values,
    w,
    permutations=999
)
print(f"Bivariate Moran's I: {bv.I:.4f}, p={bv.p_sim:.4f}")

---

Rate-Adjusted Statistics

When analyzing rates derived from counts (e.g., crime rates, disease rates), use rate-adjusted versions to account for variance instability in small populations:

from esda.moran import Moran_Rate, Moran_Local_Rate

# Rate-adjusted global Moran's I
mr = Moran_Rate(
    gdf["crime_count"].values,    # Event count (numerator)
    gdf["population"].values,      # Population (denominator)
    w,
    permutations=999
)

# Rate-adjusted LISA
mlr = Moran_Local_Rate(
    gdf["crime_count"].values,
    gdf["population"].values,
    w,
    permutations=999
)

---

References and Further Reading

Rey, S.J., Arribas-Bel, D., and Wolf, L.J. (2023). Geographic Data Science with Python. CRC Press. https://geographicdata.science/book/

• Ch. 4: Spatial weights
• Ch. 6: Global spatial autocorrelation
• Ch. 7: Local spatial autocorrelation
• Ch. 8: Point pattern analysis
• Ch. 11: Spatial regression
• Ch. 12: Spatial feature engineering

Anselin, L. (1995). "Local Indicators of Spatial Association — LISA." Geographical Analysis, 27(2), 93-115.

Anselin, L. and Rey, S.J. (2014). Modern Spatial Econometrics in Practice. GeoDa Press.

Rey, S.J. et al. (2022). "The PySAL Ecosystem: Philosophy and Implementation." Geographical Analysis, 54(3), 467-487. https://pysal.org/

PySAL API documentation:

• libpysal: https://pysal.org/libpysal/
• esda: https://pysal.org/esda/
• spreg: https://pysal.org/spreg/
• pointpats: https://pysal.org/pointpats/

GeoPandas Quickstart

Installation

Basic Install

pip install geopandas
# or
conda install -c conda-forge geopandas

geopandas 1.x automatically installs core dependencies: shapely (>=2.0), pyproj (>=3.3), pyogrio (>=0.7.2, default I/O engine), and pandas.

Recommended Extras

# Visualization and basemaps
pip install matplotlib contextily mapclassify folium

# Raster integration
pip install rasterio rasterstats rioxarray

# Spatial statistics (PySAL ecosystem)
pip install libpysal esda spreg pointpats

# GPU-accelerated visualization for large datasets
pip install lonboard

# Interactive notebooks
pip install mapwidget ipywidgets

Verify Installation

import geopandas as gpd
print(gpd.__version__)       # Should be 1.x
gpd.show_versions()          # Full dependency report (prints directly)

---

Core Concepts

GeoDataFrame = DataFrame + Geometry

A GeoDataFrame is a pandas DataFrame with a special geometry column containing Shapely geometry objects. All standard pandas operations work — plus spatial operations.

import geopandas as gpd
import pandas as pd
from shapely.geometry import Point

# A GeoDataFrame has:
# - Regular columns (name, population, etc.)
# - A geometry column (points, lines, or polygons)
# - A CRS (coordinate reference system)

GeoSeries

The geometry column is a GeoSeries — a pandas Series of Shapely geometries. It provides vectorized spatial operations:

gdf.geometry              # Access the GeoSeries
gdf.geometry.area         # Area of each geometry
gdf.geometry.centroid     # Centroid of each geometry
gdf.geometry.is_valid     # Validity check for each geometry

---

Creating GeoDataFrames

From a DataFrame with Coordinates

The most common case — a DataFrame with latitude/longitude columns:

import geopandas as gpd
import pandas as pd

df = pd.DataFrame({
    "name": ["School A", "School B", "School C"],
    "lon": [-77.036, -77.009, -76.995],
    "lat": [38.901, 38.889, 38.910],
    "enrollment": [500, 800, 650]
})

gdf = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy(df.lon, df.lat),
    crs="EPSG:4326"
)

gpd.points_from_xy(x, y) takes longitude as x and latitude as y. This is the standard mathematical convention (x=horizontal=longitude) but trips up many users who think "lat, lon". See gotchas.md for more on coordinate order.

From Shapely Geometry Objects

from shapely.geometry import Point, Polygon

gdf = gpd.GeoDataFrame(
    {"name": ["Park", "Lake"], "type": ["green", "water"]},
    geometry=[
        Polygon([(-77.05, 38.90), (-77.04, 38.90), (-77.04, 38.91), (-77.05, 38.91)]),
        Polygon([(-77.02, 38.88), (-77.01, 38.88), (-77.01, 38.89), (-77.02, 38.89)])
    ],
    crs="EPSG:4326"
)

From a GeoJSON-like Dictionary

gdf = gpd.GeoDataFrame.from_features([
    {
        "type": "Feature",
        "geometry": {"type": "Point", "coordinates": [-77.036, 38.901]},
        "properties": {"name": "Capitol", "visitors": 3000000}
    }
], crs="EPSG:4326")

Empty GeoDataFrame with Schema

from shapely.geometry import Point

gdf = gpd.GeoDataFrame(
    columns=["name", "value", "geometry"],
    geometry="geometry",
    crs="EPSG:4326"
)

---

Reading Spatial Data

Read Any Supported Vector Format

# GeoPackage (recommended modern format)
gdf = gpd.read_file("counties.gpkg")

# Shapefile (legacy — see data-io.md for limitations)
gdf = gpd.read_file("counties.shp")

# GeoJSON
gdf = gpd.read_file("counties.geojson")

# FlatGeobuf
gdf = gpd.read_file("counties.fgb")

Read with Filters (Efficient for Large Files)

# Read only features within a bounding box (west, south, east, north)
gdf = gpd.read_file("large_file.gpkg", bbox=(-78.0, 38.0, -76.0, 40.0))

# Read only specific columns
gdf = gpd.read_file("large_file.gpkg", columns=["NAME", "POP", "geometry"])

# Read limited rows
gdf = gpd.read_file("large_file.gpkg", rows=100)

# Read specific layer from multi-layer file
gdf = gpd.read_file("multi_layer.gpkg", layer="counties")

Read GeoParquet (Fastest for Analytical Workflows)

gdf = gpd.read_parquet("counties.parquet")

# With bounding box filter
gdf = gpd.read_parquet("counties.parquet", bbox=(-78.0, 38.0, -76.0, 40.0))

# With column selection
gdf = gpd.read_parquet("counties.parquet", columns=["NAME", "POP", "geometry"])

For more formats and advanced I/O, see data-io.md.

---

Writing Spatial Data

# GeoPackage (preferred)
gdf.to_file("output.gpkg", driver="GPKG")

# GeoParquet (preferred for analytical pipelines)
gdf.to_parquet("output.parquet")

# Shapefile (legacy — avoid if possible)
gdf.to_file("output.shp")

# GeoJSON
gdf.to_file("output.geojson", driver="GeoJSON")

---

Basic Inspection

gdf.head()                # First rows
gdf.shape                 # (rows, columns)
gdf.columns               # Column names
gdf.dtypes                # Column types (geometry shows as 'geometry')
gdf.crs                   # Coordinate Reference System
gdf.total_bounds          # [minx, miny, maxx, maxy] bounding box
gdf.geom_type.unique()    # Geometry types present (Point, Polygon, etc.)
gdf.geometry.is_valid.all()  # Check all geometries are valid

---

Basic Plotting

Quick Plot

# Default plot — just the geometries
gdf.plot()

# Choropleth — color by a column
gdf.plot(column="population", legend=True)

# With figure size
gdf.plot(column="population", legend=True, figsize=(12, 8))

Layered Plot

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(12, 8))
counties.plot(ax=ax, color="lightgray", edgecolor="white")
schools.plot(ax=ax, color="red", markersize=5)
ax.set_title("Schools by County")
ax.set_axis_off()
plt.tight_layout()
plt.savefig("map.png", dpi=300, bbox_inches="tight")

Interactive Map

# Built-in folium integration (opens in notebook or saves to HTML)
gdf.explore(column="population", cmap="YlOrRd", legend=True)

# Save to HTML
m = gdf.explore(column="population")
m.save("map.html")

For advanced visualization (basemaps, classification schemes, lonboard), see visualization.md.

---

Essential Spatial Operations Preview

# Reproject
gdf_projected = gdf.to_crs(epsg=5070)

# Spatial join (which polygon contains each point?)
result = gpd.sjoin(points_gdf, polygons_gdf, predicate="within")

# Buffer (create 1km buffer around points — CRS must be in meters)
gdf_buffered = gdf_projected.copy()
gdf_buffered["geometry"] = gdf_projected.buffer(1000)

# Dissolve (merge counties into states)
states = counties.dissolve(by="state_fips", aggfunc="sum")

# Overlay (intersect two polygon layers)
overlap = gpd.overlay(layer1, layer2, how="intersection")

For complete spatial operation reference, see spatial-operations.md.

---

Next Steps

• Learn about CRS and projections — essential before any geometric computation
• Master spatial operations — joins, overlays, dissolve
• Explore visualization options — static and interactive maps
• Understand file formats — choosing the right format for your workflow

References and Further Reading

Jordahl, K. et al. (2024). geopandas: Python tools for geographic data. https://geopandas.org/

Tenkanen, H., Heikinheimo, V., and Whipp, D. (2024). Introduction to Python for Geographic Data Analysis. CRC Press. https://pythongis.org/

Dorman, M., Graser, A., Nowosad, J., and Lovelace, R. (2025). Geocomputation with Python. CRC Press. https://py.geocompx.org/

Raster Integration

Working with raster data alongside geopandas — reading GeoTIFFs with rasterio, multidimensional rasters with xarray/rioxarray, zonal statistics, and raster-vector conversion.

---

Rasterio Basics

rasterio provides Python access to GDAL's raster I/O. It reads and writes GeoTIFF and other raster formats.

Reading a GeoTIFF

import rasterio
import numpy as np

with rasterio.open("elevation.tif") as src:
    # Metadata
    print(src.crs)          # CRS
    print(src.bounds)       # BoundingBox(left, bottom, right, top)
    print(src.res)          # (pixel_width, pixel_height)
    print(src.shape)        # (rows, cols)
    print(src.count)        # Number of bands
    print(src.dtypes)       # Data type per band
    print(src.nodata)       # NoData value
    print(src.transform)    # Affine transform (pixel→world coordinates)

    # Read band data as NumPy array
    band1 = src.read(1)     # Read band 1 (1-indexed)
    all_bands = src.read()  # Read all bands: shape (bands, rows, cols)

Reading a Window (Subset)

from rasterio.windows import from_bounds

with rasterio.open("large_raster.tif") as src:
    # Read only a geographic extent
    window = from_bounds(
        left=-78.0, bottom=38.0, right=-76.0, top=40.0,
        transform=src.transform
    )
    subset = src.read(1, window=window)

Writing a GeoTIFF

import rasterio
from rasterio.transform import from_bounds
import numpy as np

data = np.random.rand(100, 100).astype(np.float32)

with rasterio.open(
    "output.tif",
    mode="w",
    driver="GTiff",
    height=data.shape[0],
    width=data.shape[1],
    count=1,                    # Number of bands
    dtype=data.dtype,
    crs="EPSG:4326",
    transform=from_bounds(-78, 38, -76, 40, data.shape[1], data.shape[0]),
    nodata=-9999
) as dst:
    dst.write(data, 1)         # Write to band 1

Masking Raster by Polygon

import rasterio
from rasterio.mask import mask
import geopandas as gpd

gdf = gpd.read_file("study_area.gpkg")
geometries = gdf.geometry.values

with rasterio.open("raster.tif") as src:
    # Clip raster to polygon boundaries
    out_image, out_transform = mask(
        src,
        geometries,
        crop=True,        # Crop to geometry extent
        nodata=-9999,     # Fill outside with NoData
        all_touched=False # Only cells whose center is inside
    )
    out_meta = src.meta.copy()
    out_meta.update({
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform
    })

# Write clipped raster
with rasterio.open("clipped.tif", "w", **out_meta) as dst:
    dst.write(out_image)

Extracting Values at Points

import rasterio
import geopandas as gpd

points = gpd.read_file("sample_points.gpkg")

with rasterio.open("temperature.tif") as src:
    # Ensure same CRS
    points_proj = points.to_crs(src.crs)

    # Sample raster at each point location
    coords = [(pt.x, pt.y) for pt in points_proj.geometry]
    values = [val[0] for val in src.sample(coords)]

    points["temperature"] = values

---

xarray + rioxarray

For multidimensional rasters (multiple bands, time series, or named dimensions), xarray with the rioxarray extension provides a higher-level interface.

Installation

pip install xarray rioxarray netcdf4

Reading Rasters

import xarray as xr
import rioxarray  # Extends xarray with .rio accessor

# Read a GeoTIFF
ds = xr.open_dataarray("temperature.tif", engine="rasterio")
# or
ds = rioxarray.open_rasterio("temperature.tif")

# Inspect
print(ds.dims)       # ('band', 'y', 'x')
print(ds.rio.crs)    # CRS
print(ds.rio.bounds())  # Bounding box
print(ds.rio.resolution())  # Pixel size

CRS Operations

# Reproject
ds_proj = ds.rio.reproject("EPSG:5070")

# Set CRS (if missing)
ds = ds.rio.write_crs("EPSG:4326")

Clipping by Polygon

import geopandas as gpd

gdf = gpd.read_file("study_area.gpkg")

# Clip raster to polygon
ds_clipped = ds.rio.clip(gdf.geometry, gdf.crs)

# Clip to bounding box
ds_bbox = ds.rio.clip_box(minx=-78, miny=38, maxx=-76, maxy=40)

Writing

ds.rio.to_raster("output.tif")

NetCDF / Multi-Temporal Rasters

# Read NetCDF with time dimension
ds = xr.open_dataset("climate_data.nc")

# Select a time slice
temp_jan = ds["temperature"].sel(time="2024-01")

# Spatial operations on each time step
ds_proj = ds.rio.reproject("EPSG:5070")

---

Zonal Statistics

Zonal statistics summarize raster cell values within vector polygon boundaries — for example, computing mean elevation within each county.

Using rasterstats

from rasterstats import zonal_stats
import geopandas as gpd

counties = gpd.read_file("counties.gpkg")

# Basic zonal statistics
stats = zonal_stats(
    counties,                    # Vector polygons (GeoDataFrame or path)
    "temperature.tif",           # Raster (path or numpy array)
    stats=["mean", "min", "max", "std", "count", "median"]
)
# Returns: list of dicts, one per polygon
# [{'mean': 15.2, 'min': 8.1, 'max': 22.4, ...}, ...]

# Merge back to GeoDataFrame
import pandas as pd
stats_df = pd.DataFrame(stats)
counties_with_stats = pd.concat([counties, stats_df], axis=1)

Zonal Statistics Parameters

zonal_stats(
    vectors,                # Polygons (GeoDataFrame, path, or GeoJSON)
    raster,                 # Raster (path or ndarray + affine transform)
    stats=["mean", "sum"],  # Statistics to compute
    all_touched=False,      # Include cells that touch boundary (not just center-in)
    nodata=None,            # Override raster NoData value
    categorical=False,      # For categorical rasters: count per category
    category_map=None,      # Map raster values to labels
    band=1,                 # Which raster band
    geojson_out=False,      # Return GeoJSON features
    prefix=""               # Prefix for output column names
)

Available Statistics

Choosing the Right Statistic

Categorical Zonal Statistics

# Count pixels per land use class within each polygon
stats = zonal_stats(
    counties,
    "landuse.tif",
    categorical=True,
    category_map={1: "urban", 2: "forest", 3: "water", 4: "agriculture"}
)
# Returns: [{'urban': 450, 'forest': 1200, 'water': 50, 'agriculture': 800}, ...]

Data Completeness Check

Always report how many cells contributed to each polygon's statistics:

stats = zonal_stats(counties, "raster.tif", stats=["mean", "count", "nodata"])
stats_df = pd.DataFrame(stats)

# Flag polygons with too few cells (unreliable statistics)
stats_df["reliable"] = stats_df["count"] >= 10
print(f"Unreliable polygons: {(~stats_df['reliable']).sum()}")

---

Raster-Vector Conversion

Rasterize (Vector to Raster)

Convert vector features to a raster grid:

from rasterio.features import rasterize
from rasterio.transform import from_bounds
import numpy as np

# Define output grid
transform = from_bounds(*gdf.total_bounds, width=1000, height=1000)

# Burn values from a column
shapes = [(geom, val) for geom, val in zip(gdf.geometry, gdf["population"])]
raster = rasterize(
    shapes,
    out_shape=(1000, 1000),
    transform=transform,
    fill=0,                  # Background value
    dtype=np.float32,
    all_touched=False
)

Vectorize (Raster to Vector)

Convert raster cells to vector polygons:

from rasterio.features import shapes
import geopandas as gpd
from shapely.geometry import shape

with rasterio.open("classified.tif") as src:
    image = src.read(1)
    mask = image != src.nodata

    results = [
        {"geometry": shape(geom), "value": val}
        for geom, val in shapes(image, mask=mask, transform=src.transform)
    ]

gdf = gpd.GeoDataFrame(results, crs=src.crs)

---

CRS Matching Between Raster and Vector

Raster and vector data must be in the same CRS for zonal statistics, masking, and extraction. Always verify:

import rasterio
import geopandas as gpd

gdf = gpd.read_file("polygons.gpkg")

with rasterio.open("raster.tif") as src:
    if gdf.crs != src.crs:
        gdf = gdf.to_crs(src.crs)
    # Now safe to use together

---

References and Further Reading

Gillies, S. et al. (2024). rasterio: Fast and direct raster I/O for Python. https://rasterio.readthedocs.io/

Hoyer, S. and Hamman, J. (2017). "xarray: N-D labeled arrays and datasets in Python." Journal of Open Research Software, 5(1), 10.

Snow, A. et al. (2024). rioxarray: xarray extension for rasterio. https://corteva.github.io/rioxarray/

Perry, M. (2024). rasterstats: Summary statistics of geospatial raster datasets based on vector geometries. https://pythonhosted.org/rasterstats/

Dorman, M., Graser, A., Nowosad, J., and Lovelace, R. (2025). Geocomputation with Python, Chs. 5-6: "Raster-vector interactions" and "Reprojecting geographic data." https://py.geocompx.org/

Spatial Operations

Vector spatial operations in geopandas — joins, overlays, dissolve, clip, buffer, distance, and areal interpolation. For methodology guidance (when/why to use each operation), see the data-scientist skill's geospatial-operations.md.

---

Spatial Joins

Spatial joins connect records from two GeoDataFrames based on geographic relationships rather than shared key columns.

Basic Spatial Join

# Which county contains each school? (point-in-polygon)
result = gpd.sjoin(schools, counties, predicate="within")

# Which schools fall within each county? (polygon-contains-point)
result = gpd.sjoin(counties, schools, predicate="contains")

# Which features overlap? (most permissive)
result = gpd.sjoin(gdf1, gdf2, predicate="intersects")

Spatial Join Parameters

gpd.sjoin(
    left_df,                # Left GeoDataFrame
    right_df,               # Right GeoDataFrame
    how="inner",            # 'inner', 'left', 'right'
    predicate="intersects", # Spatial predicate
    lsuffix="left",         # Suffix for overlapping column names from left
    rsuffix="right"         # Suffix for overlapping column names from right
)

Spatial Predicates

Nearest-Neighbor Join

# Find nearest school to each census tract centroid
result = gpd.sjoin_nearest(
    tracts,                   # Left GeoDataFrame
    schools,                  # Right GeoDataFrame
    how="left",               # Keep all left features
    max_distance=10000,       # Max search distance (CRS units)
    distance_col="dist_m"     # Add column with actual distance
)

sjoin_nearest requires both GeoDataFrames in a projected CRS for meaningful distances.

Post-Join Validation

Spatial joins can produce unexpected row counts due to many-to-many relationships. Always validate:

print(f"Left rows: {len(schools)}")
print(f"Result rows: {len(result)}")
print(f"Duplicated left indices: {result.index.duplicated().sum()}")
print(f"Null join columns: {result['county_name'].isna().sum()}")

# If duplicates exist, decide how to handle:
# Option 1: Keep first match
result_dedup = result[~result.index.duplicated(keep="first")]

# Option 2: Aggregate
result_agg = result.groupby(result.index).agg({"county_name": "first", "value": "sum"})

---

Attribute Joins (Non-Spatial)

Standard pandas merge for joining by shared columns:

# Join census data to county geometries by FIPS code
counties_with_data = counties.merge(census_df, left_on="GEOID", right_on="fips", how="left")

# Note: merge returns a GeoDataFrame if the left input is a GeoDataFrame

---

Overlay Operations

Overlays combine two polygon layers to produce new geometries. Unlike spatial joins (which transfer attributes), overlays create new geometries from intersections.

gpd.overlay(
    df1,                   # First GeoDataFrame (polygons)
    df2,                   # Second GeoDataFrame (polygons)
    how="intersection",    # 'intersection', 'union', 'difference',
                           # 'symmetric_difference', 'identity'
    keep_geom_type=True,   # Filter out geometry type changes
    make_valid=True        # Auto-fix invalid geometries before overlay
)

Overlay Types

Overlay Example

# Find areas where flood zones and school districts overlap
flood_school = gpd.overlay(school_districts, flood_zones, how="intersection")

# Compute area of overlap
flood_school["overlap_area_km2"] = flood_school.to_crs(epsg=5070).area / 1e6

---

Dissolve (Aggregate Geometries)

Dissolve merges geometries by a grouping column, combining multiple features into one. It is the spatial equivalent of groupby().agg().

# Merge counties into states (sum population, merge geometries)
states = counties.dissolve(by="state_fips", aggfunc="sum")

# Multiple aggregation functions
states = counties.dissolve(
    by="state_fips",
    aggfunc={
        "population": "sum",
        "area_sq_mi": "sum",
        "median_income": "mean"
    }
)

# Dissolve all features into one (no grouping)
us_boundary = counties.dissolve()

Dissolve Parameters

gdf.dissolve(
    by=None,              # Column(s) to group by (None = dissolve all)
    aggfunc="first",      # Aggregation: 'first', 'sum', 'mean', 'min', 'max', or dict
    as_index=True,        # Use group column as index
    sort=True,            # Sort by group column
    method="unary"        # 'unary' (default), 'coverage_union' (faster if no overlaps)
)

---

Clipping

Clip features to a boundary — everything outside the mask is removed.

# Clip schools to a state boundary
schools_in_state = gpd.clip(schools, state_boundary)

# Clip a polygon layer to a bounding box
from shapely.geometry import box
bbox = box(-78.0, 38.0, -76.0, 40.0)
clipped = gpd.clip(gdf, bbox)

clip vs sjoin

• `clip`: Modifies geometries — polygons are cut at the mask boundary
• `sjoin`: Keeps original geometries — only filters which features to include

---

Buffering

Create a zone around features at a specified distance.

# Buffer points by 1 km (CRS must be in meters!)
gdf_proj = gdf.to_crs(epsg=5070)
gdf_buffered = gdf_proj.copy()
gdf_buffered["geometry"] = gdf_proj.buffer(1000)  # 1000 meters

# Variable-distance buffer (different distance per feature)
gdf_buffered["geometry"] = gdf_proj.buffer(gdf_proj["radius_m"])

# Buffer parameters
gdf_proj.buffer(
    distance=1000,       # Distance in CRS units
    resolution=16,       # Number of segments per quarter circle
    cap_style="round",   # 'round', 'flat', 'square'
    join_style="round",  # 'round', 'mitre', 'bevel'
    single_sided=False   # True for one-sided buffer (lines only)
)

Buffering in a geographic CRS (degrees) produces elliptical buffers that vary with latitude — always project first.

---

Distance

# Distance between aligned GeoSeries (element-wise)
distances = gdf1.distance(gdf2)  # Returns Series of distances in CRS units

# Distance from every feature to a single point
from shapely.geometry import Point
capitol = Point(-77.009, 38.890)
gdf["dist_to_capitol"] = gdf.to_crs(epsg=5070).distance(
    gpd.GeoSeries([capitol], crs="EPSG:4326").to_crs(epsg=5070).iloc[0]
)

Always compute distances in a projected CRS for results in meters.

---

Centroid

# Centroid of each geometry
gdf["centroid"] = gdf.to_crs(epsg=5070).centroid

# Representative point (guaranteed to be inside the polygon — centroid might not be)
gdf["rep_point"] = gdf.to_crs(epsg=5070).representative_point()

Use representative_point() for irregular or concave polygons where the centroid might fall outside the geometry.

---

Area and Length

# Area (project to equal-area CRS first)
gdf_proj = gdf.to_crs(epsg=5070)
gdf["area_m2"] = gdf_proj.area
gdf["area_km2"] = gdf_proj.area / 1e6

# Length (for LineString geometries)
gdf["length_m"] = gdf_proj.length
gdf["length_km"] = gdf_proj.length / 1e3

---

Geometry Manipulation

Simplify

Reduce geometry complexity (fewer vertices) for faster rendering or smaller file size:

# Simplify with tolerance (in CRS units)
gdf["geometry"] = gdf.simplify(tolerance=100)  # 100 meters if projected

# Preserve topology (prevents holes between adjacent polygons)
gdf["geometry"] = gdf.simplify(tolerance=100, preserve_topology=True)

Convex Hull

gdf["convex_hull"] = gdf.convex_hull

Explode (Multi to Single)

Split MultiPolygons/MultiLineStrings into individual geometries:

gdf_single = gdf.explode(index_parts=False)

Unary Union

Merge all geometries into one:

merged = gdf.geometry.union_all()  # Returns a single Shapely geometry

Bounds

# Bounding box of each geometry
bounds = gdf.bounds  # DataFrame with minx, miny, maxx, maxy columns

# Overall bounding box
total_bounds = gdf.total_bounds  # Array: [minx, miny, maxx, maxy]

---

Areal Interpolation (Boundary Mismatch)

When source and target boundaries don't align (e.g., redistributing census tract data to school districts), use the tobler package:

from tobler.area_weighted import area_interpolate

# Redistribute population from tracts to school districts
result = area_interpolate(
    source_df=tracts,                # Source polygons with data
    target_df=school_districts,      # Target polygons to receive data
    extensive_variables=["population", "housing_units"],  # Totals: distributed by area fraction
    intensive_variables=["median_income", "poverty_rate"]  # Rates: area-weighted average
)

Extensive vs intensive:

• Extensive (population, count): scales with area. 1000 people in a tract, 60% area overlap → 600 assigned.
• Intensive (rate, density): independent of area. 15% poverty rate → area-weighted average.

Confusing them is the most common areal interpolation error.

Dasymetric Refinement

Simple area-weighting assumes uniform distribution within source zones. Dasymetric methods use auxiliary data (land use rasters) for better accuracy:

from tobler.dasymetric import masked_area_interpolate

result = masked_area_interpolate(
    source_df=tracts,
    target_df=districts,
    extensive_variables=["population"],
    raster="landuse.tif",             # Binary raster: 1 = inhabited, 0 = uninhabited
)

---

References and Further Reading

Jordahl, K. et al. (2024). geopandas: Merging data. https://geopandas.org/en/stable/docs/user_guide/mergingdata.html

Jordahl, K. et al. (2024). geopandas: Set operations with overlay. https://geopandas.org/en/stable/docs/user_guide/set_operations.html

Rey, S.J., Arribas-Bel, D., and Wolf, L.J. (2023). Geographic Data Science with Python, Ch. 4: "Spatial weights." https://geographicdata.science/book/

Dorman, M., Graser, A., Nowosad, J., and Lovelace, R. (2025). Geocomputation with Python, Chs. 3-4: "Attribute data operations" and "Spatial data operations." https://py.geocompx.org/

tobler documentation. https://pysal.org/tobler/

Spatial Visualization

Making maps with geopandas and the Python visualization ecosystem — static choropleths, classification schemes, basemap tiles, interactive maps, and GPU-accelerated rendering for large datasets.

---

Static Maps with .plot()

geopandas .plot() wraps matplotlib for quick static maps.

Basic Choropleth

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(12, 8))
gdf.plot(
    column="poverty_rate",    # Column to color by
    ax=ax,
    legend=True,
    cmap="YlOrRd",            # Colormap
    edgecolor="white",
    linewidth=0.3,
    missing_kwds={"color": "lightgrey", "label": "No data"}
)
ax.set_title("Poverty Rate by County", fontsize=14)
ax.set_axis_off()
plt.tight_layout()
plt.savefig("choropleth.png", dpi=300, bbox_inches="tight")

Layered Map

fig, ax = plt.subplots(figsize=(12, 8))

# Base layer: counties
counties.plot(ax=ax, color="lightyellow", edgecolor="gray", linewidth=0.5)

# Overlay: schools colored by enrollment
schools.plot(
    ax=ax,
    column="enrollment",
    cmap="Blues",
    markersize=10,
    legend=True,
    legend_kwds={"label": "Enrollment", "shrink": 0.6}
)

# Overlay: district boundaries
districts.plot(ax=ax, facecolor="none", edgecolor="red", linewidth=1.5)

ax.set_title("Schools by Enrollment")
ax.set_axis_off()
plt.tight_layout()

Categorical Map

gdf.plot(
    column="school_type",     # Categorical column
    categorical=True,
    legend=True,
    cmap="Set2",
    legend_kwds={"loc": "lower right", "fontsize": 8, "title": "Type"}
)

Common Colormaps

Plot Parameters Quick Reference

---

Classification Schemes (mapclassify)

mapclassify provides statistical classification methods for choropleth maps. Without classification, continuous color ramps can obscure patterns.

Installation

pip install mapclassify

Using with .plot()

gdf.plot(
    column="poverty_rate",
    scheme="FisherJenks",   # Classification method
    k=5,                     # Number of classes
    cmap="YlOrRd",
    legend=True,
    legend_kwds={"loc": "lower right", "fontsize": 8}
)

Available Schemes

Direct mapclassify Usage

import mapclassify

# Compute classification
classifier = mapclassify.FisherJenks(gdf["poverty_rate"], k=5)
print(classifier.bins)        # Class breakpoints
print(classifier.counts)      # Observations per class
print(classifier.adcm)        # Absolute deviation around class medians

# Apply custom classification
gdf["class"] = mapclassify.UserDefined(
    gdf["poverty_rate"],
    bins=[5, 10, 15, 20, 30]
).yb  # .yb = class labels (0-indexed)

Choosing a Classification Scheme

What's your data distribution?
├─ Roughly uniform → EqualInterval
├─ Skewed (most values low, few high) → Quantiles or FisherJenks
├─ Heavy-tailed / power-law → HeadTailBreaks
├─ Need to show deviation from average → StdMean
├─ Domain-specific thresholds exist → UserDefined
└─ Unsure → FisherJenks (safest default)

---

Basemap Tiles (contextily)

contextily adds web map tiles (OpenStreetMap, CartoDB, etc.) as background layers.

Installation

pip install contextily

Adding Basemaps

import contextily as ctx
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(12, 8))
gdf.to_crs(epsg=3857).plot(ax=ax, alpha=0.5, edgecolor="red")

# Add basemap tiles
ctx.add_basemap(ax, source=ctx.providers.CartoDB.Positron)
ax.set_axis_off()

CRS requirement: Basemap tiles are in Web Mercator (EPSG:3857). Either reproject your data to 3857 before plotting, or pass the crs parameter:

# Option 1: Reproject data to 3857
gdf_3857 = gdf.to_crs(epsg=3857)
gdf_3857.plot(ax=ax)
ctx.add_basemap(ax)

# Option 2: Keep data in its CRS, let contextily handle reprojection
gdf.plot(ax=ax)
ctx.add_basemap(ax, crs=gdf.crs, source=ctx.providers.CartoDB.Positron)

Popular Tile Providers

Controlling Zoom Level

# Auto zoom (default)
ctx.add_basemap(ax, zoom="auto")

# Specific zoom level (higher = more detail, more tiles)
ctx.add_basemap(ax, zoom=10)

# Adjust auto zoom
ctx.add_basemap(ax, zoom="auto", zoom_adjust=1)  # One level more detail

---

Interactive Maps (folium via .explore())

geopandas .explore() creates folium-based interactive maps with pan, zoom, and hover tooltips.

Basic Interactive Map

# Quick interactive map
m = gdf.explore(
    column="poverty_rate",
    cmap="YlOrRd",
    legend=True,
    tooltip=["NAME", "poverty_rate", "population"],
    popup=True,
    tiles="CartoDB positron"
)

# Save to HTML
m.save("interactive_map.html")

Layered Interactive Map

# Start with one layer
m = counties.explore(
    column="poverty_rate",
    cmap="YlOrRd",
    name="Counties",
    legend=True,
    tooltip=["NAME", "poverty_rate"]
)

# Add more layers
schools.explore(
    m=m,                          # Pass existing map
    color="blue",
    marker_kwds={"radius": 3},
    name="Schools",
    tooltip=["school_name", "enrollment"]
)

# Add layer control
import folium
folium.LayerControl().add_to(m)

m.save("layered_map.html")

Direct Folium Usage

import folium

# Create map centered on data
center = [gdf.geometry.centroid.y.mean(), gdf.geometry.centroid.x.mean()]
m = folium.Map(location=center, zoom_start=10, tiles="CartoDB positron")

# Add choropleth
folium.Choropleth(
    geo_data=gdf.to_json(),
    data=gdf,
    columns=["GEOID", "poverty_rate"],
    key_on="feature.properties.GEOID",
    fill_color="YlOrRd",
    fill_opacity=0.7,
    line_opacity=0.2,
    legend_name="Poverty Rate (%)"
).add_to(m)

m.save("folium_map.html")

---

GPU-Accelerated Rendering (lonboard)

lonboard uses deck.gl for GPU-accelerated rendering of very large datasets (millions of features) that would overwhelm matplotlib or folium.

Installation

pip install lonboard

Quick Visualization

from lonboard import viz

# Automatic layer type detection based on geometry
m = viz(gdf)

# Save to HTML
m.to_html("large_map.html")

For column-based coloring, use the layer API (below) which gives explicit control over color mapping.

Layer-Based API

from lonboard import Map, ScatterplotLayer, PolygonLayer

# Points
layer = ScatterplotLayer.from_geopandas(
    schools_gdf,
    get_radius=500,           # Radius in meters
    radius_units="meters",
    get_fill_color=[255, 0, 0, 180],  # RGBA
    pickable=True
)

# Polygons
layer = PolygonLayer.from_geopandas(
    counties_gdf,
    get_fill_color=[200, 200, 200, 100],
    get_line_color=[50, 50, 50, 255],
    get_line_width=1,
    pickable=True
)

# Combine layers in a map
m = Map(layers=[polygon_layer, point_layer])
m.to_html("multi_layer.html")

---

Publication Maps with Cartopy

cartopy provides matplotlib-based map projections with proper geographic axes — essential for publication-quality maps requiring explicit projection control, graticules, and natural features.

Installation

pip install cartopy

Basic Map with Projection

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

fig, ax = plt.subplots(
    figsize=(12, 8),
    subplot_kw={"projection": ccrs.AlbersEqualArea(
        central_longitude=-96, central_latitude=37.5,
        standard_parallels=(29.5, 45.5)
    )}
)

# Add natural features
ax.add_feature(cfeature.STATES, linewidth=0.5, edgecolor="gray")
ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle="--")

# Plot geopandas data on the cartopy axes
# transform= tells cartopy what CRS the data is in
gdf.plot(ax=ax, column="poverty_rate", cmap="YlOrRd", legend=True,
         transform=ccrs.PlateCarree())  # PlateCarree = lon/lat (EPSG:4326)

ax.set_extent([-125, -66, 24, 50], crs=ccrs.PlateCarree())  # Continental US
ax.set_title("Poverty Rate by County", fontsize=14)
plt.tight_layout()
plt.savefig("cartopy_map.png", dpi=300, bbox_inches="tight")

Common Cartopy Projections

Key cartopy + geopandas Integration

The transform= parameter is critical: it declares what CRS the input data is in, so cartopy can reproject it to the axes projection. If your GeoDataFrame is in EPSG:4326, use transform=ccrs.PlateCarree(). If already projected, match accordingly.

---

Datashader for Massive Point Datasets

datashader rasterizes point or line data into pixel grids before rendering, enabling visualization of billions of points without browser or memory limitations.

Installation

pip install datashader

Basic Point Density

import datashader as ds
import datashader.transfer_functions as tf
import pandas as pd

# Extract coordinates (datashader works on plain DataFrames, not GeoDataFrames)
df = pd.DataFrame({"x": gdf.geometry.x, "y": gdf.geometry.y})

canvas = ds.Canvas(plot_width=800, plot_height=600)
agg = canvas.points(df, "x", "y")
img = tf.shade(agg, cmap="viridis")
img = tf.set_background(img, "black")

# Save as PNG
from datashader.utils import export_image
export_image(img, "point_density", background="black")

When to Use Each Tool

---

Multi-Panel Maps

Matplotlib Subplots

fig, axes = plt.subplots(1, 3, figsize=(18, 6))

for ax, (col, title) in zip(axes, [
    ("poverty_rate", "Poverty Rate"),
    ("median_income", "Median Income"),
    ("graduation_rate", "Graduation Rate")
]):
    gdf.plot(column=col, ax=ax, cmap="YlOrRd", legend=True, scheme="FisherJenks", k=5)
    ax.set_title(title)
    ax.set_axis_off()

plt.tight_layout()
plt.savefig("panel_maps.png", dpi=300, bbox_inches="tight")

---

Exporting Maps

Static Formats

# PNG (default for web and reports)
plt.savefig("map.png", dpi=300, bbox_inches="tight")

# SVG (scalable, for publications)
plt.savefig("map.svg", bbox_inches="tight")

# PDF (for print)
plt.savefig("map.pdf", bbox_inches="tight")

Interactive Formats

# Folium to HTML
m = gdf.explore(column="value")
m.save("map.html")

# Lonboard to HTML
from lonboard import viz
m = viz(gdf)
m.to_html("map.html")

---

References and Further Reading

Jordahl, K. et al. (2024). geopandas: Making maps. https://geopandas.org/en/stable/docs/user_guide/mapping.html

Rey, S.J., Arribas-Bel, D., and Wolf, L.J. (2023). Geographic Data Science with Python, Ch. 5: "Choropleth mapping." https://geographicdata.science/book/

Dorman, M., Graser, A., Nowosad, J., and Lovelace, R. (2025). Geocomputation with Python, Ch. 8: "Making maps with Python." https://py.geocompx.org/

contextily documentation. https://contextily.readthedocs.io/

folium documentation. https://python-visualization.github.io/folium/

lonboard documentation. https://developmentseed.org/lonboard/

mapclassify documentation. https://pysal.org/mapclassify/