This skill covers the complete vector GIS workflow in Python using GeoPandas: reading and
writing spatial files, coordinate reference system (CRS) management, spatial joins, overlay
operations, geometry processing, spatial statistics, and choropleth mapping. It also covers
the integration point between vector data (GeoPandas) and raster data (rasterio/rasterstats).
GeoPandas extends pandas DataFrames with a geometry column powered by Shapely, making it
straightforward to apply familiar tabular operations alongside spatial predicates.
Setup
# Recommended: install via conda for reliable binary dependencies
conda create -n gis python=3.11
conda activate gis
conda install -c conda-forge geopandas rasterio rasterstats matplotlib requests fiona
# Or with pip (may require system GDAL/GEOS on Linux/macOS):
pip install geopandas rasterio rasterstats matplotlib requests fiona
"""
Read a vector file (shapefile, GeoJSON, GPKG, etc.) into a GeoDataFrame.
Parameters
----------
path : str or Path
Path to the file. Supports any format recognised by Fiona/GDAL
(`.shp`, `.geojson`, `.gpkg`, `.gdb`, etc.).
crs_target : str, optional
EPSG code or PROJ string to reproject to, e.g. 'EPSG:3857'.
If None, keeps the source CRS.
bbox : tuple, optional
Spatial filter (minx, miny, maxx, maxy) in the source CRS.
Returns
-------
gpd.GeoDataFrame
Examples
--------
>>> countries = load_shapefile("ne_10m_admin_0_countries.shp", crs_target="EPSG:4326")
>>> europe = load_shapefile("ne_10m_admin_0_countries.shp",
... crs_target="EPSG:3035",
... bbox=(-25, 34, 45, 72))
"""
"""
Join each feature in gdf1 to its nearest neighbour in gdf2.
Both GeoDataFrames must share the same (projected) CRS. Distances are
reported in the units of that CRS (metres for metric projections).
Parameters
----------
gdf1 : gpd.GeoDataFrame
Left GeoDataFrame (e.g. census tracts to enrich).
gdf2 : gpd.GeoDataFrame
Right GeoDataFrame to join from (e.g. river segments).
max_distance : float, optional
If provided, only joins where the nearest feature is within this
distance. Features beyond max_distance get NaN attribute values.
Returns
-------
gpd.GeoDataFrame
gdf1 augmented with columns from gdf2 and a 'distance' column.
"""
"""
Tag each point (or polygon centroid) with the polygon it falls within.
Parameters
----------
gdf_points : gpd.GeoDataFrame
Point features to tag.
gdf_polygons : gpd.GeoDataFrame
Polygon features containing spatial attributes to transfer.
how : str
pandas-style join type: 'left', 'inner', 'right'.
Returns
-------
gpd.GeoDataFrame
"""
"""
Clip a GeoDataFrame to an axis-aligned bounding box.
Parameters
----------
gdf : gpd.GeoDataFrame
bbox : tuple
(minx, miny, maxx, maxy) in the same CRS as gdf.
Returns
-------
gpd.GeoDataFrame with geometries clipped (not just filtered).
"""
"""
Buffer features and optionally dissolve overlapping buffers.
Parameters
----------
gdf : gpd.GeoDataFrame
Must be in a projected (metric) CRS for meaningful distance buffers.
distance : float
Buffer distance in CRS units (metres for UTM/EPSG:3857).
dissolve_by : str, optional
Column to group/dissolve by. Pass None to dissolve all into one.
cap_style : int
Shapely cap style: 1=round, 2=flat, 3=square.
Returns
-------
gpd.GeoDataFrame
"""
"geometry"
if
is
not
None
False
else
True
return
def
overlay_intersection
gdf1: gpd.GeoDataFrame,
gdf2: gpd.GeoDataFrame,
"""
Return the geometric intersection of two GeoDataFrames (polygon-on-polygon).
Attribute columns from both inputs are retained.
"""
"""
Compute summary statistics (area, perimeter, centroid, value statistics)
for each feature or group of features.
Parameters
----------
gdf : gpd.GeoDataFrame
Must be in a projected CRS to get meaningful area/length values.
value_col : str
Numeric attribute column to summarise.
group_col : str, optional
Column to group by before aggregating.
Returns
-------
pd.DataFrame with columns: count, area_km2, length_km, mean, std, min, max.
"""
"""
Plot a choropleth map of a numeric column.
Parameters
----------
gdf : gpd.GeoDataFrame
col : str
Column to map.
title : str
Figure title.
cmap : str
Matplotlib colormap.
figsize : tuple
legend : bool
Whether to include a colorbar legend.
edgecolor : str
Polygon border color.
linewidth : float
Polygon border width.
output_path : str, optional
If given, save the figure here.
Returns
-------
matplotlib Figure
"""
1
1
"label"
"orientation"
"horizontal"
"pad"
0.04
"color"
"lightgrey"
"label"
"No data"
14
10
if
150
"tight"
print
f"Saved: {output_path}"
return
Example 1: Land-Use Area Statistics for Administrative Regions
Download Natural Earth administrative boundaries and a land-use GeoJSON, then compute
area statistics for each land-use class within each admin region.
"""
example_landuse_stats.py
------------------------
Calculate land-use area statistics per administrative region using
spatial overlay and geopandas aggregation.
"""import requests
import zipfile
import io
from pathlib import Path
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
from geopandas_gis import (
load_shapefile,
overlay_intersection,
compute_spatial_stats,
plot_choropleth,
)
DATA_DIR = Path("gis_data")
DATA_DIR.mkdir(exist_ok=True)
defdownload_natural_earth(scale: str = "10m", category: str = "cultural",
name: str = "admin_1_states_provinces") -> Path:
"""Download a Natural Earth shapefile and return the .shp path."""
url = (f"https://naciscdn.org/naturalearth/{scale}/{category}/"f"ne_{scale}_{name}.zip")
out_zip = DATA_DIR / f"ne_{scale}_{name}.zip"ifnot out_zip.exists():
print(f"Downloading {url} ...")
r = requests.get(url, timeout=60)
r.raise_for_status()
out_zip.write_bytes(r.content)
# Extract
out_dir = DATA_DIR / f"ne_{scale}_{name}"ifnot out_dir.exists():
with zipfile.ZipFile(out_zip) as zf:
zf.extractall(out_dir)
shp = next(out_dir.glob("*.shp"))
return shp
# ---- 1. Load administrative boundaries (example: Germany federal states) ----
states_shp = download_natural_earth("10m", "cultural", "admin_1_states_provinces")
states = load_shapefile(states_shp, crs_target="EPSG:3035")
germany = states[states["admin"] == "Germany"].copy()
print(f"Loaded {len(germany)} German federal states")
# ---- 2. Create synthetic land-use data for demonstration ----# In a real workflow, replace this with an actual land-use shapefile# (e.g. CORINE Land Cover from https://land.copernicus.eu)import numpy as np
from shapely.geometry import Polygon
rng = np.random.default_rng(42)
land_use_types = ["Forest", "Agriculture", "Urban", "Water", "Shrubland"]
# Generate random polygons inside Germany's bounding box
minx, miny, maxx, maxy = germany.total_bounds
synthetic_polygons = []
for _ inrange(500):
cx = rng.uniform(minx, maxx)
cy = rng.uniform(miny, maxy)
size = rng.uniform(5000, 50000) # 5–50 km²
poly = Polygon([
(cx - size, cy - size), (cx + size, cy - size),
(cx + size, cy + size), (cx - size, cy + size),
])
lu_type = rng.choice(land_use_types)
synthetic_polygons.append({"land_use": lu_type, "geometry": poly})
land_use = gpd.GeoDataFrame(synthetic_polygons, crs="EPSG:3035")
land_use = gpd.clip(land_use, germany.union_all())
# ---- 3. Overlay: intersect land-use with state boundaries ----
lu_by_state = overlay_intersection(land_use, germany[["NAME_1", "geometry"]])
lu_by_state["area_km2"] = lu_by_state.geometry.area / 1e6# ---- 4. Summarise: total area by state and land-use type ----
summary = (
lu_by_state.groupby(["NAME_1", "land_use"])["area_km2"]
.sum()
.unstack(fill_value=0)
.round(1)
)
print("\nLand-use area (km²) by federal state:")
print(summary.to_string())
# ---- 5. Compute forest fraction per state ----
summary["forest_fraction"] = (
summary.get("Forest", 0) / summary.sum(axis=1) * 100
)
germany = germany.merge(
summary[["forest_fraction"]].reset_index(),
left_on="NAME_1", right_on="NAME_1", how="left",
)
# ---- 6. Choropleth of forest fraction ----
fig = plot_choropleth(
germany.to_crs("EPSG:4326"),
col="forest_fraction",
title="Estimated Forest Coverage by German Federal State (%)",
cmap="Greens",
output_path="germany_forest_fraction.png",
)
plt.show()
# ---- 7. Bar chart: top 5 states by urban area ----if"Urban"in summary.columns:
top_urban = summary["Urban"].sort_values(ascending=False).head(5)
fig2, ax2 = plt.subplots(figsize=(8, 4))
top_urban.plot.barh(ax=ax2, color="steelblue")
ax2.set_xlabel("Urban area (km²)")
ax2.set_title("Top 5 German States by Urban Area")
fig2.tight_layout()
fig2.savefig("germany_urban_top5.png", dpi=150)
print("Saved: germany_urban_top5.png")
Example 2: Census Tracts Within 5 km of a River Network
Find all census tracts whose centroid falls within a 5 km buffer of a river network,
then compute population statistics for the flood-risk zone.
"""
example_river_buffer.py
-----------------------
Identify census tracts within 5 km of rivers and compute population at risk.
Uses synthetic data; swap in real shapefiles for production use.
"""import numpy as np
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
from shapely.geometry import LineString, Point
from geopandas_gis import (
buffer_and_dissolve,
spatial_join_within,
compute_spatial_stats,
plot_choropleth,
)
rng = np.random.default_rng(0)
# ---- 1. Create synthetic river network (LineStrings) ----# Replace with: gpd.read_file("rivers.shp").to_crs("EPSG:32632")defmake_river(start, n_segments=20, spread=5000):
coords = [start]
for _ inrange(n_segments):
last = coords[-1]
coords.append((last[0] + rng.uniform(1000, 8000),
last[1] + rng.uniform(-spread, spread)))
return LineString(coords)
rivers_data = [
{"river_name": "Main River", "geometry": make_river((400000, 5540000))},
{"river_name": "Rhine", "geometry": make_river((340000, 5600000), spread=8000)},
{"river_name": "Tributary A", "geometry": make_river((380000, 5560000), n_segments=10)},
]
rivers = gpd.GeoDataFrame(rivers_data, crs="EPSG:32632")
print(f"River network: {len(rivers)} segments, "f"total length {rivers.geometry.length.sum()/1000:.1f} km")
# ---- 2. Create 5 km buffer around all rivers, dissolved ----
river_buffer = buffer_and_dissolve(rivers, distance=5000)
river_buffer["zone"] = "5km_buffer"print(f"Buffer area: {river_buffer.geometry.area.sum()/1e6:.1f} km²")
# ---- 3. Create synthetic census tracts with population ----
minx, miny, maxx, maxy = rivers.total_bounds
n_tracts = 300
tract_data = []
for i inrange(n_tracts):
cx = rng.uniform(minx - 10000, maxx + 10000)
cy = rng.uniform(miny - 10000, maxy + 10000)
w = rng.uniform(1500, 4000)
from shapely.geometry import box as sbox
poly = sbox(cx - w, cy - w, cx + w, cy + w)
pop = int(rng.lognormal(8, 0.8)) # log-normal population
tract_data.append({"tract_id": i, "population": pop, "geometry": poly})
tracts = gpd.GeoDataFrame(tract_data, crs="EPSG:32632")
tracts["centroid_geom"] = tracts.geometry.centroid
# ---- 4. Find tracts whose centroid is inside the buffer ----
tract_centroids = tracts.copy()
tract_centroids["geometry"] = tracts["centroid_geom"]
in_buffer = spatial_join_within(tract_centroids, river_buffer[["zone", "geometry"]])
at_risk_ids = set(in_buffer["tract_id"].unique())
tracts["at_risk"] = tracts["tract_id"].isin(at_risk_ids)
print(f"\nTracts at risk (within 5 km): {tracts['at_risk'].sum()} / {len(tracts)}")
print(f"Population at risk: {tracts.loc[tracts['at_risk'], 'population'].sum():,}")
print(f"Total population: {tracts['population'].sum():,}")
pct = tracts.loc[tracts['at_risk'], 'population'].sum() / tracts['population'].sum() * 100print(f"Fraction at risk: {pct:.1f}%")
# ---- 5. Map ----
fig, ax = plt.subplots(figsize=(12, 10))
# Background tracts
tracts.plot(ax=ax, column="at_risk",
cmap="RdYlGn_r", alpha=0.6,
edgecolor="grey", linewidth=0.3)
# River buffer overlay
river_buffer.plot(ax=ax, color="steelblue", alpha=0.25, label="5 km buffer")
# River lines
rivers.plot(ax=ax, color="blue", linewidth=1.5, label="Rivers")
# Legend patchesfrom matplotlib.patches import Patch
legend_elements = [
Patch(facecolor="firebrick", alpha=0.7, label="At-risk tracts"),
Patch(facecolor="green", alpha=0.7, label="Safe tracts"),
Patch(facecolor="steelblue", alpha=0.4, label="5 km river buffer"),
]
ax.legend(handles=legend_elements, loc="upper left")
ax.set_title("Census Tracts Within 5 km of River Network\n"f"Population at risk: {pct:.1f}%", fontsize=13)
ax.set_axis_off()
fig.tight_layout()
fig.savefig("river_flood_risk.png", dpi=150, bbox_inches="tight")
print("\nSaved: river_flood_risk.png")
plt.show()
# ---- 6. Spatial statistics for at-risk vs safe tracts ----
stats_at_risk = compute_spatial_stats(
tracts[tracts["at_risk"]].to_crs("EPSG:32632"),
value_col="population",
)
stats_safe = compute_spatial_stats(
tracts[~tracts["at_risk"]].to_crs("EPSG:32632"),
value_col="population",
)
print("\nSpatial statistics — at-risk tracts:")
print(stats_at_risk.to_string(index=False))
print("\nSpatial statistics — safe tracts:")
print(stats_safe.to_string(index=False))
Raster-Vector Integration with rasterio
"""
raster_vector.py
----------------
Extract raster statistics (mean elevation, median NDVI, etc.) for each
polygon in a GeoDataFrame using rasterstats.
"""import numpy as np
import geopandas as gpd
import rasterio
from rasterio.mask import mask as rio_mask
from rasterstats import zonal_stats
from pathlib import Path
defzonal_statistics(
gdf: gpd.GeoDataFrame,
raster_path: str,
stats: list = None,
all_touched: bool = False,
) -> gpd.GeoDataFrame:
"""
Compute zonal statistics from a raster for each polygon.
Parameters
----------
gdf : gpd.GeoDataFrame
Polygon GeoDataFrame. Will be reprojected to match the raster CRS.
raster_path : str
Path to a GeoTIFF or other GDAL-readable raster.
stats : list, optional
Statistics to compute. Defaults to ['mean', 'std', 'min', 'max', 'count'].
all_touched : bool
If True, include all cells touching a polygon (not just cell centres inside).
Returns
-------
gpd.GeoDataFrame
Input GeoDataFrame with new columns for each requested statistic.
"""if stats isNone:
stats = ["mean", "std", "min", "max", "count"]
with rasterio.open(raster_path) as src:
raster_crs = src.crs.to_epsg()
gdf_proj = gdf.to_crs(f"EPSG:{raster_crs}")
results = zonal_stats(
gdf_proj,
raster_path,
stats=stats,
all_touched=all_touched,
geojson_out=False,
)
stats_df = gpd.GeoDataFrame(results, index=gdf.index)
return gdf.join(stats_df)
defclip_raster_to_polygon(
raster_path: str,
gdf: gpd.GeoDataFrame,
output_path: str,
) -> None:
"""
Clip a raster to the union of all polygons in a GeoDataFrame.
Parameters
----------
raster_path : str
gdf : gpd.GeoDataFrame
output_path : str
Path for the output clipped GeoTIFF.
"""with rasterio.open(raster_path) as src:
gdf_proj = gdf.to_crs(src.crs)
shapes = [geom.__geo_interface__ for geom in gdf_proj.geometry]
out_image, out_transform = rio_mask(src, shapes, crop=True)
out_meta = src.meta.copy()
out_meta.update({
"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform,
})
with rasterio.open(output_path, "w", **out_meta) as dest:
dest.write(out_image)
print(f"Clipped raster saved: {output_path}")
CRS Quick Reference
EPSG Code
Description
Use case
4326
WGS 84 geographic (lat/lon degrees)
Default GPS, web display
3857
Web Mercator (metres)
Tile-based web maps
3035
ETRS89 / LAEA Europe (metres)
Area calculations in Europe
32632
UTM zone 32N (metres)
Local metric in Central Europe
4269
NAD83 geographic
North America
5070
Albers Equal Area CONUS (metres)
Area statistics for US lower 48
Always use a projected (metric) CRS before computing area, length, or buffer distance.
Tips and Best Practices
Check CRS early: call gdf.crs immediately after loading. A None CRS means no
projection information was embedded; you must set it manually with gdf.set_crs(...).
sjoin vs sjoin_nearest: use sjoin (with predicates within, intersects,
contains) when an exact spatial relationship is required; use sjoin_nearest when
you need the closest feature regardless of whether it intersects.
Large files: use bbox parameter in gpd.read_file to read only the spatial subset
you need, avoiding loading millions of features into memory.
GeoPandas 1.0+: unary_union is deprecated in favour of union_all().
Shapely 2.0: geometry operations are significantly faster. Ensure shapely>=2.0.
Avoid apply with geometry: prefer vectorised methods (gdf.geometry.buffer,
gdf.geometry.area) over row-wise apply(lambda row: ...) for performance.
