---
name: jaechang-hits/geopandas-geospatial
source: https://app.decimal.ai/s/jaechang-hits-geopandas-geospatial@1/SKILL.md
source_sha256: 13726c434f1f
---

# GeoPandas Geospatial Analysis

## Overview

GeoPandas extends pandas with spatial operations on geometric types, combining pandas DataFrames with Shapely geometries and Fiona for file I/O. It enables reading, writing, manipulating, and visualizing geospatial vector data (points, lines, polygons) with a familiar pandas-like API.

## When to Use

- Reading and writing spatial file formats (Shapefile, GeoJSON, GeoPackage, Parquet)
- Performing spatial joins between geographic datasets (points in polygons, nearest neighbors)
- Running overlay operations (intersection, union, difference, clipping)
- Computing geometric properties (area, distance, buffer, centroid)
- Creating choropleth maps and interactive web maps
- Reprojecting data between coordinate reference systems
- Aggregating spatial features by attribute (dissolve)
- For raster data analysis, use rasterio/xarray instead
- For large-scale distributed geospatial, consider Dask-GeoPandas or Apache Sedona

## Prerequisites

```bash
pip install geopandas matplotlib
# Optional:
# pip install folium       — interactive maps
# pip install mapclassify  — classification schemes for choropleth
# pip install contextily   — basemaps
# pip install pyarrow      — faster I/O (2-4x speedup)
# pip install psycopg2 geoalchemy2  — PostGIS support
```

## Quick Start

```python
import geopandas as gpd

# Read spatial data
gdf = gpd.read_file("data.geojson")
print(f"Shape: {gdf.shape}, CRS: {gdf.crs}")
print(f"Geometry types: {gdf.geometry.geom_type.unique()}")

# Reproject, compute area, save
gdf_proj = gdf.to_crs("EPSG:3857")
gdf_proj['area_m2'] = gdf_proj.geometry.area
gdf_proj.to_file("output.gpkg")

# Quick map
gdf.plot(column='population', legend=True, figsize=(10, 8))
```

## Core API

### 1. Data I/O

```python
import geopandas as gpd

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

# Filtered reading (load only needed data)
gdf = gpd.read_file("data.gpkg", bbox=(xmin, ymin, xmax, ymax))
gdf = gpd.read_file("data.gpkg", columns=["name", "geometry"])
gdf = gpd.read_file("data.gpkg", where="population > 10000")

# Arrow acceleration (2-4x faster)
gdf = gpd.read_file("data.gpkg", use_arrow=True)

# Parquet/Feather (columnar, fast, preserves CRS)
gdf = gpd.read_parquet("data.parquet")
gdf.to_parquet("output.parquet")

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

# Write
gdf.to_file("output.gpkg")             # GeoPackage (recommended)
gdf.to_file("output.shp")              # Shapefile
gdf.to_file("output.geojson", driver="GeoJSON")
```

### 2. CRS Management

```python
# Check current CRS
print(gdf.crs)                  # e.g., EPSG:4326
print(gdf.crs.is_geographic)    # True for lat/lon
print(gdf.crs.is_projected)     # True for meters

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

# Set CRS (only when metadata missing, does NOT transform coordinates)
gdf = gdf.set_crs("EPSG:4326")

# Estimate appropriate UTM zone
utm_crs = gdf.estimate_utm_crs()
gdf_utm = gdf.to_crs(utm_crs)
```

**Common EPSG codes**:

| Code | Name | Use |
|------|------|-----|
| 4326 | WGS 84 | GPS coordinates, web data |
| 3857 | Web Mercator | Web mapping (Google/OSM tiles) |
| 326xx | UTM zones (N) | Area/distance calculations |
| 5070 | Albers Equal Area (US) | Area-preserving US maps |

### 3. Geometric Operations

```python
# Buffer (expand/erode geometry by distance)
buffered = gdf.geometry.buffer(100)      # 100 units (meters if projected)
eroded = gdf.geometry.buffer(-50)        # Negative = erosion

# Simplify (reduce complexity)
simplified = gdf.geometry.simplify(tolerance=10, preserve_topology=True)

# Centroid, convex hull, envelope
centroids = gdf.geometry.centroid
hulls = gdf.geometry.convex_hull
bounds = gdf.geometry.envelope

# Union all geometries
unified = gdf.geometry.union_all()

# Affine transformations
rotated = gdf.geometry.rotate(angle=45, origin='center')
scaled = gdf.geometry.scale(xfact=2.0, yfact=2.0)
translated = gdf.geometry.translate(xoff=100, yoff=50)

# Geometric properties
areas = gdf.geometry.area           # Use projected CRS for accuracy
lengths = gdf.geometry.length       # Perimeter for polygons
is_valid = gdf.geometry.is_valid    # Validate geometry
total = gdf.geometry.total_bounds   # [minx, miny, maxx, maxy]
```

### 4. Spatial Analysis

```python
# Spatial join (combine datasets by spatial relationship)
joined = gpd.sjoin(points_gdf, polygons_gdf, predicate='intersects')
joined = gpd.sjoin(gdf1, gdf2, predicate='within')
joined = gpd.sjoin(gdf1, gdf2, predicate='contains', how='left')

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

# Overlay operations (set-theoretic)
intersection = gpd.overlay(gdf1, gdf2, how='intersection')
union = gpd.overlay(gdf1, gdf2, how='union')
difference = gpd.overlay(gdf1, gdf2, how='difference')
sym_diff = gpd.overlay(gdf1, gdf2, how='symmetric_difference')

# Dissolve (aggregate by attribute)
dissolved = gdf.dissolve(by='region', aggfunc='sum')
dissolved = gdf.dissolve(by='region', aggfunc={'population': 'sum', 'area': 'mean'})

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

# Distance calculations (use projected CRS)
distances = gdf.geometry.distance(single_point)

# Spatial predicates
within_mask = gdf1.geometry.within(gdf2.geometry)
intersects_mask = gdf1.geometry.intersects(gdf2.geometry)
```

### 5. Visualization

```python
import matplotlib.pyplot as plt

# Basic plot
gdf.plot(figsize=(10, 8))

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

# Classification schemes (requires mapclassify)
gdf.plot(column='population', scheme='quantiles', k=5, legend=True)
gdf.plot(column='population', scheme='fisher_jenks', k=5, legend=True)

# Multi-layer map
fig, ax = plt.subplots(figsize=(12, 10))
polygons_gdf.plot(ax=ax, color='lightblue', edgecolor='black')
points_gdf.plot(ax=ax, color='red', markersize=10)
roads_gdf.plot(ax=ax, color='gray', linewidth=0.5)
ax.set_title('Multi-layer Map')
ax.set_axis_off()

# Interactive map (requires folium)
m = gdf.explore(column='population', cmap='YlOrRd', legend=True,
                tooltip=['name', 'population'])
m.save('map.html')

# Multi-layer interactive
m = gdf1.explore(color='blue', name='Layer 1')
gdf2.explore(m=m, color='red', name='Layer 2')
import folium
folium.LayerControl().add_to(m)

# Basemap (requires contextily)
import contextily as ctx
gdf_wm = gdf.to_crs(epsg=3857)
ax = gdf_wm.plot(alpha=0.5, figsize=(10, 10))
ctx.add_basemap(ax)
```

## Key Concepts

### Data Structures

- **GeoSeries**: Pandas Series of Shapely geometries with spatial methods (area, distance, buffer, etc.)
- **GeoDataFrame**: Pandas DataFrame with one or more geometry columns. One column is the "active geometry" used by spatial methods

```python
from shapely.geometry import Point

# Create from coordinates
gdf = gpd.GeoDataFrame(
    {'name': ['A', 'B'], 'value': [10, 20]},
    geometry=[Point(0, 0), Point(1, 1)],
    crs="EPSG:4326"
)

# Multiple geometry columns
gdf['centroid'] = gdf.geometry.centroid
gdf = gdf.set_geometry('centroid')  # Switch active geometry
```

### CRS Rules for Spatial Operations

- **Always check CRS** before any spatial operation: `print(gdf.crs)`
- **Match CRS** before spatial joins, overlays, or distance calculations
- **Use projected CRS** (meters) for area/distance — geographic CRS (degrees) gives wrong results
- **set_crs()** only adds metadata; **to_crs()** transforms coordinates

### Spatial Indexing

GeoPandas automatically creates spatial indexes (R-tree) for sjoin, overlay, and other spatial operations. For manual queries:

```python
sindex = gdf.sindex
possible_idx = list(sindex.intersection((xmin, ymin, xmax, ymax)))
```

## Common Workflows

### Workflow 1: Load → Transform → Analyze → Export

```python
import geopandas as gpd

# Load
gdf = gpd.read_file("parcels.shp")
print(f"CRS: {gdf.crs}, Rows: {len(gdf)}")

# Transform to projected CRS for measurements
gdf = gdf.to_crs(gdf.estimate_utm_crs())

# Analyze
gdf['area_ha'] = gdf.geometry.area / 10000  # hectares
gdf['perimeter_m'] = gdf.geometry.length
print(f"Total area: {gdf['area_ha'].sum():.1f} ha")

# Export
gdf.to_file("parcels_analyzed.gpkg")
```

### Workflow 2: Spatial Join and Aggregate

```python
# Count points per polygon
points_in_poly = gpd.sjoin(points_gdf, polygons_gdf, predicate='within')
counts = points_in_poly.groupby('index_right').agg(
    point_count=('geometry', 'size'),
    total_value=('value', 'sum')
)
result = polygons_gdf.merge(counts, left_index=True, right_index=True, how='left')
result['point_count'] = result['point_count'].fillna(0)
print(f"Polygons with points: {(result['point_count'] > 0).sum()}/{len(result)}")
```

### Workflow 3: Multi-Source Integration

```python
# Read from different sources, ensure matching CRS
roads = gpd.read_file("roads.shp")
buildings = gpd.read_file("buildings.geojson")
parcels = gpd.read_postgis("SELECT * FROM parcels", con=engine, geom_col='geom')

target_crs = roads.crs
buildings = buildings.to_crs(target_crs)
parcels = parcels.to_crs(target_crs)

# Find buildings within 50m of roads
buildings_near_roads = gpd.sjoin_nearest(
    buildings, roads, max_distance=50, distance_col='road_dist'
)
print(f"Buildings near roads: {len(buildings_near_roads)}/{len(buildings)}")
```

## Key Parameters

| Parameter | Function | Default | Effect |
|-----------|----------|---------|--------|
| `predicate` | sjoin | `'intersects'` | Spatial relationship: intersects, within, contains, touches, crosses |
| `how` | sjoin, overlay | `'inner'` | Join type: inner, left, right |
| `max_distance` | sjoin_nearest | None | Search radius limit (improves performance) |
| `k` | sjoin_nearest | 1 | Number of nearest neighbors to find |
| `tolerance` | simplify | Required | Douglas-Peucker tolerance (in CRS units) |
| `preserve_topology` | simplify | True | Prevents self-intersections |
| `resolution` | buffer | 16 | Number of segments for buffer curves |
| `aggfunc` | dissolve | 'first' | Aggregation function for non-geometry columns |
| `use_arrow` | read_file | False | Enable Arrow acceleration (2-4x faster) |
| `scheme` | plot | None | Classification: quantiles, equal_interval, fisher_jenks |
| `k` | plot (with scheme) | 5 | Number of classification bins |

## Best Practices

1. **Always check CRS** before spatial operations — mismatched CRS gives wrong or empty results
2. **Use projected CRS** for area and distance calculations — geographic CRS (degrees) is meaningless for measurements
3. **Match CRS before joins** — `gdf2 = gdf2.to_crs(gdf1.crs)` before `gpd.sjoin(gdf1, gdf2)`
4. **Validate geometries** with `.is_valid` before complex operations — invalid geometries cause silent errors
5. **Use GeoPackage** over Shapefile — no 10-char column name limit, supports multiple layers, better performance
6. **Filter during read** — use `bbox`, `columns`, `where` to load only needed data for large files
7. **Set max_distance in sjoin_nearest** — unbounded nearest-neighbor search is slow on large datasets
8. **Use `.copy()` when modifying geometry** — avoid unintended side effects on original GeoDataFrame

## Common Recipes

### Recipe: Count Points in Polygons

When to use: Aggregate point observations (sample sites, observations) by region (county, watershed).

```python
import geopandas as gpd

regions = gpd.read_file("regions.geojson")
points = gpd.read_file("observations.geojson").to_crs(regions.crs)

joined = gpd.sjoin(points, regions, how="inner", predicate="within")
counts = joined.groupby("index_right").size().rename("point_count")
regions_with_counts = regions.join(counts).fillna(0)
regions_with_counts["point_count"] = regions_with_counts["point_count"].astype(int)
print(regions_with_counts[["name", "point_count"]].sort_values("point_count", ascending=False).head())
```

### Recipe: Buffer Around Points and Dissolve Overlaps

When to use: Create service areas or catchment zones from point locations.

```python
import geopandas as gpd

facilities = gpd.read_file("facilities.geojson")
utm_crs = facilities.estimate_utm_crs()   # Project to meters
facilities_m = facilities.to_crs(utm_crs)

# 1 km buffer
buffers = facilities_m.copy()
buffers["geometry"] = facilities_m.geometry.buffer(1000)

# Dissolve overlapping buffers into single polygon
service_area = buffers.dissolve()
service_area_wgs84 = service_area.to_crs("EPSG:4326")
service_area_wgs84.to_file("service_area.geojson", driver="GeoJSON")
print(f"Service area: {service_area_wgs84.geometry.area.sum():.0f} sq degrees")
```

## Troubleshooting

| Problem | Cause | Solution |
|---------|-------|----------|
| Empty spatial join result | CRS mismatch between GeoDataFrames | Ensure matching CRS: `gdf2 = gdf2.to_crs(gdf1.crs)` |
| Wrong area/distance values | Using geographic CRS (degrees) | Reproject to projected CRS: `gdf.to_crs(gdf.estimate_utm_crs())` |
| `DriverError` on read | Missing driver or corrupt file | Check file exists; try `driver="GeoJSON"` explicitly |
| Geometry column lost after merge | Called `df.merge(gdf)` instead of `gdf.merge(df)` | Always call merge ON the GeoDataFrame |
| `TopologicalError` on overlay | Invalid geometries | Fix with `gdf.geometry = gdf.geometry.buffer(0)` |
| Slow sjoin_nearest | No distance limit on large dataset | Set `max_distance` parameter |
| Folium map blank | Geometries not in EPSG:4326 | Reproject to WGS84: `gdf.to_crs(epsg=4326).explore()` |
| Shapefile column names truncated | Shapefile 10-char limit | Use GeoPackage instead: `gdf.to_file("out.gpkg")` |

## References

- GeoPandas documentation: https://geopandas.org/
- GeoPandas GitHub: https://github.com/geopandas/geopandas
- Shapely documentation: https://shapely.readthedocs.io/
- EPSG registry: https://epsg.io/

## Related Skills

- **matplotlib-scientific-plotting** — advanced map styling and figure export
- **polars-dataframes** — high-performance tabular analysis before/after spatial operations
- **folium** — advanced interactive web mapping beyond geopandas.explore()