Add a Solar calculator with offline NASA POWER data

Estimates a solar array's typical-year production from a location and the
panels; everything else has defaults under an optional "Fine-tune" section.

The model works from 12 monthly averages: Klein's average day per month, the
Erbs monthly diffuse fraction, Collares-Pereira & Rabl / Liu & Jordan hourly
split, HDKR transposition for fixed, single-axis and dual-axis arrays, then
PVWatts' ASHRAE incidence loss, SAPM cell temperature, DC power, system
losses and optional inverter.  "Best for the year" and "Best for winter" sweep
the tilt.

Sunlight, diffuse, albedo and temperature come from NASA POWER's 2001-2020
climatology on a 1 degree grid, packed into app/public/data (3.1 MB, ~1 MB
gzipped) by scripts/build_solar_climatology.py and fetched only when the
calculator opens.  Users can replace it with their own monthly values.

Every setting lives in the URL, with defaults left out, so the Share button
shares the whole configuration.

Against PVWatts v8 for a 1 kW fixed array: Phoenix +0.6%, Seattle +8.8%,
Anchorage +8.0% annually; these are recorded as tests.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

scripts/build_solar_climatology.py

@@ -0,0 +1,170 @@
+#!/usr/bin/env python3
+"""
+Build the offline solar climatology used by the Solar calculator.
+
+Downloads NASA POWER's 2001-2020 monthly climatology for the whole globe and packs it into
+app/public/data/solar-climatology.bin, which the React app fetches when the Solar calculator opens.
+This needs the Internet and is run by a developer; WROLPi never runs it.  The output only changes
+when NASA publishes a new climatology period, so regenerating it should be rare.
+
+    python3 scripts/build_solar_climatology.py
+
+Downloaded tiles are cached (see --cache), so an interrupted run resumes where it stopped.
+
+NASA POWER data is licensed CC BY 4.0; the calculator credits it.  https://power.larc.nasa.gov/
+
+File format (all integers little-endian), decoded by app/src/components/calculators/solarData.js:
+
+    bytes 0-3   b"WSOL"
+    byte  4     format version (1)
+    byte  5     number of parameters (4)
+    bytes 6-7   reserved (0)
+    then        uint8 values, [parameter][month][row][column]
+
+Rows are 1° latitude bands centered on -89.5 ... 89.5 (south to north), columns are 1° longitude
+bands centered on -179.5 ... 179.5 (west to east).  255 means no data.  Each parameter is stored as
+round((value + offset) / step); see PARAMETERS.
+"""
+import argparse
+import json
+import math
+import pathlib
+import struct
+import sys
+import time
+import urllib.error
+import urllib.parse
+import urllib.request
+from concurrent.futures import ThreadPoolExecutor, as_completed
+
+REPO = pathlib.Path(__file__).resolve().parent.parent
+DEFAULT_OUTPUT = REPO / 'app/public/data/solar-climatology.bin'
+DEFAULT_CACHE = REPO / 'scripts/.solar_climatology_cache'
+
+API = 'https://power.larc.nasa.gov/api/temporal/climatology/regional'
+# The regional endpoint accepts at most a 10° x 10° box and one parameter per request.
+TILE = 10
+MONTHS = ['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC']
+ROWS, COLUMNS = 180, 360
+MISSING = 255
+VERSION = 1
+
+# (NASA parameter, offset, step).  Must match PARAMETERS in solarData.js, in the same order.
+PARAMETERS = [
+    ('ALLSKY_SFC_SW_DWN', 0, 0.04),  # global horizontal, kWh/m²/day, 0-10.16
+    ('ALLSKY_SFC_SW_DIFF', 0, 0.04),  # diffuse horizontal, kWh/m²/day
+    ('ALLSKY_SRF_ALB', 0, 0.004),  # surface albedo, 0-1.016
+    ('T2M', 70, 0.5),  # air temperature at 2 m, °C, -70 to +57
+]
+
+
+def tile_url(parameter, lat_min, lon_min):
+    query = urllib.parse.urlencode({
+        'parameters': parameter,
+        'community': 'RE',
+        'latitude-min': lat_min,
+        'latitude-max': lat_min + TILE,
+        'longitude-min': lon_min,
+        'longitude-max': lon_min + TILE,
+        'format': 'JSON',
+    })
+    return f'{API}?{query}'
+
+
+def fetch_tile(parameter, lat_min, lon_min, cache, pause, retries=8):
+    """Return the tile's GeoJSON, from the cache when possible.
+
+    NASA rate limits the API (HTTP 429); back off for up to 10 minutes rather than give up, since a
+    full build is about 2,600 requests.
+    """
+    path = cache / f'{parameter}_{lat_min}_{lon_min}.json'
+    if path.is_file():
+        return json.loads(path.read_text())
+
+    delay = 30
+    for attempt in range(retries):
+        try:
+            time.sleep(pause)
+            with urllib.request.urlopen(tile_url(parameter, lat_min, lon_min), timeout=180) as response:
+                data = json.loads(response.read())
+            if 'features' not in data:
+                raise ValueError(f'No features: {data.get("messages")}')
+            path.write_text(json.dumps(data))
+            return data
+        except (urllib.error.URLError, TimeoutError, ValueError) as e:
+            if attempt == retries - 1:
+                raise
+            print(f'  retrying {parameter} {lat_min},{lon_min} in {delay}s: {e}', file=sys.stderr, flush=True)
+            time.sleep(delay)
+            delay = min(delay * 2, 600)
+
+
+def bin_points(points):
+    """Average points into 1° cells.  `points` maps (lon, lat) to 12 monthly values.
+
+    Solar parameters are already on the 1° grid (one point per cell).  Temperature comes on a finer
+    0.5° x 0.625° grid, so several points average into each cell.
+    """
+    sums = {}
+    for (lon, lat), values in points.items():
+        row = min(ROWS - 1, max(0, math.floor(lat + 90)))
+        column = math.floor(lon + 180) % COLUMNS
+        cell = sums.setdefault((row, column), [[0.0, 0] for _ in MONTHS])
+        for month, value in enumerate(values):
+            if value is not None:
+                cell[month][0] += value
+                cell[month][1] += 1
+    return {cell: [s / n if n else None for s, n in months] for cell, months in sums.items()}
+
+
+def quantize(value, offset, step):
+    if value is None:
+        return MISSING
+    return max(0, min(MISSING - 1, round((value + offset) / step)))
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    parser.add_argument('--output', type=pathlib.Path, default=DEFAULT_OUTPUT)
+    parser.add_argument('--cache', type=pathlib.Path, default=DEFAULT_CACHE)
+    parser.add_argument('--workers', type=int, default=1, help='Concurrent requests; NASA rate limits the API.')
+    parser.add_argument('--pause', type=float, default=1.0, help='Seconds to wait before each request.')
+    args = parser.parse_args()
+    args.cache.mkdir(parents=True, exist_ok=True)
+
+    tiles = [(lat, lon) for lat in range(-90, 90, TILE) for lon in range(-180, 180, TILE)]
+    out = bytearray(b'WSOL' + struct.pack('<BBH', VERSION, len(PARAMETERS), 0))
+
+    for parameter, offset, step in PARAMETERS:
+        print(f'{parameter}: {len(tiles)} tiles', flush=True)
+        points = {}
+        with ThreadPoolExecutor(max_workers=args.workers) as pool:
+            futures = [pool.submit(fetch_tile, parameter, lat, lon, args.cache, args.pause) for lat, lon in tiles]
+            for done, future in enumerate(as_completed(futures), 1):
+                data = future.result()
+                header = data.get('header')
+                fill = header.get('fill_value', -999) if isinstance(header, dict) else -999
+                for feature in data['features']:
+                    lon, lat = feature['geometry']['coordinates'][:2]
+                    monthly = feature['properties']['parameter'][parameter]
+                    # Tiles share edges; a point on an edge arrives twice with the same value.
+                    points[(lon, lat)] = [None if monthly[m] == fill else monthly[m] for m in MONTHS]
+                if done % 50 == 0:
+                    print(f'  {done}/{len(tiles)}', flush=True)
+
+        cells = bin_points(points)
+        grid = bytearray([MISSING]) * (12 * ROWS * COLUMNS)
+        for (row, column), months in cells.items():
+            for month, value in enumerate(months):
+                grid[(month * ROWS + row) * COLUMNS + column] = quantize(value, offset, step)
+        missing = grid.count(MISSING)
+        print(f'  {len(points)} points into {len(cells)} cells; {missing} missing values')
+        out += grid
+
+    args.output.parent.mkdir(parents=True, exist_ok=True)
+    args.output.write_bytes(out)
+    print(f'Wrote {args.output} ({len(out):,} bytes)')
+
+
+if __name__ == '__main__':
+    main()