geopandas-spatial

GeoPandas Spatial Analysis

Geospatial data analysis using geopandas for vector data and xarray for multidimensional climate/weather datasets.

When to Use

• Loading and analyzing shapefiles, GeoJSON, GeoPackage
• Spatial joins, intersections, buffers, and dissolves
• Climate and weather data from NetCDF files
• CRS transformations and geographic projections
• Map visualization and choropleth maps
• Area, distance, and geometric calculations

When NOT to Use

• Satellite imagery classification or ML (use rasterio/torchgeo)
• Real-time GPS tracking or routing
• Interactive web map applications (use folium or deck.gl)
• General tabular data without spatial component (use pandas)

Reading Vector Data

import geopandas as gpd

# Read various vector formats
gdf = gpd.read_file("boundaries.shp")
gdf = gpd.read_file("data.geojson")
gdf = gpd.read_file("database.gpkg", layer="cities")

# Read from URL
gdf = gpd.read_file("https://example.com/regions.geojson")

# Inspect the GeoDataFrame
print(gdf.head())
print(gdf.crs)                    # coordinate reference system
print(gdf.geometry.type.unique()) # geometry types present
print(gdf.total_bounds)           # [minx, miny, maxx, maxy]

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")