Open Energy Data Portals #

import asyncio
import logging
import aiohttp
from pathlib import Path
from pydantic import BaseModel, ValidationError
from typing import Dict, Any

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

class PortalMetadata(BaseModel):
    dataset_id: str
    crs: str
    bbox: list[float]
    format: str
    last_updated: str

async def fetch_portal_metadata(api_url: str, params: Dict[str, Any]) -> PortalMetadata:
    """Asynchronously query an energy data portal catalog and validate the response schema."""
    async with aiohttp.ClientSession() as session:
        async with session.get(api_url, params=params) as response:
            response.raise_for_status()
            payload = await response.json()
            try:
                return PortalMetadata(**payload["metadata"])
            except ValidationError as e:
                raise RuntimeError(f"Invalid portal schema: {e}") from e

import geopandas as gpd
import rasterio
from rasterio.transform import from_bounds
from pyproj import CRS

def harmonize_and_validate_crs(
    gdf: gpd.GeoDataFrame,
    raster_path: str,
    target_epsg: int
) -> tuple[gpd.GeoDataFrame, rasterio.io.DatasetReader]:
    """Transform vector to target CRS and verify raster compatibility."""
    target_crs = CRS.from_epsg(target_epsg)

    if not gdf.crs:
        raise ValueError("Input GeoDataFrame lacks CRS definition. Cannot harmonize.")

    gdf_transformed = gdf.to_crs(target_crs)

    # Validate raster CRS and bounds alignment
    with rasterio.open(raster_path) as src:
        if not src.crs:
            raise ValueError("Raster dataset lacks CRS definition. Rejecting.")
        if src.crs != target_crs:
            raise RuntimeError(
                f"CRS mismatch: Raster {src.crs} != Target {target_crs}. "
                "Reproject raster before pipeline execution."
            )
        return gdf_transformed, src

import numpy as np
from rasterio.windows import Window

def generate_raster_windows(src: rasterio.io.DatasetReader, chunk_size: int = 1024):
    """Yield memory-bounded raster windows for chunked processing."""
    for col in range(0, src.width, chunk_size):
        for row in range(0, src.height, chunk_size):
            width = min(chunk_size, src.width - col)
            height = min(chunk_size, src.height - row)
            yield Window(col, row, width, height)

async def process_chunk_async(
    window: Window,
    src: rasterio.io.DatasetReader,
    constraint_mask: np.ndarray
) -> Dict[str, Any]:
    """Read a raster window, apply spatial mask, and compute suitability metrics."""
    data = src.read(1, window=window)
    transform = src.window_transform(window)

    # Mask out constrained areas (e.g., wetlands, protected lands)
    valid_pixels = data[~constraint_mask]

    return {
        "window": window,
        "mean_value": float(np.nanmean(valid_pixels)) if valid_pixels.size > 0 else np.nan,
        "valid_count": int(valid_pixels.size),
        "transform": transform
    }

def run_spatial_quality_checks(gdf: gpd.GeoDataFrame, raster_src: rasterio.io.DatasetReader) -> None:
    """Execute mandatory spatial validation gates."""
    # 1. Geometry validity
    invalid = gdf[~gdf.is_valid]
    if not invalid.empty:
        logging.warning("Repairing %d invalid geometries before overlay.", len(invalid))
        gdf.geometry = gdf.geometry.make_valid()

    # 2. Extent overlap verification
    raster_bounds = raster_src.bounds
    gdf_bounds = gdf.total_bounds
    if not (
        gdf_bounds[0] <= raster_bounds[2] and
        gdf_bounds[2] >= raster_bounds[0] and
        gdf_bounds[1] <= raster_bounds[3] and
        gdf_bounds[3] >= raster_bounds[1]
    ):
        raise ValueError("Vector and raster extents do not overlap. Aborting pipeline.")

    # 3. Null/NaN band check
    if raster_src.count > 0:
        stats = raster_src.read(1, masked=True)
        if np.all(stats.mask):
            raise RuntimeError("Raster band is entirely masked/null. Check source integrity.")

    logging.info("All spatial quality gates passed.")

async def execute_suitability_pipeline(
    api_catalog_url: str,
    vector_constraints_path: str,
    raster_irradiance_path: str,
    target_epsg: int = 32611
) -> gpd.GeoDataFrame:
    """Orchestrate full open-portal ingestion, scoring, and compliance routing."""
    logging.info("Starting suitability pipeline...")

    # 1. Ingest & validate metadata
    meta = await fetch_portal_metadata(api_catalog_url, {"format": "geojson"})
    logging.info("Catalog metadata validated: %s", meta.dataset_id)

    # 2. Load & harmonize CRS
    constraints = gpd.read_file(vector_constraints_path)
    constraints, raster_src = harmonize_and_validate_crs(constraints, raster_irradiance_path, target_epsg)

    # 3. Run quality gates
    run_spatial_quality_checks(constraints, raster_src)

    # 4. Async chunked processing
    tasks = []
    for window in generate_raster_windows(raster_src):
        tasks.append(process_chunk_async(window, raster_src, constraints.geometry.values))

    results = await asyncio.gather(*tasks)

    # 5. Aggregate & route
    composite_score = np.mean([r["mean_value"] for r in results if not np.isnan(r["mean_value"])])
    logging.info("Pipeline complete. Composite irradiance score: %.2f kWh/m²/day", composite_score)

    # Attach audit trail to vector output for compliance routing
    constraints["pipeline_score"] = composite_score
    constraints["audit_timestamp"] = pd.Timestamp.utcnow().isoformat()
    constraints["portal_version"] = meta.last_updated

    return constraints