Calculating Wind Shear Coefficients with Python #

RuntimeWarning: invalid value encountered in log followed by an all-NaN shear map — or a MemoryError during xarray broadcasting — is the failure signature this page exists to eliminate. It breaks the hub-height scaling stage of the Wind Speed & Direction Modeling workflow: the moment a paired set of height rasters is reduced to a power-law exponent (α), zero-wind cells, mismatched projections, and unchunked temporal stacks turn a one-line calculation into a corrupted resource assessment. Vertical wind speed extrapolation feeds turbine hub-height estimates and every downstream energy-yield number a lender treats as ground truth, so a silently wrong α propagates straight into capacity-factor projections without ever raising a hard error.

The power-law exponent is derived from two wind speeds measured at two heights:

The arithmetic is trivial. The production failures are not — they live in the data the formula consumes, not the formula itself.

Root-cause analysis #

Three compounding causes account for nearly every broken shear pipeline, and each maps to a distinct fix stage below:

1. Zero or negative wind speeds. ln(0) is -inf and ln(<0) is NaN. Calm periods, sensor dropouts, and nodata sentinels (often -9999) enter the logarithm directly. Aggregated across hourly or sub-hourly steps during temporal data aggregation, a single bad cell poisons the annual mean for that pixel and corrupts the downstream wake model.
2. CRS and affine divergence. Height layers sourced from different reanalysis products or LiDAR campaigns frequently carry mismatched projections or affine transforms. Without strict coordinate reference system alignment, xarray broadcasts on index position rather than geography, so cell (i, j) in the 80 m grid is divided by a different physical location in the 120 m grid. The result is a smooth-looking raster of physically impossible exponents (α > 0.6 or α < 0.0) that no exception flags.
3. Memory bloat from unchunked stacks. Loading a full multi-year stack into RAM before broadcasting triggers MemoryError on standard analytical workstations, particularly for 10+ years of ERA5 or WRF output at 100 m resolution.

Pre-flight validation #

Surface the root cause before the logarithm runs. The naive script below is the broken pattern — no spatial check, no zero handling, no chunking — and it is exactly what produces the silent corruption:

import numpy as np
import rioxarray  # xr.open_rasterio was removed in xarray 2022.06; use rioxarray instead

# Flawed approach: no CRS validation, no zero-handling, no chunking
v_low = rioxarray.open_rasterio("wind_80m.tif")
v_high = rioxarray.open_rasterio("wind_120m.tif")

alpha = np.log(v_high / v_low) / np.log(120 / 80)  # -inf / NaN on calm cells
alpha.to_netcdf("shear_coeff.nc")

The pre-flight function isolates which of the three causes is present, so a CI/CD run fails fast with a precise message instead of writing a poisoned NetCDF:

import numpy as np
import rioxarray
import xarray as xr


def preflight_shear_inputs(v_low: xr.DataArray, v_high: xr.DataArray) -> None:
    """Raise on the exact root cause before any logarithm is evaluated."""
    # Cause 2: CRS parity — both grids must share an identical EPSG/WKT
    if v_low.rio.crs != v_high.rio.crs:
        raise ValueError(
            f"CRS mismatch: {v_low.rio.crs} vs {v_high.rio.crs}. "
            "Reproject inputs to a common projected grid (e.g. EPSG:32610) first."
        )
    # Cause 2: affine transform parity — origin and resolution must align
    if not np.allclose(v_low.rio.transform(), v_high.rio.transform(), atol=1e-6):
        raise ValueError(
            "Affine transform divergence. Realign with rioxarray.reproject_match()."
        )
    # Cause 1: surface the prevalence of non-positive speeds rather than hiding it
    bad = int(((v_low <= 0) | (v_high <= 0)).sum())
    if bad:
        print(f"[preflight] {bad} cells have non-positive wind speed; "
              "these will route to terrain fallback, not NaN.")

Fix implementation #

The corrected function enforces safe logarithmic evaluation, applies dask chunking for out-of-core processing, and routes invalid cells to a terrain-classified fallback exponent. Parameter choices are justified for energy use: chunk_size=1024 matches common COG tile geometry, the [0.0, 0.50] clamp bounds α to physically plausible values, and nodata-driven calm cells fall back to a documented terrain default rather than NaN.

import numpy as np
import xarray as xr
import rioxarray


def compute_wind_shear_exponent(
    v_low_path: str,
    v_high_path: str,
    h_low: float,
    h_high: float,
    fallback_alpha: float = 0.14,   # open terrain; use 0.22 for forested
    chunk_size: int = 1024,
    min_alpha: float = 0.0,
    max_alpha: float = 0.50,
) -> xr.DataArray:
    """Power-law shear exponent with CRS alignment, zero-wind masking,
    out-of-core chunking, and deterministic terrain fallback routing."""
    # 1. Open lazily with dask chunks for memory-safe execution
    chunks = {"band": 1, "y": chunk_size, "x": chunk_size}
    v_low = rioxarray.open_rasterio(v_low_path, chunks=chunks)
    v_high = rioxarray.open_rasterio(v_high_path, chunks=chunks)

    # 2. Spatial validation (Cause 2) — fail before broadcasting
    preflight_shear_inputs(v_low, v_high)

    # 3. Safe logarithm (Cause 1) — mask non-positive speeds out of ln()
    valid = (v_low > 0) & (v_high > 0)
    log_ratio = np.log(v_high.where(valid) / v_low.where(valid)) / np.log(h_high / h_low)

    # 4. Fallback routing + physical bounds clamp
    alpha = xr.where(valid, log_ratio, fallback_alpha).clip(min_alpha, max_alpha)

    # 5. Audit-ready provenance metadata
    alpha.rio.write_crs(v_low.rio.crs, inplace=True)
    alpha.attrs.update({
        "method": "power_law_shear",
        "fallback_exponent": fallback_alpha,
        "h_low_m": h_low,
        "h_high_m": h_high,
        "alpha_bounds": [min_alpha, max_alpha],
        "processing_note": "non-positive speeds routed to terrain fallback; clipped to physical bounds",
        "spatial_validation": "CRS & affine transform parity enforced",
    })
    return alpha

Calling v.where(valid) before the division keeps the lazy dask graph intact and ensures the logarithm never sees a non-positive operand, so the RuntimeWarning disappears at its source rather than being suppressed after the fact.

Fallback routing & performance tuning #

For continental-scale or CI/CD runs where the full stack will not fit in RAM, layer these strategies on top of the core function:

• Match chunk geometry to tile size. Align chunk_size with the raster’s internal tiling (typically 256×256 or 512×512). Mismatched chunks force dask to re-block, multiplying I/O and latency.
• Compress on write. Persist with encoding={"zlib": True, "complevel": 4} to balance throughput against storage footprint on multi-terabyte shear archives.
• Stay lazy. Avoid .load(), .values, or .compute() until the final aggregation. Let xarray defer execution until .to_netcdf() so dask can fuse the logarithm, mask, and clip into one pass.
• Slice time sequentially. Process annual or seasonal slabs one at a time instead of broadcasting a full multi-year array; ingest with xr.open_mfdataset(..., parallel=True).
• Mask non-representative terrain. Power-law assumptions break over coastlines and complex relief, so clip with a validated land-use polygon — the same boundary discipline used in terrain shadow analysis pipelines — and let those cells inherit the documented fallback rather than a spurious extrapolation.

Downstream validation #

Before a shear map feeds a yield model, gate it with an assertion function suitable for a CI/CD pipeline. This catches nodata bleed, out-of-bounds exponents, and CRS loss introduced by an upstream regression:

def assert_shear_integrity(alpha: xr.DataArray, max_fallback_frac: float = 0.25) -> None:
    """CI/CD gate: fail the build if the shear map is not assessment-grade."""
    assert alpha.rio.crs is not None, "output lost its CRS"
    finite = alpha.where(np.isfinite(alpha))
    assert float(finite.min()) >= 0.0, "negative shear exponent present"
    assert float(finite.max()) <= 0.50, "non-physical shear exponent (>0.50)"
    assert int(np.isnan(alpha).sum()) == 0, "NaN bleed into shear output"

    fb = alpha.attrs.get("fallback_exponent")
    frac = float((alpha == fb).mean())
    assert frac <= max_fallback_frac, (
        f"{frac:.0%} of cells routed to fallback (> {max_fallback_frac:.0%}); "
        "inputs are too sparse to be defensible"
    )

Logging the fallback fraction as part of the provenance trail is what keeps the assessment auditable: an independent engineer reviewing the interconnection or project-finance package can see exactly how many cells were extrapolated versus measured, mirroring the metadata discipline enforced in spatial data quality validation. Pin xarray, rioxarray, and dask versions in pyproject.toml so a default-broadcasting change cannot silently shift the map between runs.

Related #

• Wind Speed & Direction Modeling — parent workflow for the hub-height field this exponent scales
• Temporal Data Aggregation — reduce hourly speeds to the annual layers shear consumes
• Terrain Shadow Analysis Pipelines — terrain masking and complex-relief boundary handling
• Coordinate Reference Systems for Energy Projects — projected-CRS enforcement before any raster math