Geopandas Spatial

CRS Transformations and Projections

# Check and set CRS
print(gdf.crs)                                # e.g., EPSG:4326 (WGS84)
gdf = gdf.set_crs("EPSG:4326")               # assign if missing

# Reproject to a different CRS
gdf_proj = gdf.to_crs("EPSG:3857")           # Web Mercator
gdf_utm = gdf.to_crs("EPSG:32633")           # UTM Zone 33N

# Area calculation (reproject to equal-area CRS first)
gdf_equal = gdf.to_crs("ESRI:54009")         # Mollweide equal-area
gdf_equal["area_km2"] = gdf_equal.geometry.area / 1e6

Spatial Operations

from shapely.geometry import Point, Polygon, box

# Create geometries
point = Point(-73.985, 40.748)
polygon = Polygon([(-74, 40.7), (-74, 40.8), (-73.9, 40.8), (-73.9, 40.7)])
bbox = box(-74.05, 40.68, -73.90, 40.82)

# Spatial joins
joined = gpd.sjoin(points_gdf, polygons_gdf, how="inner", predicate="within")

# Buffer around geometries (in CRS units)
gdf_buffered = gdf.copy()
gdf_buffered["geometry"] = gdf.geometry.buffer(1000)  # 1000m if projected CRS

# Dissolve by attribute (merge geometries)
dissolved = gdf.dissolve(by="region", aggfunc="sum")

# Overlay operations
intersection = gpd.overlay(gdf1, gdf2, how="intersection")
union = gpd.overlay(gdf1, gdf2, how="union")
difference = gpd.overlay(gdf1, gdf2, how="difference")

# Clip to bounding box or polygon
clipped = gpd.clip(gdf, mask=bbox_gdf)

# Nearest join
nearest = gpd.sjoin_nearest(points_gdf, target_gdf, how="left", distance_col="dist_m")

Map Visualization

import matplotlib.pyplot as plt

# Basic plot
ax = gdf.plot(figsize=(12, 8), edgecolor="black", linewidth=0.5)
ax.set_title("Geographic Boundaries")

# Choropleth map
fig, ax = plt.subplots(1, 1, figsize=(14, 10))
gdf.plot(column="population", cmap="YlOrRd", legend=True,
         legend_kwds={"label": "Population"}, ax=ax,
         edgecolor="gray", linewidth=0.3)
ax.set_axis_off()
plt.savefig("choropleth.pdf", bbox_inches="tight")

# Multi-layer map
fig, ax = plt.subplots(figsize=(12, 8))
polygons_gdf.plot(ax=ax, color="lightblue", edgecolor="gray")
points_gdf.plot(ax=ax, color="red", markersize=5)
lines_gdf.plot(ax=ax, color="darkblue", linewidth=1)
plt.savefig("multi_layer.pdf", bbox_inches="tight")

NetCDF and Climate Data with xarray

import xarray as xr
import numpy as np

# Open a single NetCDF file
ds = xr.open_dataset("climate.nc")
print(ds)                                     # dimensions, variables, coords
print(ds.data_vars)                           # available variables

# Open multiple files (e.g., monthly data)
ds = xr.open_mfdataset("data_*.nc", combine="by_coords")

# Select by coordinates
temp = ds["temperature"]
subset = temp.sel(lat=slice(30, 50), lon=slice(-100, -70))
single_time = temp.sel(time="2020-06-15", method="nearest")

# Time series operations
annual_mean = ds.groupby("time.year").mean(dim="time")
monthly_clim = ds.groupby("time.month").mean(dim="time")
seasonal = ds.resample(time="QS-DEC").mean()

# Spatial mean
area_avg = temp.mean(dim=["lat", "lon"])

# Save processed data
ds_subset.to_netcdf("processed_output.nc")

Combining xarray with GeoPandas

# Extract xarray values at point locations
import xarray as xr
import geopandas as gpd

ds = xr.open_dataset("climate.nc")
stations = gpd.read_file("stations.geojson")

# Sample raster values at station coordinates
for idx, row in stations.iterrows():
    val = ds["temperature"].sel(
        lat=row.geometry.y, lon=row.geometry.x, method="nearest"
    ).values
    stations.loc[idx, "temperature"] = float(val)

Best Practices

1. Always check and set CRS before spatial operations; mismatched CRS cause errors.
2. Reproject to an equal-area CRS (e.g., Mollweide) before area calculations.
3. Use projected CRS (meters) for buffer and distance operations, not EPSG:4326.
4. Use xr.open_mfdataset() with chunks= for large multi-file climate datasets.
5. Close datasets with ds.close() or use context managers for large files.
6. Prefer gpd.sjoin_nearest() over manual distance loops for point matching.
7. Use .dissolve() instead of manual groupby for merging geometries.
8. Write intermediate results to GeoPackage (.gpkg) for multi-layer support.