14 min read 3187 words

Wind Speed & Direction Modeling

Wind speed and direction modeling sits at the analytical core of the Solar & Wind Resource Modeling Workflows pipeline, and it fails in a way that scalar resource modeling never does. The specific failure mode this workflow exists to eliminate is the 0°/360° directional discontinuity: wind direction is a circular quantity, so interpolating bearings directly — averaging a 350° reading with a 10° reading and getting 180° instead of 0° — produces a vector field that points the wrong way across half the domain. A naive script does not raise an error. It returns a smooth-looking raster of completely wrong directions, the wind rose rotates, the wake model places turbine deficits in the wrong cells, and the energy yield that a lender treats as ground truth inherits a systematic bias no downstream step can detect.

Two further failure modes compound the first. Distance-based interpolation run in geographic coordinates (EPSG:4326) weights stations by degrees rather than metres, so a kilometre of east-west separation at 50° latitude counts for roughly two-thirds of the same distance north-south — the field is stretched before a single physical calculation runs. And the regular grids that bankable wind atlases demand are large enough that an unchunked scipy interpolation call materializes the full coordinate stack in RAM and triggers the out-of-memory reaper precisely on the continental runs that matter. This page builds a deterministic workflow that turns raw anemometer, mast, and LiDAR observations into an analysis-ready hub-height wind field: vectors are decomposed into orthogonal components before any interpolation, every input is forced into a projected metric frame, the grid is filled in bounded memory chunks, speeds are scaled to turbine hub height with an explicit shear law, and every output band carries the provenance an interconnection or project-finance review needs.

Why naive directional interpolation fails

Wind is a vector, but most station feeds report it as two scalars: a speed wind_speed_ms and a meteorological bearing wind_dir_deg measured clockwise from north as the direction the wind blows from. The temptation is to interpolate those two scalars independently onto the target grid. Speed interpolates cleanly because it is a true magnitude. Direction does not, because the number line wraps: 359° and 1° are two degrees apart physically but 358 degrees apart numerically. Any interpolator that treats bearing as an ordinary real number — linear, nearest-neighbour, kriging, or IDW — produces a spurious gradient wherever the field crosses north, and that boundary almost always runs straight through the prevailing-wind sector of a real site.

The fix is to decompose each observation into orthogonal U (west-to-east) and V (south-to-north) components, interpolate each component independently as an ordinary continuous field, and reconstruct speed and direction from the gridded components. Using the meteorological from-convention the components are:

and the inverse reconstruction, with the four-quadrant arctangent and a wrap into the [0, 360) range, is:

Because atan2 consumes the signed u and v directly, the discontinuity never enters the arithmetic — the wrap is applied once, at the very end, on the reconstructed bearing rather than on every interpolation weight. This is the same vector-first discipline that the child workflow on calculating wind shear coefficients with Python applies to the vertical profile.

Why bearings must be decomposed before interpolation Two side-by-side panels each take the same two wind directions, 350 degrees and 10 degrees, which lie two degrees apart around north. The left panel takes the scalar mean of the bearings and gets 180 degrees, a compass arrow pointing south, labelled wrong because the number line wraps at 360. The right panel decomposes each reading into orthogonal U and V components, interpolates the two components independently, then reconstructs the bearing with atan2 and gets 0 degrees, a compass arrow pointing north, labelled correct because the wrap never enters the arithmetic. Naive — interpolate the bearings WRONG Reading A 350° Reading B 10° scalar mean of bearings (350 + 10) / 2 = 180° number line wraps across 0° / 360° N points South Correct — decompose to U / V first CORRECT 350° → u, v u=−V·sinθ, v=−V·cosθ 10° → u, v u=−V·sinθ, v=−V·cosθ interpolate u and v separately atan2(v, u) → 0° wrap applied once, at the end only S points North

Prerequisites & data requirements

This workflow assumes a tabular station inventory and a projected analysis grid. Concretely:

  • Inputs: a CSV or Parquet table of meteorological stations carrying station_id, latitude, longitude, wind_speed_ms, and wind_dir_deg. Sourcing these from versioned, machine-readable open energy data portals keeps provenance and licensing explicit when the same artifacts later feed a permitting submission.
  • Coordinate frames: stations arrive in geographic coordinates (EPSG:4326). All interpolation runs in a projected metric frame — a UTM zone such as EPSG:32612 (UTM Zone 12N) or a regional Albers — chosen so that distances are preserved. Picking the right projection is the subject of coordinate reference systems for energy projects; for wind layout work a conformal UTM zone preserves the local angles that direction and terrain channeling depend on.
  • Geometry hygiene: invalid or out-of-range records must be quarantined before interpolation, applying the same spatial data quality validation gates the rest of the pipeline relies on.
  • Library versions: geopandas ≥ 0.14, pyproj ≥ 3.6, scipy ≥ 1.11, and rasterio ≥ 1.3. Pin them, because scipy.interpolate and rasterio occasionally change default behaviour across releases.

CRS harmonization & station ingestion

Spatial interpolation degrades rapidly when distances are computed in degrees. The first step standardizes every input into the target projected CRS and rejects records that would poison the field — negative speeds, bearings outside [0, 360], coordinates off the globe — before any geometry is constructed. The ingest gate is the cheapest place to catch these errors.

python
import geopandas as gpd
import pandas as pd
import numpy as np
import logging
from pathlib import Path
from pyproj import CRS

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")

def prepare_station_gdf(csv_path: str | Path, target_epsg: int = 32612) -> gpd.GeoDataFrame:
    """
    Ingest a wind station CSV, validate coordinates and observations, and
    transform to a projected metric CRS. Filtering runs before geometry
    construction to keep the memory footprint flat on large station feeds.
    """
    path = Path(csv_path)
    if not path.exists():
        raise FileNotFoundError(f"Station data not found: {path}")

    df = pd.read_csv(path, dtype={"station_id": str})
    required_cols = {"station_id", "latitude", "longitude", "wind_speed_ms", "wind_dir_deg"}
    missing = required_cols - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")

    # Spatial + physical validation: drop records that would corrupt the field
    valid_mask = (
        df["latitude"].between(-90, 90) &
        df["longitude"].between(-180, 180) &
        df["wind_speed_ms"].ge(0) &
        df["wind_dir_deg"].between(0, 360)
    )
    dropped = int((~valid_mask).sum())
    df = df.loc[valid_mask].copy()
    logging.info(f"Retained {len(df)} stations; quarantined {dropped} invalid records.")

    station_gdf = gpd.GeoDataFrame(
        df,
        geometry=gpd.points_from_xy(df.longitude, df.latitude),
        crs=CRS.from_epsg(4326),
    )
    station_proj = station_gdf.to_crs(epsg=target_epsg)
    logging.info(f"Transformed to EPSG:{target_epsg} | bounds: {station_proj.total_bounds}")
    return station_proj

The to_crs call is explicit and logged — never an implicit on-the-fly reprojection buried inside an analysis routine — so the EPSG decision that governs every later distance is recoverable from the run log.

Core implementation: vectorized U/V interpolation

With clean, projected stations in hand, the happy-path workflow decomposes the vectors, interpolates each component onto a regular grid, and returns the gridded U and V arrays plus the affine metadata that downstream rasterization needs. Direction never touches the interpolator; only its sine and cosine projections do.

python
import rasterio
from scipy.interpolate import griddata
from typing import Tuple

def decompose_and_interpolate(
    station_gdf: gpd.GeoDataFrame,
    grid_res_m: float = 100.0,
    method: str = "linear",
    chunk_size: int = 10_000,
) -> Tuple[np.ndarray, np.ndarray, dict]:
    """
    Decompose wind vectors into U/V, interpolate each component onto a regular
    grid, and return the gridded components plus raster metadata. Avoids the
    0/360 discontinuity by interpolating components, never bearings.
    """
    coords = np.column_stack((station_gdf.geometry.x, station_gdf.geometry.y))
    speeds = station_gdf["wind_speed_ms"].to_numpy(dtype="float64")
    dirs_rad = np.deg2rad(station_gdf["wind_dir_deg"].to_numpy(dtype="float64"))

    # Vector decomposition (meteorological "from" convention)
    u = -speeds * np.sin(dirs_rad)
    v = -speeds * np.cos(dirs_rad)

    # Regular grid spanning the station bounds at the requested resolution
    minx, miny, maxx, maxy = station_gdf.total_bounds
    cols = int(np.ceil((maxx - minx) / grid_res_m))
    rows = int(np.ceil((maxy - miny) / grid_res_m))
    xi = np.linspace(minx, maxx, cols)
    yi = np.linspace(miny, maxy, rows)
    grid_x, grid_y = np.meshgrid(xi, yi)

    # Chunked interpolation: fill the grid in bounded blocks so a continental
    # domain never materializes the full coordinate stack in RAM at once.
    u_grid = np.full(grid_x.shape, np.nan, dtype="float32")
    v_grid = np.full(grid_x.shape, np.nan, dtype="float32")
    flat = np.column_stack((grid_x.ravel(), grid_y.ravel()))
    for i in range(0, len(flat), chunk_size):
        block = flat[i:i + chunk_size]
        u_grid.ravel()[i:i + chunk_size] = griddata(coords, u, block, method=method)
        v_grid.ravel()[i:i + chunk_size] = griddata(coords, v, block, method=method)

    metadata = {
        "transform": rasterio.transform.from_origin(minx, maxy, grid_res_m, grid_res_m),
        "crs": station_gdf.crs,
        "shape": (rows, cols),
        "bounds": (minx, miny, maxx, maxy),
        "grid_res_m": grid_res_m,
    }
    return u_grid, v_grid, metadata


def reconstruct_speed_direction(
    u_grid: np.ndarray, v_grid: np.ndarray
) -> Tuple[np.ndarray, np.ndarray]:
    """Rebuild scalar speed and meteorological bearing from gridded components."""
    speed = np.hypot(u_grid, v_grid)                       # sqrt(u^2 + v^2)
    bearing = (270.0 - np.degrees(np.arctan2(v_grid, u_grid))) % 360.0
    return speed.astype("float32"), bearing.astype("float32")

Keeping the working dtype at float32 halves memory versus float64 with negligible loss for a wind field, and np.hypot avoids the intermediate overflow that a literal sqrt(u**2 + v**2) can hit on extreme gusts.

Hub-height extrapolation & wind shear

Turbine hub heights routinely exceed the mast or LiDAR measurement elevation, so the gridded components must be scaled vertically. The power-law profile relates speed at height to the reference speed via the shear exponent :

Because and scale linearly with speed, the same ratio applies to both components, so the field can be scaled by broadcasting a single scalar across the grid. The exponent itself varies with terrain roughness and atmospheric stability; deriving a defensible, site-specific rather than assuming the open-terrain default of 0.143 is the subject of the companion workflow on calculating wind shear coefficients with Python.

python
def apply_hub_height_scaling(
    u_grid: np.ndarray,
    v_grid: np.ndarray,
    alpha: float = 0.143,
    meas_height_m: float = 50.0,
    hub_height_m: float = 100.0,
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Apply the power-law shear ratio to U/V grids. NaN cells (outside the
    convex hull of the stations) are preserved as nodata, not scaled.
    """
    if not 0.0 <= alpha <= 0.5:
        logging.warning("Shear exponent %.3f outside typical [0.0, 0.5]; verify site.", alpha)

    ratio = (hub_height_m / meas_height_m) ** alpha
    u_scaled = np.where(np.isnan(u_grid), np.nan, u_grid * ratio).astype("float32")
    v_scaled = np.where(np.isnan(v_grid), np.nan, v_grid * ratio).astype("float32")
    return u_scaled, v_scaled

Local topographic acceleration and flow channeling are not captured by a uniform shear ratio; correcting for them means debiting the field with the slope and aspect masks produced by terrain and shadow analysis pipelines before the field is treated as final.

Error handling & edge cases

The three failure modes named in the problem framing each need an explicit guard rather than a hopeful assumption.

1. Directional discontinuity leaking back in. The decomposition only protects the field if nothing downstream re-interpolates the reconstructed bearing. Guard against accidental scalar handling by asserting that any directional aggregation goes through the components:

python
def circular_mean_deg(dirs_deg: np.ndarray) -> float:
    """Correct mean bearing via unit-vector averaging — never a scalar mean()."""
    rad = np.deg2rad(dirs_deg)
    s, c = np.nanmean(np.sin(rad)), np.nanmean(np.cos(rad))
    if np.hypot(s, c) < 1e-9:
        return float("nan")          # directionless: cancelling vectors
    return float((np.degrees(np.arctan2(s, c))) % 360.0)

2. CRS mismatch or a geographic grid. If the station GeoDataFrame is still in EPSG:4326 when it reaches the interpolator, every distance weight is wrong. Fail loudly instead of producing a stretched field:

python
def assert_projected(station_gdf: gpd.GeoDataFrame) -> None:
    crs = station_gdf.crs
    if crs is None:
        raise ValueError("Station CRS is undefined; refuse to interpolate.")
    if crs.is_geographic:
        raise ValueError(
            f"Interpolation requires a projected CRS; got geographic {crs.to_epsg()}. "
            "Reproject to a UTM zone (e.g. EPSG:32612) first."
        )

3. Sparse stations and empty grids. When too few stations survive the ingest gate, griddata returns an all-NaN grid for method="linear" (which only fills the convex hull). Detect the degenerate case and either fall back to nearest for the extrapolation margin or abort with a clear message:

python
def guard_station_density(station_gdf: gpd.GeoDataFrame, min_stations: int = 4) -> None:
    if len(station_gdf) < min_stations:
        raise ValueError(
            f"Only {len(station_gdf)} valid stations; need ≥ {min_stations} for a "
            "defensible linear interpolation. Widen the catchment or use nearest."
        )

Performance & scalability: async rasterization

For a continental wind atlas the interpolation is CPU-bound and the serialization is I/O-bound, and the two should not block each other. The chunked grid fill already bounds interpolation memory; the write side benefits from offloading the GeoTIFF serialization so the event loop stays responsive during large tiled writes. Note that a rasterio DatasetWriter is not safe to share across threads — the entire write loop is offloaded to a single executor thread rather than fanning blocks across a pool.

python
import asyncio
from concurrent.futures import ThreadPoolExecutor
from rasterio.windows import Window

async def write_wind_raster_async(
    u_grid: np.ndarray,
    v_grid: np.ndarray,
    metadata: dict,
    output_path: str | Path,
    block_size: int = 512,
) -> None:
    """Write U/V bands to a single tiled, compressed GeoTIFF off the event loop."""
    out_path = Path(output_path)
    out_path.parent.mkdir(parents=True, exist_ok=True)

    if not np.isfinite(u_grid).any() or not np.isfinite(v_grid).any():
        raise ValueError("Grid has no valid cells; check station density and bounds.")

    profile = {
        "driver": "GTiff",
        "dtype": "float32",
        "count": 2,
        "width": metadata["shape"][1],
        "height": metadata["shape"][0],
        "crs": metadata["crs"],
        "transform": metadata["transform"],
        "compress": "deflate",
        "nodata": np.nan,
        "blockxsize": block_size,
        "blockysize": block_size,
        "tiled": True,
    }

    def _write() -> None:
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.set_band_description(1, "wind_speed_u_ms")
            dst.set_band_description(2, "wind_speed_v_ms")
            for row in range(0, profile["height"], block_size):
                for col in range(0, profile["width"], block_size):
                    win = Window(col, row,
                                 min(block_size, profile["width"] - col),
                                 min(block_size, profile["height"] - row))
                    u_block = u_grid[win.row_off:win.row_off + win.height,
                                     win.col_off:win.col_off + win.width]
                    v_block = v_grid[win.row_off:win.row_off + win.height,
                                     win.col_off:win.col_off + win.width]
                    dst.write(np.stack([u_block, v_block]), indexes=[1, 2], window=win)

    loop = asyncio.get_running_loop()
    with ThreadPoolExecutor(max_workers=1) as executor:
        await loop.run_in_executor(executor, _write)
    logging.info("Async rasterization complete: %s", out_path)

Beyond a single grid, the usual scaling levers apply: align block_size with the GeoTIFF tile dimensions to avoid re-blocking on read, raise GDAL_CACHEMAX for write-heavy runs, and for very large domains build per-tile interpolations behind a VRT rather than one monolithic array. Most bottlenecks here come from a redundant reprojection inside a loop or an unchunked grid fill — profile before reaching for a bigger machine.

Wind Speed and Direction Modeling Pipeline A six-stage left-to-right data flow. Stage 1 ingests a station CSV or Parquet table of anemometer, mast and LiDAR observations. Stage 2 is the CRS gate, which validates records and reprojects from EPSG:4326 to a metric frame such as EPSG:32612. Stage 3 decomposes each vector into U and V components using u equals minus V sine theta and v equals minus V cosine theta. Stage 4 interpolates each component independently with chunked griddata in bounded memory; a downward branch off this stage detects an all-NaN or under-four-station grid and routes it to an abort-and-log node so no corrupt field is serialized. Stage 5 applies the power-law hub-height shear ratio. The highlighted terminal Stage 6 writes a tiled, compressed GeoTIFF carrying U and V bands plus hub-height and alpha provenance tags. 1 · Stations CSV / Parquet anemometer · LiDAR 2 · CRS Gate validate · reproject 4326 → 32612 3 · Decompose u = −V·sinθ v = −V·cosθ 4 · Interpolate chunked griddata per component 5 · Hub Scaling power-law α (z / z_ref)^α 6 · GeoTIFF async tiled write U/V bands · tags degenerate guard: all-NaN / < 4 stations abort run + log reason no corrupt field written

Validation & audit trail

A wind field is only bankable if its integrity is asserted, not assumed. Final outputs should conform to CF-Conventions and OGC GeoTIFF expectations, mirroring the quality gates applied in solar irradiance raster processing so the two technologies stay interoperable within a hybrid portfolio. The audit function below is suitable for a CI/CD gate that blocks a release when output integrity regresses.

python
import json

def audit_wind_raster(raster_path: str, expected_epsg: int,
                      climatology_p95_ms: float = 35.0) -> dict:
    """Assert band count, dtype, CRS, and physical sanity of a wind GeoTIFF."""
    with rasterio.open(raster_path) as src:
        report = {
            "band_count": src.count,
            "dtype": src.dtypes[0],
            "crs_epsg": src.crs.to_epsg() if src.crs else None,
            "descriptions": list(src.descriptions),
        }
        u = src.read(1, masked=True)
        v = src.read(2, masked=True)

    speed = np.ma.sqrt(u**2 + v**2)
    report["max_speed_ms"] = float(speed.max())

    assert report["band_count"] == 2, "Expected U and V bands"
    assert report["dtype"] == "float32", f"Want float32, got {report['dtype']}"
    assert report["crs_epsg"] == expected_epsg, f"CRS drift: {report['crs_epsg']}"
    # Statistical sanity: reconstructed speed must not exceed regional climatology
    assert report["max_speed_ms"] < climatology_p95_ms, "Speed exceeds climatology p95"

    logging.info("Audit passed: %s", json.dumps(report))
    return report

The completeness checklist for a compliant artifact is: band descriptions (wind_speed_u_ms, wind_speed_v_ms) plus the hub_height_m and shear alpha recorded as raster tags; projection consistency with the project boundary shapefile; and the statistical sanity check that reconstructed speed stays below the 95th percentile of regional climatology. Production deployments wrap the whole sequence in a configuration-driven orchestrator (Prefect or Airflow) for temporal aggregation, chunk-level retries, and metadata cataloging. Those gridded fields then feed grid-screening work, where the resource surface is cross-referenced against grid capacity buffer analysis thresholds — the metadata contract is what lets two pipelines trust each other’s outputs.

Explore this section

Back to Solar & Wind Resource Modeling Workflows