Terrain & Shadow Analysis Pipelines #

Terrain and shadow analysis is the validation layer that decides whether a resource estimate survives contact with the actual ground. Mesoscale models and satellite irradiance products assume an unobstructed sky; real sites sit in valleys, behind ridgelines, and on north-facing slopes that clip the morning and evening sun and steer the wind. This workflow is the topographic-correction stage of the broader Solar & Wind Resource Modeling Workflows pipeline: it takes a digital elevation model (DEM) and a time series of solar positions and produces the shadow-loss and slope-constraint surfaces that turn a flat-sky resource figure into a defensible, terrain-aware yield projection. The specific failure mode this stage exists to eliminate is horizon drift in misaligned DEM stacks — when the elevation grid and the irradiance grid disagree on projection, pixel registration, or vertical datum by even a fraction of a cell, every cast shadow lands in the wrong place and the resulting terrain loss is biased in a way no later step can detect.

The goal is deterministic: convert raw elevation into a binary or fractional shadow mask, plus slope and aspect derivatives, that align pixel-for-pixel with the irradiance surface they will modulate, carry explicit provenance, and quantify terrain-induced losses before the financial model consumes them. This page covers the conceptual foundation, the prerequisites, a full runnable horizon-and-shadow function, the failure modes that break naive shadow casting, the scalability patterns for high-resolution lidar-derived DEMs, and the audit trail that makes a terrain-loss figure bankable in permitting and project-finance review.

Why naive shadow casting fails #

The intuitive approach — loop over every pixel, march a ray toward the sun, and flag the pixel as shaded the moment a higher cell appears — fails on two independent axes: spatial correctness and computational cost. Both compound into the silent bias this stage must prevent.

First, spatial misregistration between the DEM and the resource grid. Shadow masks are not consumed in isolation; they multiply the direct component of an irradiance surface produced upstream by solar irradiance raster processing. If the DEM is delivered in geographic coordinates (EPSG:4326) with degree-based spacing while the irradiance grid is in a projected metric system such as EPSG:32612, the horizon angles computed from degree distances are wrong by a latitude-dependent factor, and the mask is registered half a cell — tens of metres on the ground — away from the irradiance pixels it is supposed to dim. Enforcing coordinate reference system alignment into one projected, metric target is the precondition for every subsequent angle calculation.

Second, vertical datum and unit mismatch. Horizon elevation is an angle built from a rise over a run. If the horizontal run is in metres but the vertical rise is in feet, or the DEM mixes an ellipsoidal height with an orthometric (geoid) reference partway through a mosaic, the computed terrain angle is systematically wrong. A 30% vertical scaling error from a foot/metre confusion turns a true 4° horizon into a 5.2° horizon and over-reports morning shadow loss across the whole site.

Third, per-pixel ray marching does not scale. A naive nested loop is $O(N\times T\times R)$ — for $N$ pixels, $T$ timestamps, and $R$ steps along each ray. A 4000×4000 lidar tile evaluated hourly across a year is on the order of $10^{13}$ ray steps, which is why production pipelines precompute a per-cell horizon profile once and reuse it for every timestamp. The angle to the local horizon in a given azimuth direction does not change with time; only the sun moves. Separating the time-invariant horizon profile from the time-varying solar position collapses the cost by orders of magnitude.

A clean pipeline therefore decouples spatial validation, horizon profiling, and temporal shadow evaluation into discrete, testable stages so that projection drift, datum mismatch, or registration error is caught and rejected before any shadow is cast, rather than discovered after the yield model has already absorbed the bias.

Prerequisites and data requirements #

Before running the workflow, pin the inputs and the environment so terrain results are reproducible across a portfolio:

• Library versions: rasterio>=1.3, numpy>=1.24, pyproj>=3.5, and pvlib>=0.10 for solar position. GDAL underpins rasterio; keep it >=3.6. Datum transforms should defer to the pyproj documentation and the EPSG registry rather than hand-coded proj strings.
• Target CRS: a single projected, metric CRS matching the project’s UTM zone (for example EPSG:32612 for the US Mountain West). Store the EPSG integer, never an unqualified “UTM 12N” string, and confirm the DEM shares this CRS with the irradiance grid produced by solar irradiance raster processing.
• Input geometry: a single-band float DEM (GeoTIFF or Cloud-Optimized GeoTIFF) of ground or surface elevation in metres, with a defined CRS, a declared nodata value, and a documented vertical datum (e.g. NAVD88 orthometric). Lidar-derived DSMs capture vegetation and structures that cast real shadow; bare-earth DTMs do not — pick the one that matches the obstruction question. Surfaces missing a CRS must be rejected at ingest; see spatial data quality validation for the cleaning patterns this stage assumes upstream.
• Solar geometry: per-timestamp solar azimuth and elevation for the site latitude/longitude, typically from pvlib.solarposition. Only timestamps with positive solar elevation matter; the sun below the horizon produces a trivially fully-shaded frame.
• Source provenance: DEM acquisition date, source portal, product version, and ground sample distance, sourced from one of the documented open energy data portals so the terrain-loss output’s lineage is auditable.

A pixel is shaded at a given instant when the sun sits below the local terrain horizon in the sun’s azimuth direction. With $h(\theta )$ the horizon elevation angle toward azimuth $\theta$ and ${\alpha }_s$ the solar elevation:

The horizon angle itself, for a neighbour cell at planimetric distance $d$ and elevation difference $\Delta z$, is:

Because $d$ must be a true metric distance, this formula is only valid once the DEM is in a projected CRS — the algebraic restatement of why the spatial validation gate runs first.

Core implementation #

The loader below enforces a metric projection and a sane vertical range, then generates memory-safe windows for chunked execution. It rejects geographic CRS inputs early, because degree-based runs make the arctan horizon angle meaningless. Variable names are energy-specific throughout.

import logging
from pathlib import Path
from typing import List, Tuple

import numpy as np
import rasterio
from rasterio.windows import Window

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

TARGET_EPSG = 32612          # UTM Zone 12N — must match the irradiance grid
CHUNK_PX = 1024              # memory-safe tile dimension
PLAUSIBLE_ELEV_M = (-430.0, 8849.0)  # Dead Sea floor to Everest — sanity bounds


def validate_and_window_dem(dem_path: Path,
                            target_epsg: int = TARGET_EPSG,
                            chunk_px: int = CHUNK_PX) -> Tuple[dict, List[Window]]:
    """Validate DEM CRS, datum range, and resolution; return metadata + windows."""
    with rasterio.open(dem_path) as src:
        if src.crs is None:
            raise ValueError(f"{dem_path.name}: DEM lacks a defined CRS — assign before processing.")
        if src.crs.is_geographic:
            raise RuntimeError(
                f"{dem_path.name}: geographic CRS {src.crs} detected. Reproject to a metric "
                f"projection (EPSG:{target_epsg}) so horizon arctan runs are true distances."
            )
        if src.crs.to_epsg() != target_epsg:
            raise RuntimeError(
                f"{dem_path.name}: CRS EPSG:{src.crs.to_epsg()} != irradiance grid EPSG:{target_epsg}; "
                f"shadow mask would be misregistered against the resource surface."
            )

        dem = src.read(1, masked=True)
        lo, hi = float(dem.min()), float(dem.max())
        if lo < PLAUSIBLE_ELEV_M[0] or hi > PLAUSIBLE_ELEV_M[1]:
            raise ValueError(
                f"{dem_path.name}: elevation range [{lo:.0f}, {hi:.0f}] m implausible — "
                f"check vertical datum / unit (feet vs metres)."
            )

        windows = [
            Window(col, row,
                   min(chunk_px, src.width - col),
                   min(chunk_px, src.height - row))
            for row in range(0, src.height, chunk_px)
            for col in range(0, src.width, chunk_px)
        ]
        logging.info("%s: %d windows, native res %.1f m, EPSG:%d.",
                     dem_path.name, len(windows), abs(src.res[0]), target_epsg)
        return src.meta.copy(), windows

With validation in place, the next function precomputes the time-invariant horizon profile and evaluates the time-varying shadow state. The horizon is sampled along a set of azimuth bearings by marching outward in metric steps and tracking the cumulative maximum elevation angle — the standard “ray-sweep” formulation — and a per-timestamp mask is then a vectorized comparison of solar elevation against the horizon angle for the sun’s bearing.

def compute_horizon_profile(elevation_m: np.ndarray,
                            cell_size_m: float,
                            azimuth_deg: float,
                            max_distance_m: float = 5000.0) -> np.ndarray:
    """Time-invariant horizon angle (radians) toward one azimuth, via cumulative-max ray sweep."""
    rows, cols = elevation_m.shape
    horizon = np.zeros((rows, cols), dtype=np.float32)

    # Per-step pixel offsets along the azimuth bearing (0deg = North, clockwise).
    az = np.radians(azimuth_deg)
    step_row = -np.cos(az)   # north is negative row direction
    step_col = np.sin(az)
    n_steps = int(max_distance_m / cell_size_m)

    base_r, base_c = np.mgrid[0:rows, 0:cols]
    for step in range(1, n_steps + 1):
        sample_r = np.round(base_r + step * step_row).astype(int)
        sample_c = np.round(base_c + step * step_col).astype(int)
        inside = (sample_r >= 0) & (sample_r < rows) & (sample_c >= 0) & (sample_c < cols)
        rr = np.clip(sample_r, 0, rows - 1)
        cc = np.clip(sample_c, 0, cols - 1)

        delta_z = elevation_m[rr, cc] - elevation_m       # rise toward the neighbour
        run = step * cell_size_m                          # true metric distance
        angle = np.where(inside, np.arctan2(delta_z, run), 0.0)
        horizon = np.maximum(horizon, angle.astype(np.float32))
    return horizon


def shadow_mask_for_timestamp(horizon_angle_rad: np.ndarray,
                              solar_elevation_deg: float) -> np.ndarray:
    """Binary shadow mask: 1 where the sun sits at/below the local terrain horizon."""
    if solar_elevation_deg <= 0:
        return np.ones_like(horizon_angle_rad, dtype=np.uint8)   # sun below horizon
    solar_elev_rad = np.radians(solar_elevation_deg)
    return (solar_elev_rad <= horizon_angle_rad).astype(np.uint8)

The two functions split the cost exactly where it matters: compute_horizon_profile runs once per azimuth bin per tile (the expensive sweep), while shadow_mask_for_timestamp is a single vectorized comparison cheap enough to call for every hour of a test year.

Error handling and edge cases #

The failure modes named above need explicit, testable guards rather than a blanket try/except.

Geographic CRS or DEM/irradiance mismatch reaching the sweep. validate_and_window_dem already rejects an undefined CRS, a geographic CRS, and an EPSG that disagrees with the irradiance grid. This is the single most important guard, because a misregistered mask produces a plausible-looking but wrong terrain loss. Never let a degree-spaced DEM into the arctan horizon step.

Vertical datum / unit confusion. A foot-valued DEM tagged as metres, or an ellipsoidal-vs-orthometric splice, slips past a CRS check because the horizontal CRS is valid. Guard the vertical axis independently with a physical range test and a slope sanity check:

def assert_vertical_sanity(elevation_m: np.ndarray, cell_size_m: float) -> None:
    """Catch foot/metre and datum-splice errors before they bias horizon angles."""
    gy, gx = np.gradient(elevation_m, cell_size_m)
    max_slope_deg = float(np.degrees(np.arctan(np.hypot(gy, gx).max())))
    if max_slope_deg > 85.0:
        raise ValueError(
            f"Max slope {max_slope_deg:.0f} deg implies a vertical-unit or datum error "
            f"(feet read as metres inflates rise ~3.28x)."
        )

Nodata bleed into the horizon sweep. Voids in lidar DEMs (water bodies, occlusions) arrive as a sentinel such as -9999. Left unmasked, a single void cell injects a spurious cliff that casts a kilometre of false shadow. Replace nodata with np.nan before the sweep and treat NaN neighbours as non-occluding:

elevation_m = np.where(elevation_m == src.nodata, np.nan, elevation_m)
# In the sweep, NaN deltas yield NaN angles; np.fmax ignores them:
horizon = np.fmax(horizon, np.nan_to_num(angle, nan=-np.inf))

Edge truncation on tile borders. A ridge just outside a 1024-px window still shades pixels inside it. Process windows with an overlap halo of ceil(max_distance_m / cell_size_m) pixels and crop the halo after the sweep, so shadows cast by off-tile terrain are still captured.

Performance and scalability #

High-resolution DEMs routinely exceed RAM once paired with a multi-temporal solar-position array, so scaling is about bounding memory and overlapping I/O, not buying more of either:

• Profile once, evaluate many. Cache the per-azimuth horizon profile per tile; quantise the sun’s azimuth to a fixed set of bins (e.g. 1° or 2°) and reuse the nearest cached profile across every timestamp. This is the change that turns an intractable $O(N\times T\times R)$ sweep into a tractable one.
• Windowed reads, tiled writes, overlap halo. Read with rasterio windows, write tiled LZW-compressed output, and carry the halo above so cross-tile shadows survive chunking. Peak memory then scales with one tile plus its halo, not the whole DEM.
• Async over tiles, threads within a sweep. Coarse concurrency belongs at the tile level via asyncio and a semaphore; the same async pattern used in wind speed and direction modeling for independent temporal slices applies here, with each coroutine owning one DEM window so disk latency overlaps across cores.
• GDAL cache tuning. Set GDAL_CACHEMAX (e.g. 512) to bound the block cache during batch runs; a runaway cache, not the elevation data, is the usual cause of memory exhaustion.
• Store fractional, not just binary, masks. Aggregating hourly binary masks to a float32 mean per cell yields a fractional shaded-time surface that modulates the direct beam smoothly. Once aligned, these surfaces feed straight into temporal data aggregation for monthly and seasonal terrain-loss reduction.

The async orchestration below dispatches one coroutine per DEM window, computes the temporal shadow stack, and streams a fractional shaded-time tile to disk.

import asyncio
from concurrent.futures import ThreadPoolExecutor


async def run_async_shadow_pipeline(dem_path: Path,
                                    windows: List[Window],
                                    timestamps,                 # iterable of (azimuth_deg, elev_deg)
                                    out_dir: Path,
                                    max_concurrency: int = 3) -> None:
    """Cast shadows tile-by-tile, overlapping disk I/O with the CPU-bound sweep."""
    out_dir.mkdir(parents=True, exist_ok=True)
    loop = asyncio.get_running_loop()
    semaphore = asyncio.Semaphore(max_concurrency)

    def _process(window: Window, idx: int) -> None:
        with rasterio.open(dem_path) as src:
            elevation_m = src.read(1, window=window).astype(np.float32)
            elevation_m = np.where(elevation_m == src.nodata, np.nan, elevation_m)
            transform = src.window_transform(window)
            cell_size_m = abs(src.res[0])

        accum = np.zeros(elevation_m.shape, dtype=np.float32)
        cache: dict = {}
        for az_deg, elev_deg in timestamps:
            az_bin = round(az_deg)                       # 1-degree azimuth cache key
            if az_bin not in cache:
                cache[az_bin] = compute_horizon_profile(elevation_m, cell_size_m, az_bin)
            accum += shadow_mask_for_timestamp(cache[az_bin], elev_deg)
        shaded_fraction = (accum / max(len(list(timestamps)), 1)).astype(np.float32)

        with rasterio.open(out_dir / f"shaded_{idx}.tif", "w", driver="GTiff",
                           height=window.height, width=window.width, count=1,
                           dtype="float32", crs=f"EPSG:{TARGET_EPSG}", transform=transform,
                           nodata=np.nan, tiled=True, blockxsize=256, blockysize=256,
                           compress="lzw") as dst:
            dst.write(shaded_fraction, 1)

    async def _bounded(window: Window, idx: int) -> None:
        async with semaphore:
            with ThreadPoolExecutor(max_workers=1) as ex:
                await loop.run_in_executor(ex, _process, window, idx)

    await asyncio.gather(*(_bounded(w, i) for i, w in enumerate(windows)))

Validation and audit trail #

A terrain-loss figure is only bankable if its integrity is asserted automatically and its provenance is embedded. Every output should pass a post-processing gate suitable for a CI/CD pipeline: CRS and pixel-alignment verification against the irradiance grid, value-range checks on the fractional mask, and embedded solar-geometry and DEM-source metadata.

def assert_shadow_integrity(shaded_path: Path,
                            target_epsg: int = TARGET_EPSG) -> None:
    """CI/CD gate: fail the build if a shadow surface is non-compliant."""
    with rasterio.open(shaded_path) as out:
        assert out.crs.to_epsg() == target_epsg, "Shadow mask CRS not aligned to irradiance grid."
        assert out.dtypes[0] == "float32", "Unexpected dtype; expected float32 fractional mask."
        frac = out.read(1, masked=True)
        assert frac.count() > 0, "No valid pixels — possible disjoint extent or all-nodata tile."
        assert 0.0 <= float(np.ma.min(frac)) and float(np.ma.max(frac)) <= 1.0, \
            "Fractional shaded time outside [0, 1] — aggregation or nodata defect."
    with rasterio.open(shaded_path, "r+") as out:
        out.update_tags(
            DEM_SOURCE="USGS 3DEP lidar DTM",
            VERTICAL_DATUM="NAVD88",
            SOLAR_MODEL="pvlib.solarposition",
            HORIZON_MAX_DIST_M="5000",
            CRS_EPSG=str(target_epsg),
            QA_STATUS="passed",
        )

Pixel alignment is the non-negotiable invariant: a shadow mask is only meaningful when it multiplies the same pixels of the irradiance surface it was built to dim. Enforce an explicit tolerance — for example ±0.5 m for UTM-projected assets — and log the DEM-to-irradiance affine residual so the terrain loss is auditable for regulatory submission and project-finance due diligence. Embedding the DEM source, vertical datum, solar model, and horizon search radius via update_tags keeps the lineage attached to the file rather than to a notebook that will not survive the project. With alignment proven and provenance written, these surfaces become unified terrain-constraint layers — combining shadow loss with the slope and aspect derivatives detailed in Automating hillshade and slope analysis for wind turbine siting — that feed layout optimization and the proximity screens in grid capacity buffer analysis. For windowed raster I/O specifics, consult the Rasterio documentation on windowed reading and writing; for the async orchestration, the asyncio task scheduling guide.

Related #

• Solar & Wind Resource Modeling Workflows — the parent pipeline this terrain-correction stage feeds.
• Automating Hillshade and Slope Analysis for Wind Turbine Siting — the slope/aspect derivatives that combine with shadow loss into terrain-constraint layers.
• Solar Irradiance Raster Processing — the irradiance grid that shadow masks must align to pixel-for-pixel.
• Wind Speed & Direction Modeling — shares the async windowed-evaluation pattern for terrain-aware wind fields.
• Temporal Data Aggregation — reducing hourly shadow stacks to monthly and seasonal terrain-loss statistics.
• Coordinate Reference Systems for Energy Projects — the projection and datum foundations this stage enforces.