Solar Irradiance Raster Processing #

Accurate solar resource assessment hinges on rigorous raster ingestion, spatial normalization, and radiometric resampling long before any yield model or financial pro forma runs. This workflow is the data-preparation stage of the broader Solar & Wind Resource Modeling Workflows pipeline: raw satellite-derived and reanalysis products — NSRDB, PVGIS, NASA POWER, CAMS — almost never arrive on a uniform coordinate reference system, resolution, or temporal cadence, and feeding them straight into a model silently corrupts the capacity factor a lender treats as ground truth. The specific failure mode this stage exists to eliminate is CRS drift in multi-source raster stacks: when two irradiance surfaces disagree on projection, pixel registration, or grid origin by even a fraction of a cell, every downstream spatial join with parcel boundaries, terrain masks, or transmission corridors inherits a systematic radiometric bias that no later step can detect or repair.

The goal of the stage is deterministic: turn heterogeneous Global Horizontal Irradiance (GHI), Direct Normal Irradiance (DNI), and Diffuse Horizontal Irradiance (DHI) inputs into a single analysis-ready grid that preserves radiometric integrity, carries explicit provenance, and aligns pixel-for-pixel with every other layer in the assessment. This page covers the conceptual foundation, the prerequisites, a full runnable processing function, the three failure modes that break naive pipelines, the scalability patterns for continental archives, and the audit trail that makes the output defensible in permitting and project-finance review.

Why naive raster stacking fails #

The intuitive approach — load every GeoTIFF, call rasterio.merge or stack the arrays, and average — fails because raster algebra has no tolerance for spatial disagreement. Three independent mismatches compound into the bias this stage must prevent.

First, projection mismatch. Irradiance rasters are typically distributed in geographic coordinates (EPSG:4326) with degree-based pixel spacing, while project development requires a projected metric system — a UTM zone such as EPSG:32610 or a state plane CRS — for accurate area, buffer, and distance work. Operating in degrees stretches north–south distance against east–west distance by a latitude-dependent factor, so a “1 km” buffer around a substation interconnection point computed on an unprojected irradiance grid is not actually 1 km. Enforcing coordinate reference system alignment into a single projected target is the precondition for everything else.

Second, pixel registration and grid origin divergence. NASA POWER ships a 0.5° grid with center-registered pixels; PVGIS and CAMS often deliver edge-registered cells or finer 0.01° grids. Mixing registration conventions shifts pixel centers by half a cell, and a half-pixel shift on a continental GHI surface translates into tens of kilometres of misplaced irradiance. The affine transforms must be reconciled to a common origin before any resampling.

Third, radiometrically wrong resampling. Resampling is not a single operation — the correct kernel depends on what the band represents. Nearest-neighbor preserves categorical masks (cloud flags, land/water) but destroys continuous fields; bilinear preserves continuous irradiance without inventing new extremes; cubic convolution sharpens high-frequency temporal derivatives at the cost of overshoot near coastlines and terrain edges. Applying nearest-neighbor to a GHI surface, or bilinear to a quality flag, silently corrupts the data while producing a file that looks valid.

A clean pipeline therefore decouples ingestion, transformation, and validation into discrete, testable stages so that drift, registration mismatch, or resolution disagreement is caught and rejected before resource aggregation rather than discovered after the financial model has already consumed the bias.

Prerequisites and data requirements #

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

• Library versions: rasterio>=1.3, numpy>=1.24, pyproj>=3.5. GDAL underpins rasterio; keep it >=3.6 so calculate_default_transform honours resolution correctly. Authoritative 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 chosen for the project’s UTM zone (for example EPSG:32610 for the US West Coast). Always store the EPSG integer, never an unqualified “UTM 10N” string.
• Input geometry: single-band or multi-band GeoTIFF / Cloud-Optimized GeoTIFF surfaces of GHI, DNI, or DHI in W/m² or kWh/m²/day, each carrying a defined CRS and nodata value. Surfaces missing a CRS must be rejected at ingest, not assumed — see spatial data quality validation for the cleaning patterns this stage assumes upstream.
• Source provenance: acquisition date range, source portal, and product version for each input, sourced from one of the documented open energy data portals so the output’s lineage is auditable.
• Sanity bounds: daily-mean GHI realistically falls between 0 and roughly 12 kWh/m²/day depending on latitude and atmospheric clarity. A useful normalization is the dimensionless clearness index, the ratio of surface to extraterrestrial irradiance:

where ${\mathrm{GHI}}_0$ is top-of-atmosphere horizontal irradiance. Any pixel with $K_t>1$ after processing is physically impossible and signals a unit, scaling, or resampling defect.

Core implementation #

The function below ingests a list of GHI rasters, validates each one, reprojects and resamples it into the target metric CRS with windowed I/O, and returns per-file QA statistics. It uses bilinear resampling to preserve radiometric continuity, writes tiled LZW-compressed float32 GeoTIFFs with nodata=NaN, and orchestrates independent files concurrently with asyncio so I/O latency overlaps across the portfolio. Variable names are energy-specific throughout.

import asyncio
import logging
from pathlib import Path
from typing import Dict, List

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.crs import CRS
from concurrent.futures import ThreadPoolExecutor

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


def validate_spatial_metadata(src: rasterio.DatasetReader,
                              target_epsg: int,
                              target_res_m: float) -> None:
    """Reject irradiance rasters that violate alignment preconditions."""
    if src.crs is None or not src.crs.is_valid:
        raise ValueError(f"{src.name}: source CRS is undefined or invalid.")
    if src.bounds.left >= src.bounds.right or src.bounds.bottom >= src.bounds.top:
        raise ValueError(f"{src.name}: degenerate raster bounds detected.")
    if src.nodata is None:
        logging.warning("%s: no nodata declared; NaN fill will be assumed.", src.name)
    # Flag a >5% native-vs-target resolution gap so resampling is intentional.
    native_res_m = abs(src.res[0]) * (111_320 if src.crs.is_geographic else 1)
    if abs(native_res_m - target_res_m) / target_res_m > 0.05:
        logging.info("%s: native ~%.0f m differs >5%% from target %.0f m; resampling enforced.",
                     src.name, native_res_m, target_res_m)


def reproject_irradiance(src_path: Path,
                         dst_path: Path,
                         target_epsg: int,
                         target_res_m: float,
                         chunk_size: int = 2048) -> Dict[str, float]:
    """Reproject one GHI raster to a metric grid with windowed, radiometric resampling."""
    target_crs = CRS.from_epsg(target_epsg)
    with rasterio.open(src_path) as src:
        validate_spatial_metadata(src, target_epsg, target_res_m)

        dst_transform, dst_width, dst_height = calculate_default_transform(
            src.crs, target_crs, src.width, src.height, *src.bounds,
            resolution=target_res_m,
        )
        profile = src.profile | {
            "driver": "GTiff", "crs": target_crs, "transform": dst_transform,
            "width": dst_width, "height": dst_height, "dtype": "float32",
            "nodata": np.nan, "tiled": True, "blockxsize": chunk_size,
            "blockysize": chunk_size, "compress": "lzw",
        }

        with rasterio.open(dst_path, "w", **profile) as dst:
            for band in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, band),
                    destination=rasterio.band(dst, band),
                    src_crs=src.crs,
                    dst_crs=target_crs,
                    dst_transform=dst_transform,
                    # Bilinear preserves continuous irradiance without inventing extremes.
                    resampling=Resampling.bilinear,
                    num_threads=4,
                )

    with rasterio.open(dst_path) as out:
        ghi_array = out.read(1, masked=True)
        return {
            "valid_px": int(ghi_array.count()),
            "mean_ghi_kwh_m2": float(np.ma.mean(ghi_array)),
            "max_ghi_kwh_m2": float(np.ma.max(ghi_array)),
            "crs_aligned": out.crs.to_epsg() == target_epsg,
        }


async def run_irradiance_pipeline(src_paths: List[Path],
                                  dst_dir: Path,
                                  target_epsg: int = 32610,
                                  target_res_m: float = 1000.0,
                                  max_concurrency: int = 3) -> Dict[str, Dict[str, float]]:
    """Orchestrate async, memory-bounded irradiance processing across a portfolio."""
    dst_dir.mkdir(parents=True, exist_ok=True)
    loop = asyncio.get_running_loop()
    semaphore = asyncio.Semaphore(max_concurrency)  # cap concurrent disk I/O
    results: Dict[str, Dict[str, float]] = {}

    with ThreadPoolExecutor(max_workers=max_concurrency) as executor:
        async def _process(src_path: Path) -> None:
            dst_path = dst_dir / f"aligned_{src_path.name}"
            async with semaphore:
                logging.info("Processing %s", src_path.name)
                stats = await loop.run_in_executor(
                    executor, reproject_irradiance,
                    src_path, dst_path, target_epsg, target_res_m,
                )
            results[src_path.name] = stats

        await asyncio.gather(*(_process(p) for p in src_paths))
    return results

The reprojection itself is delegated to rasterio.warp.reproject, which streams the warp through GDAL’s windowed engine rather than materializing the full source array — the single most important detail for keeping continental archives inside a workstation’s RAM budget. Running each file in a thread-pool executor behind a semaphore lets independent surfaces overlap their I/O without thrashing the disk.

Error handling and edge cases #

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

Undefined or geographic CRS reaching the warp. A surface with no CRS, or one left in EPSG:4326 when the target is metric, must be stopped at validation. validate_spatial_metadata already raises on a missing CRS; reject silent geographic inputs before they contaminate a metric stack:

with rasterio.open(src_path) as src:
    if src.crs is None:
        raise ValueError(f"{src_path.name}: refusing to process — no CRS declared.")
    if src.crs.is_geographic and target_epsg not in (4326,):
        logging.warning("%s: geographic source (EPSG:%s) reprojecting to EPSG:%s.",
                        src_path.name, src.crs.to_epsg(), target_epsg)

Disjoint or non-overlapping extents. Two surfaces that do not share a footprint produce an empty or all-nodata result that downstream code happily averages to NaN. Test for overlap before merging, the same ValueError: Input shapes do not overlap raster condition handled in depth in Stacking NASA POWER and PVGIS rasters in rasterio:

from rasterio.coords import disjoint_bounds

if disjoint_bounds(src_a.bounds, src_b.bounds):
    raise ValueError("Irradiance surfaces have non-overlapping extents; check source tiling.")

Implicit reprojection memory spike. Calling src.read() on a full-resolution continental array before warping is the classic MemoryError trigger on a 32 GB workstation. The core function avoids it by handing band references to reproject so GDAL streams windows; never read the whole array eagerly when a windowed path exists. For surfaces large enough that even the warp output strains memory, cap GDAL’s cache and process in blocks rather than raising the cache ceiling.

Performance and scalability #

Scaling from a single feasibility site to a regional portfolio is a question of bounding memory and saturating I/O, not buying more RAM:

• Windowed reads and tiled writes. Tiled output (blockxsize/blockysize of 512–2048) lets every downstream consumer read the same windows the warp wrote, keeping peak memory proportional to one tile rather than the whole grid.
• Async over files, threads within a file. Coarse concurrency belongs at the file level via asyncio and a semaphore; fine-grained parallelism belongs inside the warp via num_threads. Nesting them lets a portfolio run overlap disk latency while each file still uses every core.
• GDAL cache tuning. Set GDAL_CACHEMAX (e.g. 512) to bound block cache; a runaway cache, not the data, is usually what exhausts memory during batch runs.
• Match the resampling kernel to the band. Bilinear for continuous GHI/DNI, nearest for categorical masks, average when downsampling to a coarser monthly grid — choosing average for downsampling avoids the aliasing that bilinear introduces when the target cell spans many source cells.
• Pre-flight the resolution gap. Logging a >5% native-vs-target gap before the run surfaces silent over- or under-sampling that would otherwise only appear as a quiet bias in the aggregated capacity factor. Once aligned, these surfaces feed directly into temporal data aggregation for monthly and seasonal reduction.

Validation and audit trail #

A processed raster 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: extent and pixel-alignment verification, nodata consistency, and the physical sanity bounds from the prerequisites.

def assert_irradiance_integrity(dst_path: Path, target_epsg: int) -> None:
    """CI/CD gate: fail the build if a processed GHI surface is non-compliant."""
    with rasterio.open(dst_path) as out:
        assert out.crs.to_epsg() == target_epsg, "CRS not aligned to target."
        assert out.dtypes[0] == "float32", "Unexpected dtype; expected float32."
        ghi = out.read(1, masked=True)
        assert ghi.count() > 0, "No valid pixels — possible disjoint extent."
        assert float(np.ma.max(ghi)) <= 13.0, "GHI exceeds physical ceiling (kWh/m²/day)."
        # Embed ISO 19115-style provenance directly in the GeoTIFF header.
    with rasterio.open(dst_path, "r+") as out:
        out.update_tags(
            SOURCE="NSRDB/PVGIS",
            PROCESSING="bilinear-reproject",
            CRS_EPSG=str(target_epsg),
            QA_STATUS="passed",
        )

Pixel alignment is the non-negotiable invariant: when surfaces are later stacked or intersected with transmission line and substation mapping layers or with a grid capacity buffer analysis, sub-pixel shifts compound into siting errors. Enforce an explicit tolerance — for example ±0.5 m for UTM-projected assets — and log transformation residuals so the output is auditable for regulatory submission and project-finance due diligence. Embedding acquisition timestamps and source version 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 safe inputs for cross-validation workflows such as Stacking NASA POWER and PVGIS rasters in rasterio and for uncertainty quantification across providers.

Related #

• Solar & Wind Resource Modeling Workflows — the parent pipeline this stage feeds.
• Stacking NASA POWER and PVGIS rasters in rasterio — resolving overlap and registration errors when merging multi-source surfaces.
• Temporal Data Aggregation — reducing aligned irradiance stacks to monthly and seasonal statistics.
• Terrain & Shadow Analysis Pipelines — horizon masking that requires identical pixel alignment to the irradiance grid.
• Coordinate Reference Systems for Energy Projects — the projection foundations this stage enforces.
• Open Energy Data Portals — sourcing and provenance for NSRDB, PVGIS, and NASA POWER inputs.