Automating hillshade and slope analysis for wind turbine siting #

When scaling terrain preprocessing for utility-scale wind development, automated hillshade and slope extraction frequently fails at the intersection of memory limits, coordinate reference system (CRS) unit mismatches, and windowed raster boundary discontinuities. The symptom is consistent: pipelines execute successfully on test extents but crash with MemoryError or produce artificial 90° slope cliffs at tile seams when processing regional 1m LiDAR or 30m SRTM DEMs covering 500+ km². This guide isolates the root causes, provides a minimal reproducible failure, and delivers a production-grade, chunked pipeline with explicit CRS validation, memory-safe routing, and compliance-ready validation steps.

Root-Cause Diagnostics #

Three interdependent failures dominate terrain preprocessing pipelines in renewable energy workflows:

1. CRS Unit Mismatch: Geographic CRS (e.g., EPSG:4326) expresses horizontal distances in degrees while vertical elevation remains in meters. Gradient calculations (dz/dx, dz/dy) become mathematically invalid, yielding slope values that collapse to near-zero or spike to 90° depending on latitude. Gradient derivatives assume uniform horizontal spacing, which geographic projections violate.
2. Windowed Boundary Artifacts: Derivative filters require neighboring pixels. Standard rasterio windowing without padding creates hard seams where slope gradients drop to zero or invert, directly impacting micro-siting algorithms and access road routing within Terrain & Shadow Analysis Pipelines.
3. Memory Fragmentation: Loading full DEM arrays into memory before applying numpy.gradient or scipy.ndimage filters exhausts RAM on standard CI/CD runners or cloud instances, causing silent kernel kills or MemoryError exceptions. Unmanaged array duplication during intermediate calculations compounds the issue.

Minimal Reproducible Failure #

The following pattern demonstrates the exact failure mode observed in production environments:

import rasterio
import numpy as np

# Fails on large DEMs: loads entire array, ignores CRS units, no window padding
with rasterio.open("dem_1m.tif") as src:
    dem = src.read(1)
    # Assumes 1:1 horizontal/vertical units (false for EPSG:4326)
    grad_y, grad_x = np.gradient(dem.astype(np.float32))
    slope = np.arctan(np.sqrt(grad_x**2 + grad_y**2)) * (180 / np.pi)

This approach produces three immediate defects:

• MemoryError on DEMs > 4GB due to unchunked array allocation
• Slope artifacts at tile boundaries if chunked without overlap
• Incorrect slope magnitudes when horizontal units are degrees

Production-Grade Pipeline Implementation #

The corrected pipeline enforces projected CRS validation, applies windowed processing with explicit overlap padding, and computes slope using Horn’s algorithm with vertical exaggeration. Hillshade generation is decoupled to allow parallel execution and memory-safe routing.

import rasterio
import numpy as np
from rasterio.windows import Window
from rasterio.enums import Resampling
import logging
import os

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

def validate_crs(src):
    """Ensure CRS is projected with linear units in meters."""
    if not src.crs.is_projected:
        raise ValueError("Input DEM must be in a projected CRS (e.g., EPSG:326xx). "
                         "Geographic CRS will produce invalid slope gradients.")
    if src.crs.linear_units.lower() not in ("metre", "meter", "m"):
        raise ValueError(f"CRS linear unit '{src.crs.linear_units}' is unsupported. "
                         "Reproject to a meter-based system before processing.")

def compute_horn_slope(dem_3x3, cell_size_m):
    """Horn's slope algorithm for 3x3 neighborhoods."""
    dzdx = (dem_3x3[0, 0] + 2*dem_3x3[1, 0] + dem_3x3[2, 0] -
            dem_3x3[0, 2] - 2*dem_3x3[1, 2] - dem_3x3[2, 2]) / (8.0 * cell_size_m)
    dzdy = (dem_3x3[2, 0] + 2*dem_3x3[2, 1] + dem_3x3[2, 2] -
            dem_3x3[0, 0] - 2*dem_3x3[0, 1] - dem_3x3[0, 2]) / (8.0 * cell_size_m)
    return np.degrees(np.arctan(np.sqrt(dzdx**2 + dzdy**2)))

def compute_hillshade(slope_deg, aspect_deg, sun_azimuth=315.0, sun_altitude=45.0):
    """Standard analytical hillshade model."""
    slope_rad = np.radians(slope_deg)
    aspect_rad = np.radians(aspect_deg)
    sun_az_rad = np.radians(sun_azimuth)
    sun_alt_rad = np.radians(sun_altitude)

    hs = 255.0 * (
        np.sin(sun_alt_rad) * np.cos(slope_rad) +
        np.cos(sun_alt_rad) * np.sin(slope_rad) * np.cos(sun_az_rad - aspect_rad)
    )
    return np.clip(hs, 0, 255).astype(np.uint8)

def process_dem_chunked(input_path, slope_out, hillshade_out, chunk_size=2048, overlap=1):
    """Memory-safe, chunked terrain analysis with boundary padding."""
    with rasterio.open(input_path) as src:
        validate_crs(src)

        meta = src.meta.copy()
        meta.update(dtype=rasterio.float32, count=1, compress="deflate", tiled=True)
        meta.update(driver="GTiff", nodata=np.nan)

        hs_meta = meta.copy()
        hs_meta.update(dtype=rasterio.uint8)

        with rasterio.open(slope_out, "w", **meta) as dst_slope, \
             rasterio.open(hillshade_out, "w", **hs_meta) as dst_hs:

            for row in range(0, src.height, chunk_size):
                for col in range(0, src.width, chunk_size):
                    # Define window with overlap for seamless stitching
                    win = Window(col, row, chunk_size, chunk_size)
                    win_padded = Window(col - overlap, row - overlap,
                                        chunk_size + 2*overlap, chunk_size + 2*overlap)

                    # Read with padding to handle boundaries
                    dem_chunk = src.read(1, window=win_padded, bound='pad', fill_value=np.nan)

                    # Compute slope and aspect (aspect requires atan2)
                    grad_y, grad_x = np.gradient(dem_chunk, src.res[1], src.res[0])
                    slope_chunk = compute_horn_slope(dem_chunk, src.res[0])
                    aspect_chunk = np.degrees(np.arctan2(grad_y, grad_x)) % 360.0

                    # Generate hillshade
                    hs_chunk = compute_hillshade(slope_chunk, aspect_chunk)

                    # Trim overlap before writing
                    slope_write = slope_chunk[overlap:-overlap, overlap:-overlap]
                    hs_write = hs_chunk[overlap:-overlap, overlap:-overlap]

                    # Align write window to original grid
                    write_win = Window(col, row, slope_write.shape[1], slope_write.shape[0])
                    dst_slope.write(slope_write, 1, window=write_win)
                    dst_hs.write(hs_write, 1, window=write_win)

                    logging.info(f"Processed window: ({row}, {col})")

if __name__ == "__main__":
    # Example execution
    process_dem_chunked("dem_projected.tif", "slope.tif", "hillshade.tif")

Spatial Debugging & Validation Protocol #

Production terrain pipelines require deterministic validation before downstream siting models consume the outputs. Implement the following checks to guarantee spatial integrity and audit readiness:

1. CRS & Resolution Verification: Always assert src.crs.is_projected and verify src.res matches expected ground sampling distance. If the DEM originates from cloud storage, validate that the internal tiling scheme aligns with your chunk size to prevent redundant I/O. Refer to PROJ documentation on coordinate transformations for authoritative CRS validation patterns.
2. Boundary Seam Inspection: After chunked execution, run a quick edge-difference check: np.diff(slope_array, axis=0) and np.diff(slope_array, axis=1). Values exceeding 5° at internal tile boundaries indicate insufficient overlap or padding misalignment. The overlap=1 parameter in the pipeline above guarantees Horn’s 3x3 kernel receives complete neighborhoods.
3. Memory Optimization & Fallback Routing: Monitor peak RSS during execution. If chunk_size triggers swap usage, implement a compliance-safe fallback: reduce chunk_size to 1024 or 512, or switch to rasterio.vrt for out-of-core processing. For CI/CD environments with strict memory caps, pre-convert DEMs to Cloud-Optimized GeoTIFF (COG) with internal tiling and overviews to minimize random disk access.
4. Audit-Ready Metadata Logging: Append processing parameters (sun azimuth, altitude, vertical exaggeration, CRS EPSG, chunk dimensions, and software versions) to the output TIFF tags using rasterio’s update_tags(). This satisfies environmental compliance audits and ensures reproducible resource assessments across development phases.

Downstream Integration & Workflow Alignment #

Validated slope and hillshade rasters serve as foundational inputs for exclusion zoning, turbine foundation stability modeling, and access road gradient constraints. When integrated into broader Solar & Wind Resource Modeling Workflows, these terrain derivatives must align temporally and spatially with wind speed interpolation grids and microclimate correction layers.

Ensure downstream modules consume the outputs with explicit nodata masking to prevent slope-driven routing algorithms from traversing void regions or water bodies. For large-scale project portfolios, decouple hillshade rendering from slope computation: slope drives engineering constraints, while hillshade primarily supports stakeholder visualization and environmental impact reporting. This separation reduces computational overhead and maintains strict compliance with data provenance standards.

By enforcing projected CRS validation, implementing overlap-aware windowing, and embedding memory-safe chunking, terrain preprocessing pipelines achieve deterministic behavior across regional extents. The resulting outputs are immediately compatible with siting optimization engines, grid interconnection studies, and environmental permitting workflows.