geopandas-spatial

Files (1)

• SKILL.md

SKILL.md (first lines)

# 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

```python
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

```python
# 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

```python
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

```python
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

```python
import xarray as xr
import numpy as np

From beita6969/ScienceClaw, MIT.

More from beita6969/ScienceClaw