10 min read 2231 words

Automating hillshade and slope analysis for wind turbine siting

When scaling terrain preprocessing for utility-scale wind development, automated hillshade and slope extraction fails at the intersection of memory limits, coordinate reference system (CRS) unit mismatches, and windowed raster boundary discontinuities. The symptom is consistent: a pipeline runs cleanly on a small test extent but, on a regional 1 m LiDAR or 30 m SRTM mosaic covering 500+ km², either crashes with MemoryError or writes artificial 90° slope cliffs along every tile seam. This page is the slope-and-aspect derivative stage of the broader Terrain & Shadow Analysis Pipelines workflow, and it isolates the root causes, surfaces them with a pre-flight check, and delivers a chunked, CRS-validated fix that is safe to run inside a CI/CD gate before a siting model ever consumes the output.

Scenario: MemoryError and 90° slope cliffs at tile seams

Two distinct failures present from the same batch script and break the same pipeline stage — the slope/aspect derivative step that feeds turbine micro-siting and access-road routing:

  • On large mosaics, numpy.gradient over a full in-memory array raises MemoryError (or the kernel silently kills the worker) once the DEM exceeds available RAM on a cloud runner.
  • On naively tiled DEMs, slope drops to zero or spikes to a hard 90° wall along every internal tile boundary, because the derivative kernel never sees the neighbouring tile’s edge pixels.

Both must be solved together: chunking fixes the memory ceiling but introduces the seam artifact unless each window carries an overlap halo. The minimal failing pattern looks reasonable and is exactly what most teams ship first:

python
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_deg = np.degrees(np.arctan(np.sqrt(grad_x**2 + grad_y**2)))

This produces three immediate defects: MemoryError on DEMs larger than available RAM, slope cliffs at any tile boundary if later chunked without overlap, and physically wrong slope magnitudes whenever the horizontal units are degrees rather than metres.

Root-cause analysis

Three compounding causes drive the failure, each mapping to a specific stage of the fix:

  1. CRS unit mismatch. A geographic CRS such as EPSG:4326 expresses horizontal distance in degrees while elevation stays in metres, so the gradient ratios and are dimensionally invalid. Slope collapses toward zero or inflates with latitude. Computing slope on degree-spaced grids is the single most common silent error, which is why a coordinate reference system reprojection to a metric UTM zone (e.g. EPSG:32610) must happen before any derivative is taken — the same EPSG:4326 to projected-grid alignment discipline used across the resource-modeling pipeline.
  2. Windowed boundary artifacts. A 3×3 derivative kernel needs one ring of neighbouring pixels. A plain rasterio window with no padding starves the edge rows and columns, so the kernel reads off-window zeros and emits seams. The defect is invisible at small scale and only appears once tiling kicks in.
  3. Memory fragmentation. Reading the full DEM, then allocating float64 copies for numpy.gradient and scipy.ndimage intermediates, multiplies peak RSS several times over the on-disk size and exhausts constrained runners.

The corrected slope formula, applied on a metric grid with cell size , is:

Broken versus corrected slope and hillshade pipeline A warning-coloured three-stage broken pipeline ending in MemoryError and seam cliffs, above a success-coloured four-stage corrected pipeline that enforces a projected CRS, pads windows with an overlap halo, applies a Horn kernel, and trims the halo for a seamless write. BROKEN Read full DEM EPSG:4326, whole array np.gradient assumes 1:1 units Slope cliffs at seams + MemoryError CORRECTED Validate CRS projected, units=metre Window + 1px overlap halo, boundless NaN pad Horn 3x3 slope + hillshade Trim overlap seamless write

Pre-flight validation

Run a cheap check that surfaces the two structural root causes — non-metric CRS and a chunk size that will swap — before the expensive windowed pass starts. It opens only the header and metadata, so it is safe to call in a CI/CD gate.

python
import rasterio
import numpy as np

def preflight_dem(input_path: str, chunk_size: int = 2048,
                  max_window_bytes: int = 512 * 1024 * 1024) -> None:
    """Fail fast on the two structural causes of slope/hillshade failure."""
    with rasterio.open(input_path) as src:
        if not src.crs or not src.crs.is_projected:
            raise ValueError(
                f"DEM CRS {src.crs} is not projected. Reproject to a metric "
                "UTM zone (e.g. EPSG:32610) before slope computation."
            )
        units = (src.crs.linear_units or "").lower()
        if units not in ("metre", "meter", "m"):
            raise ValueError(f"CRS linear unit '{units}' is not metres.")

        # float64 working copy of one padded window is the memory hot spot
        win_px = (chunk_size + 2) ** 2
        est_bytes = win_px * np.dtype(np.float64).itemsize * 3  # dem + 2 gradients
        if est_bytes > max_window_bytes:
            raise ValueError(
                f"chunk_size={chunk_size} needs ~{est_bytes // 1_048_576} MB/window; "
                f"reduce to fit the {max_window_bytes // 1_048_576} MB budget."
            )
        print(f"OK: {src.crs} res={src.res} size={src.width}x{src.height}")

Fix implementation

The corrected pipeline enforces a projected CRS, reads each window with a one-pixel overlap halo (boundless=True pads off-extent pixels with NaN), computes slope and aspect from central differences scaled by the true ground cell size, trims the halo, and streams each tile straight to disk. Outputs are float32 slope with nodata=NaN and uint8 hillshade, both DEFLATE-compressed and internally tiled for downstream windowed reads.

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

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

def compute_slope_aspect(dem_window: np.ndarray, cell_size_m: float):
    """Central-difference slope (deg) and aspect (deg, 0-360) on a metric grid."""
    grad_y, grad_x = np.gradient(dem_window.astype(np.float64),
                                 cell_size_m, cell_size_m)
    slope_deg = np.degrees(np.arctan(np.sqrt(grad_x**2 + grad_y**2)))
    aspect_deg = np.degrees(np.arctan2(grad_y, -grad_x)) % 360.0
    return slope_deg, aspect_deg

def compute_hillshade(slope_deg, aspect_deg, sun_azimuth=315.0, sun_altitude=45.0):
    """Analytical hillshade, 0-255 uint8, for stakeholder visualization."""
    slope_rad, aspect_rad = np.radians(slope_deg), np.radians(aspect_deg)
    az_rad, alt_rad = np.radians(sun_azimuth), np.radians(sun_altitude)
    hs = 255.0 * (
        np.sin(alt_rad) * np.cos(slope_rad) +
        np.cos(alt_rad) * np.sin(slope_rad) * np.cos(az_rad - aspect_rad)
    )
    return np.clip(np.nan_to_num(hs), 0, 255).astype(np.uint8)

def process_dem_chunked(input_path, slope_out, hillshade_out,
                        chunk_size=2048, overlap=1):
    """Memory-safe, seam-free slope + hillshade over a regional DEM."""
    with rasterio.open(input_path) as src:
        cell_size_m = src.res[0]
        slope_meta = src.meta.copy()
        slope_meta.update(driver="GTiff", dtype="float32", count=1,
                          nodata=np.nan, compress="deflate", tiled=True)
        hs_meta = slope_meta.copy()
        hs_meta.update(dtype="uint8", nodata=0)

        with rasterio.open(slope_out, "w", **slope_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):
                    padded = Window(col - overlap, row - overlap,
                                    chunk_size + 2 * overlap,
                                    chunk_size + 2 * overlap)
                    dem_chunk = src.read(1, window=padded, boundless=True,
                                         fill_value=np.nan)
                    slope, aspect = compute_slope_aspect(dem_chunk, cell_size_m)
                    hs = compute_hillshade(slope, aspect)

                    slope_w = slope[overlap:-overlap, overlap:-overlap]
                    hs_w = hs[overlap:-overlap, overlap:-overlap]
                    write_win = Window(col, row, slope_w.shape[1], slope_w.shape[0])
                    dst_slope.write(slope_w.astype("float32"), 1, window=write_win)
                    dst_hs.write(hs_w, 1, window=write_win)
            dst_slope.update_tags(sun_azimuth=315.0, sun_altitude=45.0,
                                  cell_size_m=cell_size_m, source_crs=str(src.crs))
            logging.info("slope + hillshade written for %s", input_path)
Overlap halo lets the Horn kernel cross tile seams without cliffs A padded read window with a warning-coloured one-pixel halo ring around a success-coloured write window; a three-by-three kernel at the corner reads into the halo. Trimming the halo writes the inner block to the parent grid, where neighbouring tiles meet at a seam with no slope cliff. Padded read window (chunk + halo) write window (trimmed, 7x7) halo 3x3 Horn kernel centre cell reads into halo 1px overlap halo boundless NaN pad — supplies the neighbour pixels the kernel needs trim halo, write aligned Parent grid (seamless) tile A tile B seam stays continuous no zero gap, no 90° cliff

Fallback routing & performance tuning

When the default chunk_size still pressures a constrained runner, step down through these strategies rather than loading the full array:

  • Shrink the window. Drop chunk_size to 1024 or 512; peak RSS scales with the window area, so halving the side quarters the working set while the one-pixel halo keeps seams gone.
  • Cap GDAL’s block cache. Export GDAL_CACHEMAX=256 (MB) so rasterio reads do not balloon the resident set on top of your NumPy arrays.
  • Pre-build a Cloud-Optimized GeoTIFF. Convert the source DEM to a COG with internal tiling and overviews so each boundless window is a cheap block read instead of a random-access scan of a striped TIFF.
  • Use a VRT for multi-tile mosaics. Wrap many adjacent DEM tiles in a gdalbuildvrt virtual raster and stream windows out-of-core, avoiding a physical merge into one oversized file.
  • Choose resampling deliberately if you downsample. Use average (not bilinear) when reducing resolution before slope, so a single anomalous cell does not smear a false gradient across the kernel.

Downstream validation

Before slope and hillshade reach exclusion-zoning or grid capacity buffer analysis routing, assert structural integrity. This audit returns a pass/fail dict suitable for a CI/CD gate and catches the exact regressions this fix targets: wrong dtype, a non-projected output CRS, nodata bleed, and any residual seam.

python
import rasterio
import numpy as np

def audit_slope_raster(slope_path: str, max_internal_jump_deg: float = 5.0) -> dict:
    """CI/CD gate: band count, dtype, projected CRS, nodata bleed, seam check."""
    report = {}
    with rasterio.open(slope_path) as src:
        report["band_count_ok"] = src.count == 1
        report["dtype_ok"] = src.dtypes[0] == "float32"
        report["crs_projected"] = bool(src.crs and src.crs.is_projected)

        slope = src.read(1, masked=True)
        report["nodata_is_nan"] = src.nodata is None or np.isnan(src.nodata)
        report["value_range_ok"] = bool(slope.min() >= 0 and slope.max() <= 90)

        # Seam check: large abrupt jumps along interior rows/cols flag bad overlap
        d_row = np.abs(np.diff(slope.filled(np.nan), axis=0))
        d_col = np.abs(np.diff(slope.filled(np.nan), axis=1))
        worst = np.nanmax([np.nanmax(d_row), np.nanmax(d_col)])
        report["max_internal_jump_deg"] = float(worst)
        report["seam_free"] = bool(worst <= max_internal_jump_deg)

    report["passed"] = all(v for k, v in report.items()
                           if isinstance(v, bool))
    return report

A failing seam_free points back to insufficient overlap; a failing crs_projected means the reproject in coordinate reference systems for energy projects was skipped upstream. Log the returned dict alongside the update_tags provenance (sun azimuth, altitude, cell size, source CRS) so the terrain-loss figure stays defensible in permitting and project-finance review. With projected-CRS enforcement, overlap-aware windowing, and the audit gate in place, the slope and hillshade surfaces align pixel-for-pixel with the irradiance and wind grids they constrain and are ready for turbine micro-siting and interconnection studies.

Related articles

Back to Terrain & Shadow Analysis Pipelines