diff --git a/.gitignore b/.gitignore index 655c907f..a3ab9bc5 100644 --- a/.gitignore +++ b/.gitignore @@ -150,3 +150,6 @@ app/cypress/videos release/latest.json perf_results/ scripts/api_perf_hosts.json + +# Tile cache for scripts/build_solar_climatology.py +/scripts/.solar_climatology_cache/ diff --git a/app/public/data/solar-climatology.bin b/app/public/data/solar-climatology.bin new file mode 100644 index 00000000..9a5bcb68 Binary files /dev/null and b/app/public/data/solar-climatology.bin differ diff --git a/app/src/components/Calculators.js b/app/src/components/Calculators.js index 99878cde..da82243f 100644 --- a/app/src/components/Calculators.js +++ b/app/src/components/Calculators.js @@ -14,6 +14,7 @@ import {FoodStorageCalculator} from "./calculators/FoodStorageCalculator"; import {OneTimePadCalculator} from "./calculators/OneTimePadCalculator"; import {VinDecoderCalculator} from "./calculators/VinDecoderCalculator"; import {SpeedTestCalculator} from "./calculators/SpeedTestCalculator"; +import {SolarCalculator} from "./calculators/SolarCalculator"; // Calculators are shown in these titled groups, in this order. Each calculator carries a // `group` matching one of these names; any calculator whose group is missing/unknown falls @@ -35,6 +36,7 @@ export const useCalculators = () => { {key: 'ratio', icon: 'th large', button: 'Ratio', group: 'Engineering', contents: }, {key: 'electrical', icon: 'lightning', button: 'Electrical', group: 'Engineering', contents: }, {key: 'drive', icon: 'cogs', button: 'Drive Ratio', group: 'Engineering', contents: }, + {key: 'solar', icon: 'sun', button: 'Solar', group: 'Engineering', contents: }, {key: 'water', icon: 'tint', button: 'Water Storage', group: 'Preparedness & Storage', contents: }, {key: 'ration', icon: 'food', button: 'Food Storage', group: 'Preparedness & Storage', contents: }, {key: 'antenna', icon: 'signal', button: 'Antenna', group: 'Radio & Communications', contents: }, diff --git a/app/src/components/calculators/SolarCalculator.js b/app/src/components/calculators/SolarCalculator.js new file mode 100644 index 00000000..dfc966c7 --- /dev/null +++ b/app/src/components/calculators/SolarCalculator.js @@ -0,0 +1,541 @@ +import React from "react"; +import {useSearchParams} from "react-router"; +import { + Accordion, + Button, + Grid, + Group, + Header, + Loader, + Message, + NumberInput, + Panel, + Select, + Statistic, + StatisticGroup, + Table, + TextInput, + Toggle, +} from "../ui"; +import {InfoPopup, roundDigits} from "../Common"; +import { + DEFAULT_LOSSES, + defaultAzimuth, + estimateProduction, + formatMonthly, + INVERTER_NOMINAL_EFFICIENCY, + LOSS_CATEGORIES, + MODULE_TYPES, + MONTH_NAMES, + MOUNTINGS, + optimalTilt, + parseCoordinate, + parseMonthly, + totalLoss, + TRACKING, +} from "./solar"; +import {CLIMATOLOGY_SOURCE, lookupClimatology, useClimatology} from "./solarData"; + +// Every input lives in the URL query so the Share button (which shares window.location.href) shares +// the whole configuration. Values equal to their default are left out to keep the URL, and its QR +// code, small. Keys are short and must stay stable; shared links depend on them. +const DEFAULT_PANEL_WATTS = 400; +const DEFAULT_PANEL_COUNT = 1; + +const COMPASS = [ + {value: '0', label: 'N'}, {value: '45', label: 'NE'}, {value: '90', label: 'E'}, {value: '135', label: 'SE'}, + {value: '180', label: 'S'}, {value: '225', label: 'SW'}, {value: '270', label: 'W'}, {value: '315', label: 'NW'}, +]; + +const fToC = f => (f - 32) * 5 / 9; +const cToF = c => c * 9 / 5 + 32; + +const numberOr = (value, fallback) => { + if (value === null || value === undefined || value === '') { + return fallback; + } + const n = Number(value); + return Number.isFinite(n) ? n : fallback; +}; + +const fmt = (value, digits = 2) => (Number.isFinite(value) ? `${roundDigits(value, digits)}` : '—'); + +// Read and write the calculator's query parameters. Writes replace the history entry, so typing does +// not fill the Back button with one entry per keystroke. +export function useSolarParams() { + const [searchParams, setSearchParams] = useSearchParams(); + + const get = key => searchParams.get(key); + + // `updates` maps keys to values; null, undefined and '' remove the key. + const set = updates => setSearchParams(prev => { + const next = new URLSearchParams(prev); + Object.entries(updates).forEach(([key, value]) => { + if (value === null || value === undefined || value === '') { + next.delete(key); + } else { + next.set(key, `${value}`); + } + }); + return next; + }, {replace: true}); + + // Store a number, or remove it when it equals the default. + const setNumber = (key, value, fallback) => + set({[key]: value === '' || Number(value) === fallback ? null : value}); + + return {searchParams, get, set, setNumber}; +} + +// Turn the query parameters into the model's inputs. Invalid values fall back to defaults, so a +// hand-edited or truncated link never crashes the page. +export function systemFromParams(get) { + const latitude = parseCoordinate(get('lat'), 'lat'); + const longitude = parseCoordinate(get('lon'), 'lon'); + const fahrenheit = get('tu') === 'f'; + + // Own sunlight data replaces the NASA data when present. + const ownData = get('ghi') !== null; + const ghi = parseMonthly(get('ghi')).map(v => (v !== null && v >= 0 ? v : null)); + const temperatureInput = parseMonthly(get('temp')); + const temperature = temperatureInput.map(v => (v === null ? null : (fahrenheit ? fToC(v) : v))); + + const panelWatts = Math.max(0, numberOr(get('pw'), DEFAULT_PANEL_WATTS)); + const panelCount = Math.max(0, Math.round(numberOr(get('pn'), DEFAULT_PANEL_COUNT))); + const kwp = panelWatts * panelCount / 1000; + + const tracking = TRACKING[get('track')] ? get('track') : 'fixed'; + const mounting = MOUNTINGS[get('mount')] ? get('mount') : 'open'; + const gamma = numberOr(get('gamma'), MODULE_TYPES.standard.gamma); + const defaultTilt = latitude === null ? 0 : Math.round(Math.abs(latitude)); + const tilt = Math.min(90, Math.max(0, numberOr(get('tilt'), defaultTilt))); + const azimuth = numberOr(get('az'), defaultAzimuth(latitude ?? 0)); + // Null means "use the monthly albedo from the data", which includes seasonal snow. + const albedoInput = numberOr(get('alb'), null); + const albedo = albedoInput === null ? null : Math.min(1, Math.max(0, albedoInput)); + + const lossValues = parseMonthly(get('losses')); + const losses = LOSS_CATEGORIES.map((c, i) => { + const v = lossValues[i]; + return v !== null && v >= 0 && v <= 100 ? v : c.default; + }); + + const inverterOn = get('inv') === '1'; + const acRating = numberOr(get('iac'), roundDigits(kwp / 1.2, 2)); + const efficiency = numberOr(get('ieff'), INVERTER_NOMINAL_EFFICIENCY); + + return { + latitude, longitude, fahrenheit, ownData, ghi, temperatureInput, temperature, panelWatts, panelCount, + kwp, tracking, mounting, gamma, defaultTilt, tilt, azimuth, albedo, losses, inverterOn, acRating, + efficiency, + }; +} + +// The sunlight inputs for the model: the user's own values, or the NASA data at the location. +export function sunlightFor(s, grid) { + if (s.ownData) { + return {ghi: s.ghi, diffuse: null, temperature: s.temperature, albedo: null}; + } + if (!grid || s.latitude === null || s.longitude === null) { + return null; + } + return lookupClimatology(grid, s.latitude, s.longitude); +} + +// What stops an estimate, or an empty list when it can run. +export function missingInputs(s, sunlight) { + const missing = []; + if (s.latitude === null) { + missing.push('a latitude'); + } + if (!s.ownData && s.longitude === null) { + missing.push('a longitude'); + } + if (s.ownData && s.ghi.some(v => v === null)) { + missing.push('sunlight for all 12 months'); + } + if (!(s.kwp > 0)) { + missing.push('panel watts and count'); + } + if (missing.length === 0 && (!sunlight || sunlight.ghi.some(v => v === null))) { + missing.push('a location with sunlight data'); + } + return missing; +} + +const modelInputs = (s, sunlight) => ({ + latitude: s.latitude, + ghi: sunlight.ghi, + diffuse: sunlight.diffuse, + temperature: sunlight.temperature, + albedo: s.albedo ?? sunlight.albedo ?? undefined, + kwp: s.kwp, + tilt: s.tilt, + azimuth: s.azimuth, + tracking: s.tracking, + mounting: s.mounting, + gamma: s.gamma, + losses: s.losses, + inverter: s.inverterOn ? {acRating: s.acRating, efficiency: s.efficiency} : null, +}); + +const compassLabel = azimuth => { + const exact = COMPASS.find(c => Number(c.value) === azimuth); + return exact ? exact.label : `${azimuth}°`; +}; + +// Small explanatory text at the top of a collapsed section. +const Hint = ({children}) =>

{children}

; + +const ESTIMATE_INFO = 'A typical-year estimate from monthly averages. Real years vary; expect roughly ±10–15% ' + + 'over a year and more in any single month. Shading from trees and buildings is not modeled beyond the ' + + 'Shading loss.'; + +function LocationSection({get, set, s, climatology}) { + const latError = get('lat') && s.latitude === null ? 'Latitude must be between -90 and 90' : null; + const lonError = get('lon') && s.longitude === null ? 'Longitude must be between -180 and 180' : null; + + let source; + if (s.ownData) { + source = 'Using your own sunlight data.'; + } else if (climatology.loading) { + source = ; + } else if (climatology.error) { + source = + Enter your own monthly sunlight under "Use my own sunlight data" below. + ; + } else { + source = `Sunlight and temperature: ${CLIMATOLOGY_SOURCE}.`; + } + + return <> +
1. Location
+ + + set({lat: e.target.value})}/> + + + set({lon: e.target.value})}/> + + +
{source}
+ ; +} + +function PanelsSection({setNumber, s}) { + return <> +
2. Your panels
+ + + setNumber('pw', v, DEFAULT_PANEL_WATTS)}/> + + + setNumber('pn', v, DEFAULT_PANEL_COUNT)}/> + + + + {fmt(s.kwp * 1000, 0)} W total ({fmt(s.kwp)} kW) + + + + ; +} + +function Results({result, s}) { + const unit = s.inverterOn ? 'kWh AC' : 'kWh DC'; + const rows = result.months.map(m => + {m.month} + {fmt(m.energyDaily)} + {fmt(m.energy, 0)} + {fmt(m.poa)} + {fmt(m.dayLength, 1)} h + ); + + return +
Estimate
+ + + + + + + + + + Month + {unit}/day + {unit}/month + Sun hours on panel + Day length + + + {rows} +
+ ; +} + +function AngleSection({setNumber, s, sunlight, ready}) { + const compass = COMPASS.find(c => Number(c.value) === s.azimuth)?.value || null; + const applyOptimal = goal => setNumber('tilt', optimalTilt(modelInputs(s, sunlight), goal), s.defaultTilt); + const azimuthDefault = defaultAzimuth(s.latitude ?? 0); + + if (s.tracking !== 'fixed') { + return + A tracking array turns to follow the sun, so it has no fixed angle. Change Tracking under Panel details + to set one. + ; + } + + return <> + + Leave these alone unless your panels are already mounted. The default (tilted to your latitude, facing + the equator) is close to the best for the year. "Best for winter" suits off-grid systems that must get + through the darkest month. + + + + setNumber('tilt', v, s.defaultTilt)}/> + + + + + + + + setNumber('az', v, azimuthDefault)}/> + ({value, label: m.label})), + {value: 'custom', label: 'Custom (from datasheet)'}, + ]} + value={moduleType} allowDeselect={false} + onChange={v => v !== 'custom' && setNumber('gamma', MODULE_TYPES[v].gamma, + MODULE_TYPES.standard.gamma)}/> + + + setNumber('gamma', v, MODULE_TYPES.standard.gamma)}/> + + + ({value, label: t.label}))} + value={s.tracking} allowDeselect={false} + onChange={v => set({track: v === 'fixed' ? null : v})}/> + + + set({alb: v === '' ? null : v})}/> + + + ; +} + +function LossesSection({set, s}) { + const setLoss = (index, value) => { + const next = [...s.losses]; + next[index] = value === '' ? LOSS_CATEGORIES[index].default : Number(value); + const isDefault = next.every((v, i) => v === DEFAULT_LOSSES[i]); + set({losses: isDefault ? null : formatMonthly(next)}); + }; + + return <> + + Real systems lose some power to dirt, wiring, shading and more. These are PVWatts' defaults. Raise Shading + if trees or buildings shade your panels, or Snow if snow sits on them for weeks. + + + {LOSS_CATEGORIES.map((c, i) => + setLoss(i, v)}/> + )} + + + ; +} + +function InverterSection({set, setNumber, s}) { + return <> + + Leave this off if your panels charge batteries through a charge controller. Turn it on to see AC output + from a grid-tie or hybrid inverter. + + set({inv: e.currentTarget.checked ? '1' : null})}/> + {s.inverterOn && + + setNumber('iac', v, roundDigits(s.kwp / 1.2, 2))}/> + + + setNumber('ieff', v, INVERTER_NOMINAL_EFFICIENCY)}/> + + } + ; +} + +function OwnDataSection({set, s, nasa}) { + const unit = s.fahrenheit ? '°F' : '°C'; + + // Start from the NASA values, when available, so the user edits rather than types all 24. + const enable = on => { + if (!on) { + set({ghi: null, temp: null}); + return; + } + const ghi = nasa ? nasa.ghi.map(v => (v === null ? '' : roundDigits(v, 2))) : Array(12).fill(''); + const temp = nasa + ? nasa.temperature.map(v => (v === null ? '' : roundDigits(s.fahrenheit ? cToF(v) : v, 1))) + : Array(12).fill(''); + // An all-blank list is still written for sunlight, so the section stays on for typing. + set({ghi: formatMonthly(ghi), temp: nasa ? formatMonthly(temp) : null}); + }; + + const setMonth = (key, values, index, value) => { + const next = [...values]; + next[index] = value === '' ? null : value; + set({[key]: formatMonthly(next)}); + }; + + const rows = MONTH_NAMES.map((name, i) => + {name} + + setMonth('ghi', s.ghi, i, v)}/> + + + setMonth('temp', s.temperatureInput, i, v)}/> + + ); + + return <> + + Only needed if you have better local measurements than the built-in data, for example in mountains, where + a 1° grid cell (about 69 miles across) can average very different terrain. Sunlight is global horizontal + irradiance in kWh/m²/day, the same as "peak sun hours" on flat ground. Blank temperatures use 20 °C. + + + enable(e.currentTarget.checked)}/> + {s.ownData && set({tu: e.currentTarget.checked ? 'f' : null})}/>} + + {s.ownData && + + + Month + Sunlight (kWh/m²/day) + Avg. temp ({unit}) + + + {rows} +
} + ; +} + +export function SolarCalculator() { + const {searchParams, get, set, setNumber} = useSolarParams(); + const climatology = useClimatology(); + const queryString = searchParams.toString(); + // eslint-disable-next-line react-hooks/exhaustive-deps + const s = React.useMemo(() => systemFromParams(get), [queryString]); + const nasa = React.useMemo( + () => (climatology.grid && s.latitude !== null && s.longitude !== null + ? lookupClimatology(climatology.grid, s.latitude, s.longitude) : null), + [climatology.grid, s.latitude, s.longitude]); + const sunlight = s.ownData ? sunlightFor(s, null) : nasa; + const missing = missingInputs(s, sunlight); + const ready = missing.length === 0; + const result = React.useMemo( + () => (ready ? estimateProduction(modelInputs(s, sunlight)) : null), + // eslint-disable-next-line react-hooks/exhaustive-deps + [s, sunlight, ready]); + + const angleTitle = s.tracking === 'fixed' + ? `Panel angle: ${s.tilt}° tilt, facing ${compassLabel(s.azimuth)}` + : 'Panel angle: follows the sun'; + const detailsTitle = `Panel details: ${TRACKING[s.tracking].label.toLowerCase()}, ` + + `${MOUNTINGS[s.mounting].label.toLowerCase()}`; + const waitingForData = climatology.loading && !s.ownData; + + return
+
Solar
+

+ Estimate what solar panels will produce over a typical year. Enter a location and your panels; everything + else has sensible defaults. Settings are saved in the page address, so you can bookmark or share them. +

+ + + + + {ready && } + {!ready && !waitingForData && + Enter {missing.join(', ')} to see an estimate.} + +
Fine-tune (optional)
+ + + {angleTitle} + + + + + + {detailsTitle} + + + + System losses: {fmt(totalLoss(s.losses))}% + + + + Inverter: {s.inverterOn ? 'on (AC output)' : 'off (DC output)'} + + + + Use my own sunlight data: {s.ownData ? 'on' : 'off'} + + + +
; +} diff --git a/app/src/components/calculators/SolarCalculator.test.js b/app/src/components/calculators/SolarCalculator.test.js new file mode 100644 index 00000000..58d68abf --- /dev/null +++ b/app/src/components/calculators/SolarCalculator.test.js @@ -0,0 +1,190 @@ +import React from "react"; +import {fireEvent, screen, within} from "@testing-library/react"; +import {renderWithProviders} from "../../test-utils"; +import {missingInputs, SolarCalculator, sunlightFor, systemFromParams} from "./SolarCalculator"; +import {CalculatorsPage} from "../Calculators"; +import {DEFAULT_LOSSES, estimateProduction, MODULE_TYPES} from "./solar"; +import {decodeClimatology, forgetClimatology, lookupClimatology} from "./solarData"; +import {mockClimatologyFetch, SUNNY_GHI, uniformClimatology} from "./solarData.fixtures"; +import {roundDigits} from "../Common"; + +const BASE = '/more/calculators?calc=solar'; +const DENVER = `${BASE}&lat=40&lon=-105`; + +const params = query => { + const searchParams = new URLSearchParams(query); + return key => searchParams.get(key); +}; + +const currentParams = () => new URLSearchParams(window.location.search); + +// The estimate appears once the data file has loaded and decoded, which can be slow on a busy machine. +const findEstimate = (unit = 'DC') => screen.findByText(`kWh ${unit} per year`, {}, {timeout: 5000}); + +// The annual estimate the page should show for a location, from the same uniform test data. +const expectedAnnual = (overrides = {}) => { + const grid = decodeClimatology(uniformClimatology()); + const sunlight = lookupClimatology(grid, 40, -105); + return estimateProduction({latitude: 40, kwp: 0.4, tilt: 40, azimuth: 180, ...sunlight, ...overrides}).annual; +}; + +beforeEach(() => { + forgetClimatology(); + mockClimatologyFetch(); +}); + +describe('systemFromParams', () => { + test('defaults when the URL is empty', () => { + const s = systemFromParams(params('')); + expect(s.latitude).toBeNull(); + expect(s.kwp).toBeCloseTo(0.4, 6); + expect(s.tracking).toBe('fixed'); + expect(s.mounting).toBe('open'); + expect(s.gamma).toBe(MODULE_TYPES.standard.gamma); + expect(s.losses).toEqual(DEFAULT_LOSSES); + expect(s.inverterOn).toBe(false); + expect(s.ownData).toBe(false); + expect(s.albedo).toBeNull(); + }); + + test('tilt and azimuth default from the latitude', () => { + expect(systemFromParams(params('lat=40.4')).tilt).toBe(40); + expect(systemFromParams(params('lat=40.4')).azimuth).toBe(180); + expect(systemFromParams(params('lat=-33.9')).tilt).toBe(34); + expect(systemFromParams(params('lat=-33.9')).azimuth).toBe(0); + }); + + test('invalid values fall back instead of crashing', () => { + const s = systemFromParams(params('lat=999&track=bogus&mount=x&tilt=abc&losses=1,2,abc&alb=7&pn=-3')); + expect(s.latitude).toBeNull(); + expect(s.tracking).toBe('fixed'); + expect(s.mounting).toBe('open'); + expect(s.tilt).toBe(0); + expect(s.losses.slice(0, 3)).toEqual([1, 2, DEFAULT_LOSSES[2]]); + expect(s.albedo).toBe(1); + expect(s.kwp).toBe(0); + }); + + test('own Fahrenheit temperatures are converted to Celsius', () => { + const s = systemFromParams(params('ghi=5&tu=f&temp=32,212')); + expect(s.ownData).toBe(true); + expect(s.temperature[0]).toBeCloseTo(0, 6); + expect(s.temperature[1]).toBeCloseTo(100, 6); + expect(s.temperatureInput[0]).toBe(32); + }); +}); + +describe('sunlight and missing inputs', () => { + const grid = decodeClimatology(uniformClimatology()); + + test('the NASA data is used unless the user gave their own', () => { + const nasa = sunlightFor(systemFromParams(params('lat=40&lon=-105')), grid); + expect(nasa.ghi[5]).toBeCloseTo(SUNNY_GHI[5], 1); + expect(nasa.diffuse[5]).toBeCloseTo(SUNNY_GHI[5] * 0.3, 1); + + const own = sunlightFor(systemFromParams(params(`lat=40&lon=-105&ghi=${Array(12).fill(1).join(',')}`)), grid); + expect(own.ghi).toEqual(Array(12).fill(1)); + expect(own.diffuse).toBeNull(); + }); + + test('a location is enough with the NASA data', () => { + const s = systemFromParams(params('lat=40&lon=-105')); + expect(missingInputs(s, sunlightFor(s, grid))).toEqual([]); + const noLon = systemFromParams(params('lat=40')); + expect(missingInputs(noLon, sunlightFor(noLon, grid))).toEqual(['a longitude']); + }); + + test('own data needs all 12 months but no longitude', () => { + const s = systemFromParams(params('lat=40&ghi=1,2,3')); + expect(missingInputs(s, sunlightFor(s, grid))).toEqual(['sunlight for all 12 months']); + }); +}); + +describe('SolarCalculator', () => { + test('is linked from the Engineering group of the calculators page', () => { + renderWithProviders(, {route: '/more/calculators'}); + const engineering = screen.getByRole('heading', {name: 'Engineering'}).closest('.wrolpi-header').parentElement; + const link = within(engineering).getByRole('link', {name: /Solar/}); + expect(link).toHaveAttribute('href', '/more/calculators?calc=solar'); + }); + + test('asks for a location', async () => { + renderWithProviders(, {route: BASE}); + expect(await screen.findByText('Enter a latitude, a longitude to see an estimate.', {}, {timeout: 5000})) + .toBeInTheDocument(); + }); + + test('a location alone gives an estimate from the built-in data', async () => { + renderWithProviders(, {route: DENVER}); + expect(await findEstimate()).toBeInTheDocument(); + expect(screen.getByText(`${roundDigits(expectedAnnual(), 0)}`)).toBeInTheDocument(); + expect(screen.getByText(/Sunlight and temperature: NASA POWER/)).toBeInTheDocument(); + expect(global.fetch).toHaveBeenCalledWith(expect.stringMatching(/\/data\/solar-climatology\.bin$/)); + }); + + test('typing writes the URL and leaves defaults out', async () => { + renderWithProviders(, {route: BASE}); + fireEvent.change(screen.getByLabelText('Latitude'), {target: {value: '40°0\'54"N'}}); + expect(currentParams().get('lat')).toBe('40°0\'54"N'); + expect(currentParams().get('calc')).toBe('solar'); + + fireEvent.change(screen.getByLabelText('Panels'), {target: {value: '4'}}); + expect(currentParams().get('pn')).toBe('4'); + fireEvent.change(screen.getByLabelText('Panels'), {target: {value: '1'}}); + expect(currentParams().has('pn')).toBe(false); + }); + + test('an invalid latitude is flagged', () => { + renderWithProviders(, {route: `${BASE}&lat=95`}); + expect(screen.getByText('Latitude must be between -90 and 90')).toBeInTheDocument(); + }); + + test('the inverter switches the results to AC', async () => { + renderWithProviders(, {route: DENVER}); + await findEstimate(); + fireEvent.click(screen.getByText('Inverter: off (DC output)')); + fireEvent.click(screen.getByLabelText('Include an inverter (AC output)')); + expect(currentParams().get('inv')).toBe('1'); + expect(screen.getByText('kWh AC per year')).toBeInTheDocument(); + }); + + test('best-for-the-year replaces a poor tilt', async () => { + renderWithProviders(, {route: `${DENVER}&tilt=10`}); + await findEstimate(); + expect(screen.getByText('Panel angle: 10° tilt, facing S')).toBeInTheDocument(); + fireEvent.click(screen.getByText('Panel angle: 10° tilt, facing S')); + fireEvent.click(await screen.findByRole('button', {name: 'Best for the year'})); + // A result of 40 equals the default (the latitude), so it is left out of the URL. + const tilt = Number(currentParams().get('tilt') ?? 40); + expect(tilt).toBeGreaterThanOrEqual(30); + expect(tilt).toBeLessThanOrEqual(44); + }); + + test('own sunlight data starts from the built-in values and replaces them', async () => { + renderWithProviders(, {route: DENVER}); + await findEstimate(); + fireEvent.click(screen.getByText('Use my own sunlight data: off')); + fireEvent.click(screen.getByLabelText('Use my own sunlight data')); + + const ghi = currentParams().get('ghi').split(',').map(Number); + ghi.forEach((v, i) => expect(v).toBeCloseTo(SUNNY_GHI[i], 1)); + expect(screen.getByText('Using your own sunlight data.')).toBeInTheDocument(); + expect(screen.getByLabelText('Jun sunlight')).toHaveValue(`${ghi[5]}`); + + // Halving June's sunlight lowers the estimate. + fireEvent.change(screen.getByLabelText('Jun sunlight'), {target: {value: `${ghi[5] / 2}`}}); + expect(currentParams().get('ghi').split(',')[5]).toBe(`${ghi[5] / 2}`); + + fireEvent.click(screen.getByLabelText('Use my own sunlight data')); + expect(currentParams().has('ghi')).toBe(false); + expect(currentParams().has('temp')).toBe(false); + }); + + test('a failed download explains how to continue', async () => { + forgetClimatology(); + global.fetch = jest.fn(() => Promise.resolve({ok: false, status: 404})); + renderWithProviders(, {route: DENVER}); + expect(await screen.findByText('Sunlight data could not be loaded', {}, {timeout: 5000})).toBeInTheDocument(); + expect(screen.getByText('Enter a location with sunlight data to see an estimate.')).toBeInTheDocument(); + }); +}); diff --git a/app/src/components/calculators/solar.js b/app/src/components/calculators/solar.js new file mode 100644 index 00000000..06663ddb --- /dev/null +++ b/app/src/components/calculators/solar.js @@ -0,0 +1,472 @@ +// Pure calculation functions for the Solar calculator. No React here; SolarCalculator.js is the UI. +// +// The model estimates a typical year from 12 monthly averages, the standard approach when hourly +// typical-year weather files are not available: +// +// 1. Sun geometry for Klein's representative day of each month (Duffie & Beckman ch. 1). +// 2. Monthly global horizontal irradiation (GHI) is split into diffuse and beam (Erbs monthly +// correlation), then spread across the day (Collares-Pereira & Rabl for total, Liu & Jordan +// for diffuse). +// 3. Each step is transposed onto the panel with the HDKR sky model, for a fixed, single-axis or +// dual-axis array. +// 4. PVWatts' models turn plane-of-array irradiance into energy: ASHRAE angle-of-incidence loss, +// SAPM cell temperature, linear temperature-corrected DC power, multiplicative system losses, +// and an optional inverter efficiency curve with clipping. +// +// Angles are in degrees at the API boundary. Azimuths are compass bearings (0 = north, 90 = east, +// 180 = south, 270 = west); internally the Duffie & Beckman convention is used (0 = south, east +// negative, west positive). Energy is kWh, irradiation kWh/m², irradiance W/m², power kW. + +export const SOLAR_CONSTANT = 1361; // W/m² + +const RAD = Math.PI / 180; +const sin = deg => Math.sin(deg * RAD); +const cos = deg => Math.cos(deg * RAD); +const tan = deg => Math.tan(deg * RAD); +const acos = x => Math.acos(Math.min(1, Math.max(-1, x))) / RAD; +const clamp = (x, lo, hi) => Math.min(hi, Math.max(lo, x)); + +export const MONTH_NAMES = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']; +export const DAYS_IN_MONTH = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]; + +// Klein's recommended average day of each month: the day whose extraterrestrial radiation is +// closest to the month's mean. Day of the year, 1-365. +export const KLEIN_DAYS = [17, 47, 75, 105, 135, 162, 198, 228, 258, 288, 318, 344]; + +// --------------------------------------------------------------------------- +// Sun geometry +// --------------------------------------------------------------------------- + +// Solar declination (Cooper), degrees. +export function declination(dayOfYear) { + return 23.45 * sin(360 * (284 + dayOfYear) / 365); +} + +// Sunset hour angle on a horizontal surface, degrees. 0 during polar night, 180 during polar day. +export function sunsetHourAngle(latitude, decl) { + return acos(-tan(latitude) * tan(decl)); +} + +// Orbit eccentricity correction to the solar constant. +export function eccentricity(dayOfYear) { + return 1 + 0.033 * cos(360 * dayOfYear / 365); +} + +// Daily extraterrestrial irradiation on a horizontal surface, kWh/m²/day (D&B eq. 1.10.3). +export function extraterrestrialDaily(latitude, dayOfYear) { + const decl = declination(dayOfYear); + const ws = sunsetHourAngle(latitude, decl); + const wh = (24 / Math.PI) * SOLAR_CONSTANT * eccentricity(dayOfYear) + * (cos(latitude) * cos(decl) * sin(ws) + (Math.PI * ws / 180) * sin(latitude) * sin(decl)); + return Math.max(0, wh / 1000); +} + +// Day length in hours. +export function dayLength(latitude, decl) { + return 2 * sunsetHourAngle(latitude, decl) / 15; +} + +// Sun elevation at solar noon, degrees. Negative means the sun does not rise. +export function noonElevation(latitude, decl) { + return 90 - Math.abs(latitude - decl); +} + +// Sun position for an hour angle (degrees, 15° per hour from solar noon, afternoon positive). +// Returns cosine of the zenith angle and the azimuth in the D&B convention (0 = south). +export function sunPosition(latitude, decl, hourAngle) { + const cosZenith = cos(latitude) * cos(decl) * cos(hourAngle) + sin(latitude) * sin(decl); + const sinZenith = Math.sqrt(Math.max(0, 1 - cosZenith * cosZenith)); + let azimuth = 0; + if (sinZenith > 1e-9 && Math.abs(cos(latitude)) > 1e-9) { + azimuth = acos((cosZenith * sin(latitude) - sin(decl)) / (sinZenith * cos(latitude))); + azimuth = hourAngle < 0 ? -azimuth : azimuth; + } + return {cosZenith, zenith: acos(cosZenith), azimuth}; +} + +// Cosine of the angle of incidence on a surface of tilt `tilt` and D&B azimuth `surfaceAzimuth`. +export function cosIncidence(zenith, sunAzimuth, tilt, surfaceAzimuth) { + return cos(zenith) * cos(tilt) + sin(zenith) * sin(tilt) * cos(sunAzimuth - surfaceAzimuth); +} + +// Compass bearing (0 = north) to the D&B azimuth convention (0 = south, east negative). +export function compassToSolarAzimuth(compass) { + let a = ((compass % 360) + 360) % 360 - 180; + return a === -180 ? 180 : a; +} + +// Default panel azimuth: face the equator. +export function defaultAzimuth(latitude) { + return latitude < 0 ? 0 : 180; +} + +// --------------------------------------------------------------------------- +// Monthly irradiation to sub-daily steps +// --------------------------------------------------------------------------- + +// Monthly average diffuse fraction Hd/H from the clearness index (Erbs, Klein & Duffie 1982). +// The correlation was fit for 0.3 <= KT <= 0.8, so KT is clamped to that range. +export function erbsMonthlyDiffuseFraction(kt, sunsetAngle) { + const k = clamp(kt, 0.3, 0.8); + if (sunsetAngle <= 81.4) { + return 1.391 - 3.560 * k + 4.189 * k * k - 2.137 * k * k * k; + } + return 1.311 - 3.022 * k + 3.427 * k * k - 1.821 * k * k * k; +} + +// Ratio of hourly to daily total irradiation at hour angle `w` (Collares-Pereira & Rabl). +export function totalHourlyRatio(w, sunsetAngle) { + if (Math.abs(w) >= sunsetAngle) { + return 0; + } + const a = 0.409 + 0.5016 * sin(sunsetAngle - 60); + const b = 0.6609 - 0.4767 * sin(sunsetAngle - 60); + return diffuseHourlyRatio(w, sunsetAngle) * (a + b * cos(w)); +} + +// Ratio of hourly to daily diffuse irradiation at hour angle `w` (Liu & Jordan). +export function diffuseHourlyRatio(w, sunsetAngle) { + if (Math.abs(w) >= sunsetAngle) { + return 0; + } + const denominator = sin(sunsetAngle) - (Math.PI * sunsetAngle / 180) * cos(sunsetAngle); + return (Math.PI / 24) * (cos(w) - cos(sunsetAngle)) / denominator; +} + +// Number of evenly spaced steps between sunrise and sunset. Steps adapt to the day length, so a +// short polar-edge day is still sampled. +export const STEPS_PER_DAY = 48; + +// The average day of a month as a list of steps, each with its hour angle, duration (hours), and +// the horizontal irradiances (W/m²) for global, beam and diffuse. The two hourly ratios are +// normalized over the day so the steps add back up to exactly the monthly GHI and diffuse given. +export function averageDaySteps(latitude, monthIndex, ghi, diffuse = null) { + const n = KLEIN_DAYS[monthIndex]; + const decl = declination(n); + const ws = sunsetHourAngle(latitude, decl); + const h0 = extraterrestrialDaily(latitude, n); + if (!(ghi > 0) || !(ws > 0) || !(h0 > 0)) { + return {decl, sunsetAngle: ws, kt: null, diffuseFraction: null, steps: []}; + } + + const kt = ghi / h0; + const hd = diffuse != null && diffuse >= 0 + ? Math.min(diffuse, ghi) + : ghi * erbsMonthlyDiffuseFraction(kt, ws); + + const dw = 2 * ws / STEPS_PER_DAY; + const dt = dw / 15; + const raw = []; + let totalSum = 0, diffuseSum = 0; + for (let i = 0; i < STEPS_PER_DAY; i++) { + const w = -ws + (i + 0.5) * dw; + const rt = totalHourlyRatio(w, ws); + const rd = diffuseHourlyRatio(w, ws); + raw.push({w, rt, rd}); + totalSum += rt * dt; + diffuseSum += rd * dt; + } + + const steps = raw.map(({w, rt, rd}) => { + // kWh/m² per hour is numerically kW/m², so x1000 gives the step's mean irradiance in W/m². + const global = totalSum > 0 ? (rt / totalSum) * ghi * 1000 : 0; + const diff = Math.min(global, diffuseSum > 0 ? (rd / diffuseSum) * hd * 1000 : 0); + return {hourAngle: w, hours: dt, global, diffuse: diff, beam: global - diff}; + }); + return {decl, sunsetAngle: ws, kt, diffuseFraction: hd / ghi, steps}; +} + +// --------------------------------------------------------------------------- +// Panel orientation and plane-of-array irradiance +// --------------------------------------------------------------------------- + +export const TRACKING = { + fixed: {label: 'Fixed'}, + single: {label: 'Single-axis (N–S axis)'}, + dual: {label: 'Dual-axis'}, +}; + +// Single-axis trackers are limited to this rotation either side of flat. +export const SINGLE_AXIS_LIMIT = 60; + +// Panel tilt and D&B azimuth for a step. `tilt`/`azimuth` are only used for fixed arrays. +export function panelOrientation(tracking, sun, tilt, azimuth) { + if (tracking === 'dual') { + return {tilt: sun.zenith, azimuth: sun.azimuth}; + } + if (tracking === 'single') { + // Horizontal north-south axis rotating east-west (Marion & Dobos ideal rotation). + const rotation = clamp( + Math.atan(Math.tan(sun.zenith * RAD) * Math.sin(sun.azimuth * RAD)) / RAD, + -SINGLE_AXIS_LIMIT, SINGLE_AXIS_LIMIT); + return {tilt: Math.abs(rotation), azimuth: rotation >= 0 ? 90 : -90}; + } + return {tilt, azimuth}; +} + +// Plane-of-array irradiance for one step with the HDKR (Hay-Davies-Klucher-Reindl) sky model. +// Returns the beam (direct plus circumsolar), sky diffuse and ground-reflected parts in W/m², and +// the angle of incidence. +export function poaIrradiance({global, beam, diffuse}, sun, orientation, albedo, dayOfYear) { + const {tilt, azimuth} = orientation; + const cosTheta = cosIncidence(sun.zenith, sun.azimuth, tilt, azimuth); + + // Beam normal irradiance, capped at the extraterrestrial value; near sunrise the division by a + // tiny cos(zenith) would otherwise blow up. + const extraNormal = SOLAR_CONSTANT * eccentricity(dayOfYear); + const cosZ = Math.max(sun.cosZenith, 0.01745); + const dni = Math.min(beam / cosZ, extraNormal); + const rb = Math.max(0, cosTheta) / cosZ; + + const anisotropy = clamp(dni / extraNormal, 0, 1); + const horizonFactor = global > 0 ? Math.sqrt(beam / global) : 0; + const beamPoa = dni * Math.max(0, cosTheta) + diffuse * anisotropy * rb; + const skyPoa = diffuse * (1 - anisotropy) * ((1 + cos(tilt)) / 2) + * (1 + horizonFactor * Math.pow(sin(tilt / 2), 3)); + const groundPoa = global * albedo * (1 - cos(tilt)) / 2; + return {beam: beamPoa, sky: skyPoa, ground: groundPoa, total: beamPoa + skyPoa + groundPoa, cosTheta}; +} + +// ASHRAE incidence angle modifier for the glass cover; applied to the beam part. +export function ashraeIam(cosTheta, b0 = 0.05) { + if (cosTheta <= 0) { + return 0; + } + return clamp(1 - b0 * (1 / cosTheta - 1), 0, 1); +} + +// --------------------------------------------------------------------------- +// Temperature, DC power, losses, inverter +// --------------------------------------------------------------------------- + +// Sandia (SAPM) cell temperature coefficients by mounting. +export const MOUNTINGS = { + open: {label: 'Open rack (ground or pole)', a: -3.56, b: -0.075, deltaT: 3}, + roof: {label: 'Close roof mount', a: -2.98, b: -0.0471, deltaT: 1}, + flush: {label: 'Insulated back / flush', a: -2.81, b: -0.0455, deltaT: 0}, +}; + +// PVWatts v8 temperature coefficients of power, %/°C. +export const MODULE_TYPES = { + standard: {label: 'Standard (crystalline)', gamma: -0.37}, + premium: {label: 'Premium (high efficiency)', gamma: -0.35}, + thinFilm: {label: 'Thin film', gamma: -0.32}, +}; + +export const DEFAULT_WIND_SPEED = 1; // m/s; deliberately low, so cell temperatures are not optimistic. + +// SAPM cell temperature, °C. +export function cellTemperature(poa, ambient, mounting = 'open', windSpeed = DEFAULT_WIND_SPEED) { + const m = MOUNTINGS[mounting] || MOUNTINGS.open; + const moduleTemp = poa * Math.exp(m.a + m.b * windSpeed) + ambient; + return moduleTemp + (poa / 1000) * m.deltaT; +} + +// Daily air temperature swing, used when only a monthly mean is known: the mean +/- this amount. +export const DEFAULT_DAILY_SWING = 5; + +// Air temperature at an hour angle, peaking mid-afternoon (3 pm). +export function ambientTemperature(mean, hourAngle, swing = DEFAULT_DAILY_SWING) { + return mean + swing * cos(hourAngle - 45); +} + +// DC power, kW, for an array of `kwp` kilowatts (STC rating) and temperature coefficient `gamma` (%/°C). +export function dcPower(effectiveIrradiance, cellTemp, kwp, gamma) { + return Math.max(0, kwp * (effectiveIrradiance / 1000) * (1 + (gamma / 100) * (cellTemp - 25))); +} + +// PVWatts system loss categories and their default percentages. +export const LOSS_CATEGORIES = [ + {key: 'soiling', label: 'Soiling', default: 2}, + {key: 'shading', label: 'Shading', default: 3}, + {key: 'snow', label: 'Snow', default: 0}, + {key: 'mismatch', label: 'Mismatch', default: 2}, + {key: 'wiring', label: 'Wiring', default: 2}, + {key: 'connections', label: 'Connections', default: 0.5}, + {key: 'lid', label: 'Light-induced degradation', default: 1.5}, + {key: 'nameplate', label: 'Nameplate rating', default: 1}, + {key: 'age', label: 'Age', default: 0}, + {key: 'availability', label: 'Availability', default: 3}, +]; + +export const DEFAULT_LOSSES = LOSS_CATEGORIES.map(c => c.default); + +// Combine loss percentages multiplicatively, as PVWatts does. Returns a percentage. +export function totalLoss(percentages) { + const kept = percentages.reduce((acc, p) => acc * (1 - clamp(Number(p) || 0, 0, 100) / 100), 1); + return 100 * (1 - kept); +} + +export const INVERTER_NOMINAL_EFFICIENCY = 96; +const INVERTER_REFERENCE_EFFICIENCY = 0.9637; + +// AC power, kW, from the PVWatts inverter model: part-load efficiency curve, clipped at the AC rating. +export function inverterOutput(pdc, acRating, nominalEfficiency = INVERTER_NOMINAL_EFFICIENCY) { + if (!(pdc > 0) || !(acRating > 0)) { + return 0; + } + const etaNom = nominalEfficiency / 100; + const pdc0 = acRating / etaNom; + const zeta = pdc / pdc0; + const eta = (etaNom / INVERTER_REFERENCE_EFFICIENCY) * (-0.0162 * zeta - 0.0059 / zeta + 0.9858); + return Math.min(Math.max(0, eta * pdc), acRating); +} + +// --------------------------------------------------------------------------- +// The whole model +// --------------------------------------------------------------------------- + +export const DEFAULT_ALBEDO = 0.2; + +/** + * Estimate monthly and annual production. + * + * @param {Object} system + * @param {number} system.latitude degrees, north positive + * @param {number[]} system.ghi 12 monthly average GHI values, kWh/m²/day + * @param {number[]} [system.diffuse] 12 monthly diffuse values, kWh/m²/day; estimated when absent + * @param {number[]} [system.temperature] 12 monthly mean air temperatures, °C; 20 when absent + * @param {number[]} [system.albedo] ground reflectance, one value or 12 monthly values + * @param {number} system.kwp array STC rating, kW + * @param {number} [system.tilt] degrees from horizontal (fixed arrays) + * @param {number} [system.azimuth] compass degrees (fixed arrays) + * @param {string} [system.tracking] 'fixed' | 'single' | 'dual' + * @param {string} [system.mounting] key of MOUNTINGS + * @param {number} [system.gamma] temperature coefficient of power, %/°C + * @param {number[]} [system.losses] system loss percentages, see LOSS_CATEGORIES + * @param {Object} [system.inverter] {acRating (kW), efficiency (%)}; null for DC output + */ +export function estimateProduction({ + latitude, + ghi, + diffuse = null, + temperature = null, + albedo = DEFAULT_ALBEDO, + kwp, + tilt = 0, + azimuth = defaultAzimuth(latitude), + tracking = 'fixed', + mounting = 'open', + gamma = MODULE_TYPES.standard.gamma, + losses = DEFAULT_LOSSES, + inverter = null, + }) { + const lossFactor = 1 - totalLoss(losses) / 100; + const surfaceAzimuth = compassToSolarAzimuth(azimuth); + + const months = MONTH_NAMES.map((name, m) => { + const n = KLEIN_DAYS[m]; + const givenAlbedo = Array.isArray(albedo) ? albedo[m] : albedo; + const monthAlbedo = Number.isFinite(givenAlbedo) ? givenAlbedo : DEFAULT_ALBEDO; + const meanTemp = temperature && Number.isFinite(temperature[m]) ? temperature[m] : 20; + const day = averageDaySteps(latitude, m, ghi[m], diffuse ? diffuse[m] : null); + + let poaDaily = 0, energyDaily = 0; + for (const step of day.steps) { + const sun = sunPosition(latitude, day.decl, step.hourAngle); + if (sun.cosZenith <= 0) { + continue; + } + const orientation = panelOrientation(tracking, sun, tilt, surfaceAzimuth); + const poa = poaIrradiance(step, sun, orientation, monthAlbedo, n); + const effective = poa.beam * ashraeIam(poa.cosTheta) + poa.sky + poa.ground; + const tc = cellTemperature(poa.total, ambientTemperature(meanTemp, step.hourAngle), mounting); + let power = dcPower(effective, tc, kwp, gamma) * lossFactor; + if (inverter) { + power = inverterOutput(power, inverter.acRating, inverter.efficiency); + } + poaDaily += poa.total * step.hours / 1000; + energyDaily += power * step.hours; + } + + return { + month: name, + ghi: ghi[m], + kt: day.kt, + diffuseFraction: day.diffuseFraction, + // Plane-of-array kWh/m²/day: the honest "peak sun hours" for this panel. + poa: poaDaily, + energyDaily, + energy: energyDaily * DAYS_IN_MONTH[m], + dayLength: dayLength(latitude, day.decl), + noonElevation: noonElevation(latitude, day.decl), + }; + }); + + const annual = months.reduce((acc, m) => acc + m.energy, 0); + const worst = months.reduce((a, b) => (b.energyDaily < a.energyDaily ? b : a)); + const best = months.reduce((a, b) => (b.energyDaily > a.energyDaily ? b : a)); + return { + months, + annual, + specificYield: kwp > 0 ? annual / kwp : null, + capacityFactor: kwp > 0 ? annual / (kwp * 8760) : null, + worst, + best, + totalLoss: totalLoss(losses), + }; +} + +// Sweep fixed tilts 0-90° at the given azimuth and return the best one. `goal` is 'annual' (most +// energy in a year) or 'winter' (most energy in the worst month, which sizes off-grid systems). +export function optimalTilt(system, goal = 'annual') { + let best = {tilt: 0, value: -1}; + for (let tilt = 0; tilt <= 90; tilt++) { + const result = estimateProduction({...system, tilt, tracking: 'fixed'}); + const value = goal === 'winter' ? result.worst.energyDaily : result.annual; + if (value > best.value + 1e-9) { + best = {tilt, value}; + } + } + return best.tilt; +} + +// --------------------------------------------------------------------------- +// Input parsing +// --------------------------------------------------------------------------- + +// Parse latitude/longitude text: decimal ("40.015", "-105.27") or degrees-minutes-seconds +// ("40°0'54\"N", "105 16 12 W"). Returns null when unparseable or out of range. +export function parseCoordinate(text, kind = 'lat') { + if (text === null || text === undefined) { + return null; + } + const s = `${text}`.trim().toUpperCase(); + if (!s) { + return null; + } + const hemisphere = s.match(/[NSEW]/); + const parts = s.replace(/[NSEW]/g, ' ').match(/-?\d+(\.\d+)?/g); + if (!parts || parts.length > 3) { + return null; + } + const [d, m = 0, sec = 0] = parts.map(Number); + if (m < 0 || m >= 60 || sec < 0 || sec >= 60) { + return null; + } + let value = Math.abs(d) + m / 60 + sec / 3600; + if (d < 0 || s.trim().startsWith('-')) { + value = -value; + } + if (hemisphere) { + const h = hemisphere[0]; + if ((kind === 'lat' && (h === 'E' || h === 'W')) || (kind === 'lon' && (h === 'N' || h === 'S'))) { + return null; + } + if (h === 'S' || h === 'W') { + value = -Math.abs(value); + } + } + const limit = kind === 'lat' ? 90 : 180; + return Math.abs(value) <= limit ? value : null; +} + +// Parse a comma separated list of 12 monthly numbers. Blank entries become null. +export function parseMonthly(text) { + const values = (text || '').split(',').map(v => (v.trim() === '' ? null : Number(v))); + return MONTH_NAMES.map((_, i) => (Number.isFinite(values[i]) ? values[i] : null)); +} + +export function formatMonthly(values) { + return values.map(v => (v === null || v === undefined || v === '' ? '' : `${v}`)).join(','); +} diff --git a/app/src/components/calculators/solar.test.js b/app/src/components/calculators/solar.test.js new file mode 100644 index 00000000..d24ed14a --- /dev/null +++ b/app/src/components/calculators/solar.test.js @@ -0,0 +1,310 @@ +import { + ambientTemperature, + ashraeIam, + averageDaySteps, + cellTemperature, + compassToSolarAzimuth, + DAYS_IN_MONTH, + dcPower, + declination, + DEFAULT_LOSSES, + defaultAzimuth, + diffuseHourlyRatio, + erbsMonthlyDiffuseFraction, + estimateProduction, + extraterrestrialDaily, + formatMonthly, + inverterOutput, + KLEIN_DAYS, + optimalTilt, + parseCoordinate, + parseMonthly, + poaIrradiance, + sunPosition, + sunsetHourAngle, + totalHourlyRatio, + totalLoss, +} from "./solar"; + +// A sunny 40°N GHI profile, kWh/m²/day, with a clearness index of about 0.6 year-round (Denver-like). +const GHI_40N = [2.6, 3.5, 4.7, 5.8, 6.5, 7.2, 7.0, 6.2, 5.2, 4.0, 2.9, 2.3]; +const TEMP_40N = [0, 2, 6, 10, 15, 21, 24, 23, 18, 11, 4, 0]; + +const system = (overrides = {}) => ({ + latitude: 40, + ghi: GHI_40N, + temperature: TEMP_40N, + kwp: 1, + tilt: 40, + azimuth: 180, + ...overrides, +}); + +const integrate = (ratio, ws) => { + // Midpoint rule over the day, ratios are per hour. + const n = 2000, dw = 2 * ws / n; + let sum = 0; + for (let i = 0; i < n; i++) { + sum += ratio(-ws + (i + 0.5) * dw, ws) * dw / 15; + } + return sum; +}; + +describe('sun geometry', () => { + test('declination matches Cooper at the solstices and equinox', () => { + expect(declination(172)).toBeCloseTo(23.45, 1); // June 21 + expect(declination(355)).toBeCloseTo(-23.45, 1); // December 21 + expect(Math.abs(declination(81))).toBeLessThan(0.5); // March 22 + // Klein's February day (D&B Table 1.6.1 lists -13.0). + expect(declination(KLEIN_DAYS[1])).toBeCloseTo(-12.95, 1); + }); + + test('sunset hour angle is 90° at the equator and handles polar day and night', () => { + expect(sunsetHourAngle(0, 20)).toBeCloseTo(90, 5); + expect(sunsetHourAngle(40, 0)).toBeCloseTo(90, 5); + expect(sunsetHourAngle(80, 20)).toBe(180); + expect(sunsetHourAngle(80, -20)).toBe(0); + }); + + test('extraterrestrial daily irradiation matches Duffie & Beckman Example 1.10.1', () => { + // 43°N on April 15: H0 = 33.8 MJ/m² = 9.39 kWh/m². + expect(extraterrestrialDaily(43, 105)).toBeCloseTo(33.8 / 3.6, 1); + expect(extraterrestrialDaily(80, 355)).toBe(0); // polar night + }); + + test('sun is due south at solar noon and east in the morning (northern hemisphere)', () => { + const noon = sunPosition(40, 0, 0); + expect(noon.azimuth).toBeCloseTo(0, 5); + expect(noon.zenith).toBeCloseTo(40, 5); + expect(sunPosition(40, 0, -45).azimuth).toBeLessThan(0); + expect(sunPosition(40, 0, 45).azimuth).toBeGreaterThan(0); + }); + + test('compass azimuth converts to the Duffie & Beckman convention', () => { + expect(compassToSolarAzimuth(180)).toBe(0); + expect(compassToSolarAzimuth(90)).toBe(-90); + expect(compassToSolarAzimuth(270)).toBe(90); + expect(compassToSolarAzimuth(0)).toBe(180); + expect(compassToSolarAzimuth(360)).toBe(180); + expect(defaultAzimuth(40)).toBe(180); + expect(defaultAzimuth(-33)).toBe(0); + }); +}); + +describe('monthly irradiation to steps', () => { + test('Erbs monthly diffuse fraction follows both branches and clamps KT', () => { + expect(erbsMonthlyDiffuseFraction(0.5, 80)).toBeCloseTo(0.3911, 4); + expect(erbsMonthlyDiffuseFraction(0.5, 100)).toBeCloseTo(1.311 - 1.511 + 0.85675 - 0.227625, 4); + expect(erbsMonthlyDiffuseFraction(0.95, 80)).toBeCloseTo(erbsMonthlyDiffuseFraction(0.8, 80), 10); + expect(erbsMonthlyDiffuseFraction(0.1, 80)).toBeCloseTo(erbsMonthlyDiffuseFraction(0.3, 80), 10); + }); + + test('Liu & Jordan diffuse ratio integrates to one day', () => { + for (const ws of [60, 90, 120, 180]) { + expect(integrate(diffuseHourlyRatio, ws)).toBeCloseTo(1, 3); + } + }); + + test('Collares-Pereira & Rabl total ratio integrates to about one day', () => { + // The correlation is empirical, so it only approximately conserves the daily total; the + // model normalizes it, but a large error here would mean a transcription mistake. + for (const ws of [60, 90, 120]) { + expect(Math.abs(integrate(totalHourlyRatio, ws) - 1)).toBeLessThan(0.03); + } + }); + + test('average day steps add back up to the monthly GHI and diffuse', () => { + const day = averageDaySteps(40, 5, 7.2); + const sum = key => day.steps.reduce((acc, s) => acc + s[key] * s.hours, 0) / 1000; + expect(sum('global')).toBeCloseTo(7.2, 6); + expect(sum('diffuse')).toBeCloseTo(7.2 * day.diffuseFraction, 6); + expect(day.steps.every(s => s.beam >= 0 && s.diffuse >= 0)).toBe(true); + // Peak at solar noon. + const peak = day.steps.reduce((a, b) => (b.global > a.global ? b : a)); + expect(Math.abs(peak.hourAngle)).toBeLessThan(5); + }); + + test('a supplied diffuse value is used instead of the Erbs estimate', () => { + expect(averageDaySteps(40, 5, 7.2, 1.8).diffuseFraction).toBeCloseTo(0.25, 6); + }); + + test('polar night and missing GHI produce no steps', () => { + expect(averageDaySteps(80, 11, 0.1).steps).toEqual([]); + expect(averageDaySteps(40, 5, null).steps).toEqual([]); + }); +}); + +describe('plane of array', () => { + test('a flat panel receives exactly the horizontal irradiance', () => { + const day = averageDaySteps(40, 5, 7.2); + const decl = day.decl; + for (const step of day.steps) { + const sun = sunPosition(40, decl, step.hourAngle); + if (sun.zenith > 85) continue; // the DNI cap applies right at sunrise + const poa = poaIrradiance(step, sun, {tilt: 0, azimuth: 0}, 0.2, KLEIN_DAYS[5]); + expect(poa.total).toBeCloseTo(step.global, 6); + } + }); + + test('a vertical panel sees half the ground', () => { + const step = {global: 500, beam: 0, diffuse: 500}; + const sun = sunPosition(40, 0, 0); + const poa = poaIrradiance(step, sun, {tilt: 90, azimuth: 0}, 0.2, 81); + expect(poa.ground).toBeCloseTo(500 * 0.2 / 2, 6); + }); + + test('ASHRAE incidence angle modifier', () => { + expect(ashraeIam(1)).toBe(1); + expect(ashraeIam(0.5)).toBeCloseTo(0.95, 6); // 60° + expect(ashraeIam(0)).toBe(0); + expect(ashraeIam(0.01)).toBe(0); // never negative + }); +}); + +describe('temperature, power, losses, inverter', () => { + test('SAPM cell temperature, open rack glass/polymer', () => { + // 1000 * exp(-3.56 - 0.075) + 25 = 51.38, plus deltaT 3. + expect(cellTemperature(1000, 25, 'open', 1)).toBeCloseTo(54.38, 2); + expect(cellTemperature(1000, 25, 'flush', 1)).toBeGreaterThan(cellTemperature(1000, 25, 'open', 1)); + expect(cellTemperature(0, 10, 'open')).toBe(10); + }); + + test('air temperature peaks at 3 pm and averages to the mean', () => { + expect(ambientTemperature(20, 45)).toBeCloseTo(25, 6); + expect(ambientTemperature(20, -135)).toBeCloseTo(15, 6); + }); + + test('DC power scales with irradiance and loses power when hot', () => { + expect(dcPower(1000, 25, 1, -0.37)).toBeCloseTo(1, 6); + expect(dcPower(500, 25, 2, -0.37)).toBeCloseTo(1, 6); + expect(dcPower(1000, 35, 1, -0.37)).toBeCloseTo(0.963, 6); + expect(dcPower(0, 25, 1, -0.37)).toBe(0); + }); + + test('PVWatts default losses combine to 14.08%', () => { + expect(totalLoss(DEFAULT_LOSSES)).toBeCloseTo(14.08, 2); + expect(totalLoss([])).toBe(0); + expect(totalLoss([50, 50])).toBeCloseTo(75, 6); + }); + + test('PVWatts inverter is at nominal efficiency at full load and clips', () => { + // At zeta = 1 the curve evaluates to exactly the nominal efficiency. + const ac = 1, pdc0 = ac / 0.96; + expect(inverterOutput(pdc0, ac) / pdc0).toBeCloseTo(0.96, 4); + expect(inverterOutput(5, ac)).toBe(ac); + expect(inverterOutput(0, ac)).toBe(0); + // Part load is less efficient. + expect(inverterOutput(pdc0 * 0.05, ac) / (pdc0 * 0.05)).toBeLessThan(0.96); + }); +}); + +describe('estimateProduction', () => { + test('a 1 kWp south-facing array at 40°N gives a plausible annual yield', () => { + const result = estimateProduction(system()); + // A clear 40°N site yields roughly 1500-1900 kWh/kWp DC after 14% losses. + expect(result.annual).toBeGreaterThan(1500); + expect(result.annual).toBeLessThan(1900); + // Performance ratio (energy over plane-of-array sun-hours): losses plus temperature. + const poaAnnual = result.months.reduce((acc, m, i) => acc + m.poa * DAYS_IN_MONTH[i], 0); + expect(result.annual / poaAnnual).toBeGreaterThan(0.75); + expect(result.annual / poaAnnual).toBeLessThan(0.88); + expect(result.specificYield).toBeCloseTo(result.annual, 6); + expect(result.months).toHaveLength(12); + expect(result.totalLoss).toBeCloseTo(14.08, 2); + }); + + test('tilting toward the equator gains in winter and loses in summer', () => { + const flat = estimateProduction(system({tilt: 0})); + const tilted = estimateProduction(system({tilt: 40})); + expect(tilted.months[11].poa).toBeGreaterThan(flat.months[11].poa * 1.5); + expect(tilted.months[5].poa).toBeLessThan(flat.months[5].poa); + expect(flat.months[5].poa).toBeCloseTo(GHI_40N[5], 1); + }); + + test('tracking beats fixed, and dual-axis beats single-axis', () => { + const fixed = estimateProduction(system()).annual; + const single = estimateProduction(system({tracking: 'single'})).annual; + const dual = estimateProduction(system({tracking: 'dual'})).annual; + expect(single).toBeGreaterThan(fixed); + expect(dual).toBeGreaterThan(single); + // Dual-axis is commonly 30-45% over a well-tilted fixed array. + expect(dual / fixed).toBeLessThan(1.6); + }); + + test('facing away from the equator produces less', () => { + const south = estimateProduction(system()).annual; + const north = estimateProduction(system({azimuth: 0})).annual; + expect(north).toBeLessThan(south * 0.7); + }); + + test('the southern hemisphere mirrors the northern', () => { + // Same GHI shifted 6 months, panel facing north. + const shifted = [...GHI_40N.slice(6), ...GHI_40N.slice(0, 6)]; + const north = estimateProduction(system({temperature: null})); + const south = estimateProduction(system({latitude: -40, ghi: shifted, azimuth: 0, temperature: null})); + expect(south.annual / north.annual).toBeGreaterThan(0.97); + expect(south.annual / north.annual).toBeLessThan(1.03); + }); + + test('hot weather reduces output', () => { + const cool = estimateProduction(system({temperature: Array(12).fill(0)})).annual; + const hot = estimateProduction(system({temperature: Array(12).fill(35)})).annual; + expect(hot).toBeLessThan(cool * 0.95); + }); + + test('the inverter costs a few percent and clips an oversized array', () => { + const dc = estimateProduction(system()).annual; + const ac = estimateProduction(system({inverter: {acRating: 1, efficiency: 96}})).annual; + expect(ac / dc).toBeGreaterThan(0.93); + expect(ac / dc).toBeLessThan(0.97); + const clipped = estimateProduction(system({inverter: {acRating: 0.3, efficiency: 96}})).annual; + expect(clipped).toBeLessThan(ac * 0.7); + }); + + test('polar night months produce nothing without errors', () => { + const result = estimateProduction(system({latitude: 75, tilt: 75})); + expect(result.months[11].energy).toBe(0); + expect(Number.isFinite(result.annual)).toBe(true); + }); + + test('missing months count as zero', () => { + const ghi = [...GHI_40N]; + ghi[0] = null; + expect(estimateProduction(system({ghi})).months[0].energy).toBe(0); + }); + + test('the best annual tilt is near the latitude; the best winter tilt is steeper', () => { + const annual = optimalTilt(system()); + const winter = optimalTilt(system(), 'winter'); + expect(annual).toBeGreaterThanOrEqual(30); + expect(annual).toBeLessThanOrEqual(44); + expect(winter).toBeGreaterThan(annual + 8); + }); +}); + +describe('input parsing', () => { + test('decimal and DMS coordinates', () => { + expect(parseCoordinate('40.015')).toBeCloseTo(40.015, 6); + expect(parseCoordinate('-105.27', 'lon')).toBeCloseTo(-105.27, 6); + expect(parseCoordinate('40°0\'54"N')).toBeCloseTo(40.015, 6); + expect(parseCoordinate('105 16 12 W', 'lon')).toBeCloseTo(-105.27, 6); + expect(parseCoordinate('33.9 S')).toBeCloseTo(-33.9, 6); + expect(parseCoordinate('-0.5')).toBeCloseTo(-0.5, 6); + }); + + test('invalid coordinates are null', () => { + expect(parseCoordinate('')).toBeNull(); + expect(parseCoordinate(null)).toBeNull(); + expect(parseCoordinate('91')).toBeNull(); + expect(parseCoordinate('181', 'lon')).toBeNull(); + expect(parseCoordinate('40 E')).toBeNull(); + expect(parseCoordinate('abc')).toBeNull(); + expect(parseCoordinate('40 75 0')).toBeNull(); + }); + + test('monthly lists round-trip and tolerate blanks', () => { + expect(parseMonthly('1,2,,4')).toEqual([1, 2, null, 4, null, null, null, null, null, null, null, null]); + expect(parseMonthly(null)).toHaveLength(12); + expect(parseMonthly(formatMonthly(GHI_40N))).toEqual(GHI_40N); + }); +}); diff --git a/app/src/components/calculators/solarData.fixtures.js b/app/src/components/calculators/solarData.fixtures.js new file mode 100644 index 00000000..3079df27 --- /dev/null +++ b/app/src/components/calculators/solarData.fixtures.js @@ -0,0 +1,48 @@ +import {cellIndex, COLUMNS, PARAMETERS, ROWS} from "./solarData"; + +// Build a climatology file for tests, the same format scripts/build_solar_climatology.py writes. +// `valueAt(key, month, latitude, longitude)` returns the value for a cell center, or null for no data. +export function encodeClimatology(valueAt) { + const header = 8; + const bytes = new Uint8Array(header + PARAMETERS.length * 12 * ROWS * COLUMNS); + bytes.set([...'WSOL'].map(c => c.charCodeAt(0)), 0); + bytes[4] = 1; + bytes[5] = PARAMETERS.length; + PARAMETERS.forEach(({key, offset, step}, p) => { + for (let month = 0; month < 12; month++) { + for (let row = 0; row < ROWS; row++) { + for (let column = 0; column < COLUMNS; column++) { + const value = valueAt(key, month, row - 89.5, column - 179.5); + bytes[header + cellIndex(p, month, row, column)] = value === null + ? 255 + : Math.max(0, Math.min(254, Math.round((value + offset) / step))); + } + } + } + }); + return bytes.buffer; +} + +// A plain climatology: the same sunny 40°N profile everywhere. +export const SUNNY_GHI = [2.6, 3.5, 4.7, 5.8, 6.5, 7.2, 7.0, 6.2, 5.2, 4.0, 2.9, 2.3]; +export const MILD_TEMPERATURE = [0, 2, 6, 10, 15, 21, 24, 23, 18, 11, 4, 0]; + +// Encoding the 3 MB file is slow, so it is built once and shared; nothing writes to it. +let uniform = null; +export const uniformClimatology = () => { + if (!uniform) { + uniform = encodeClimatology((key, month) => ({ + ghi: SUNNY_GHI[month], + diffuse: SUNNY_GHI[month] * 0.3, + albedo: 0.2, + temperature: MILD_TEMPERATURE[month], + })[key]); + } + return uniform; +}; + +// Make `fetch` return a climatology file. +export function mockClimatologyFetch(buffer = uniformClimatology()) { + global.fetch = jest.fn(() => Promise.resolve({ok: true, arrayBuffer: () => Promise.resolve(buffer)})); + return global.fetch; +} diff --git a/app/src/components/calculators/solarData.js b/app/src/components/calculators/solarData.js new file mode 100644 index 00000000..2c8fda84 --- /dev/null +++ b/app/src/components/calculators/solarData.js @@ -0,0 +1,133 @@ +import React from "react"; + +// Offline solar climatology for the Solar calculator: NASA POWER's 2001-2020 monthly averages on a +// 1° grid, built by scripts/build_solar_climatology.py. The format is described in that script and +// must match PARAMETERS below. The file is only fetched when the calculator opens. + +export const CLIMATOLOGY_URL = `${process.env.PUBLIC_URL || ''}/data/solar-climatology.bin`; +export const CLIMATOLOGY_SOURCE = 'NASA POWER 2001–2020 monthly averages (CC BY 4.0)'; + +const MAGIC = 'WSOL'; +const VERSION = 1; +const HEADER_BYTES = 8; +export const ROWS = 180; +export const COLUMNS = 360; +const MISSING = 255; + +// Stored as round((value + offset) / step), in this order. +export const PARAMETERS = [ + {key: 'ghi', offset: 0, step: 0.04}, + {key: 'diffuse', offset: 0, step: 0.04}, + {key: 'albedo', offset: 0, step: 0.004}, + {key: 'temperature', offset: 70, step: 0.5}, +]; + +const PARAMETER_BYTES = 12 * ROWS * COLUMNS; + +// Validate a downloaded file and return its value grid. +export function decodeClimatology(buffer) { + const bytes = new Uint8Array(buffer); + const magic = String.fromCharCode(...bytes.slice(0, 4)); + if (magic !== MAGIC) { + throw new Error('Not a solar climatology file'); + } + if (bytes[4] !== VERSION || bytes[5] !== PARAMETERS.length) { + throw new Error(`Unsupported solar climatology version ${bytes[4]}`); + } + if (bytes.length !== HEADER_BYTES + PARAMETERS.length * PARAMETER_BYTES) { + throw new Error('Solar climatology file is truncated'); + } + return bytes.subarray(HEADER_BYTES); +} + +// Index of a stored value. +export const cellIndex = (parameter, month, row, column) => + parameter * PARAMETER_BYTES + (month * ROWS + row) * COLUMNS + column; + +/** + * Monthly averages at a location, bilinearly interpolated between the four surrounding cell centers. + * Cells with no data are skipped. Longitude wraps across the antimeridian; latitude is clamped at the + * outermost cell centers. + * + * @returns {{ghi: number[], diffuse: number[], albedo: number[], temperature: number[]}} 12 values each, + * null for a month with no data + */ +export function lookupClimatology(grid, latitude, longitude) { + const y = Math.min(ROWS - 1, Math.max(0, latitude + 90 - 0.5)); + const row0 = Math.floor(y); + const row1 = Math.min(ROWS - 1, row0 + 1); + const fy = y - row0; + + const x = longitude + 180 - 0.5; + const xFloor = Math.floor(x); + const fx = x - xFloor; + const column0 = ((xFloor % COLUMNS) + COLUMNS) % COLUMNS; + const column1 = (column0 + 1) % COLUMNS; + + const corners = [ + [row0, column0, (1 - fy) * (1 - fx)], + [row0, column1, (1 - fy) * fx], + [row1, column0, fy * (1 - fx)], + [row1, column1, fy * fx], + ]; + + const result = {}; + PARAMETERS.forEach(({key, offset, step}, p) => { + result[key] = Array.from({length: 12}, (_, month) => { + let sum = 0, weight = 0; + for (const [row, column, w] of corners) { + const raw = grid[cellIndex(p, month, row, column)]; + if (raw !== MISSING && w > 0) { + sum += (raw * step - offset) * w; + weight += w; + } + } + return weight > 0 ? sum / weight : null; + }); + }); + return result; +} + +let climatologyPromise = null; + +// Fetch and decode the file once per page load. A failed fetch is not cached, so it can be retried. +export function loadClimatology() { + if (!climatologyPromise) { + climatologyPromise = fetch(CLIMATOLOGY_URL) + .then(response => { + if (!response.ok) { + throw new Error(`Could not load solar data (${response.status})`); + } + return response.arrayBuffer(); + }) + .then(decodeClimatology) + .catch(error => { + climatologyPromise = null; + throw error; + }); + } + return climatologyPromise; +} + +// Drop the loaded file, so the next load fetches again. Tests use this to simulate a failed download. +export function forgetClimatology() { + climatologyPromise = null; +} + +// The decoded grid, or null while loading. `error` is set when the file could not be loaded. +export function useClimatology() { + const [grid, setGrid] = React.useState(null); + const [error, setError] = React.useState(null); + + React.useEffect(() => { + let cancelled = false; + loadClimatology() + .then(g => !cancelled && setGrid(g)) + .catch(e => !cancelled && setError(e)); + return () => { + cancelled = true; + }; + }, []); + + return {grid, error, loading: !grid && !error}; +} diff --git a/app/src/components/calculators/solarData.test.js b/app/src/components/calculators/solarData.test.js new file mode 100644 index 00000000..c3ee7e26 --- /dev/null +++ b/app/src/components/calculators/solarData.test.js @@ -0,0 +1,141 @@ +import fs from "fs"; +import path from "path"; +import {decodeClimatology, lookupClimatology} from "./solarData"; +import {estimateProduction} from "./solar"; +import {encodeClimatology, uniformClimatology} from "./solarData.fixtures"; + +describe('decodeClimatology', () => { + test('accepts a valid file', () => { + expect(decodeClimatology(uniformClimatology()).length).toBe(4 * 12 * 180 * 360); + }); + + test('rejects a file that is not ours, a newer version, or truncated', () => { + const good = new Uint8Array(uniformClimatology()); + const wrongMagic = good.slice(); + wrongMagic[0] = 'X'.charCodeAt(0); + expect(() => decodeClimatology(wrongMagic.buffer)).toThrow('Not a solar climatology file'); + + const newer = good.slice(); + newer[4] = 2; + expect(() => decodeClimatology(newer.buffer)).toThrow('Unsupported'); + + expect(() => decodeClimatology(good.slice(0, 1000).buffer)).toThrow('truncated'); + }); +}); + +describe('lookupClimatology', () => { + // A field that varies with longitude, so interpolation is visible: GHI = 3 + lon/100. + const sloped = decodeClimatology(encodeClimatology((key, month, lat, lon) => ({ + ghi: 3 + lon / 100, + diffuse: 1, + albedo: lat > 60 ? 0.8 : 0.2, + temperature: -10 + month, + })[key])); + + test('returns a cell center value exactly', () => { + const values = lookupClimatology(sloped, 40.5, -104.5); + expect(values.ghi[0]).toBeCloseTo(3 - 1.045, 1); + expect(values.temperature.map(v => Math.round(v))).toEqual([-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1]); + expect(values.diffuse[6]).toBeCloseTo(1, 1); + }); + + test('interpolates between cell centers', () => { + // Halfway between the cells centered on -105.5 and -104.5. + const mid = lookupClimatology(sloped, 40.5, -105).ghi[0]; + const west = lookupClimatology(sloped, 40.5, -105.5).ghi[0]; + const east = lookupClimatology(sloped, 40.5, -104.5).ghi[0]; + expect(mid).toBeCloseTo((west + east) / 2, 6); + // Interpolates in latitude too: halfway between albedo 0.2 and 0.8 cells. + expect(lookupClimatology(sloped, 60, 0).albedo[0]).toBeCloseTo(0.5, 2); + }); + + test('wraps across the antimeridian', () => { + // Between the cells centered on 179.5 (ghi 4.795) and -179.5 (ghi 1.205). + expect(lookupClimatology(sloped, 0.5, 180).ghi[0]).toBeCloseTo((4.795 + 1.205) / 2, 1); + expect(lookupClimatology(sloped, 0.5, -180).ghi[0]).toBeCloseTo((4.795 + 1.205) / 2, 1); + }); + + test('clamps at the poles', () => { + expect(lookupClimatology(sloped, 90, 0.5).ghi[0]).toBeCloseTo(3.005, 1); + expect(lookupClimatology(sloped, -90, 0.5).ghi[0]).toBeCloseTo(3.005, 1); + }); + + test('skips cells with no data and returns null when all are missing', () => { + const holes = decodeClimatology(encodeClimatology((key, month, lat, lon) => + (lon > 0 ? null : {ghi: 5, diffuse: 1, albedo: 0.2, temperature: 10}[key]))); + // Straddling lon 0: only the western cells have data. + expect(lookupClimatology(holes, 10.5, 0).ghi[0]).toBeCloseTo(5, 1); + expect(lookupClimatology(holes, 10.5, 90).ghi[0]).toBeNull(); + }); +}); + +// The real data file, when it has been built. The expected values come from NASA POWER's point API +// (2001-2020 climatology), so these confirm the build script and decoder agree end to end. +const REAL_FILE = path.join(__dirname, '../../../public/data/solar-climatology.bin'); +const describeReal = fs.existsSync(REAL_FILE) ? describe : describe.skip; + +describeReal('the built climatology file', () => { + let grid; + beforeAll(() => { + grid = decodeClimatology(fs.readFileSync(REAL_FILE).buffer.slice(0)); + }); + + // NASA POWER's monthly GHI for two cells, copied from its API (2001-2020 climatology). A cell + // center is returned exactly, give or take the file's 0.04 kWh/m²/day step. + const expectCell = (actual, expected) => actual.forEach((v, i) => expect(Math.abs(v - expected[i])) + .toBeLessThanOrEqual(0.021)); + + test('matches NASA POWER at the cell near Denver', () => { + const values = lookupClimatology(grid, 40.5, -104.5); + expectCell(values.ghi, + [2.4089, 3.2969, 4.6241, 5.6225, 6.3463, 7.2588, 7.0428, 6.2326, 5.2351, 3.7486, 2.6717, 2.0921]); + expect(values.diffuse.every((v, i) => v > 0 && v < values.ghi[i])).toBe(true); + expect(values.temperature[6]).toBeGreaterThan(values.temperature[0] + 15); + }); + + test('matches NASA POWER at the cell near Sydney, with the seasons reversed', () => { + const values = lookupClimatology(grid, -33.5, 151.5); + expectCell(values.ghi, + [6.4908, 5.6611, 4.8137, 3.8191, 2.9803, 2.357, 2.772, 3.6725, 4.8242, 5.6698, 6.2074, 6.5594]); + }); + + // PVWatts v8 results for a 1 kW standard-module, open-rack, fixed array with the default 14.08% + // losses, DC output in kWh, fetched from developer.nlr.gov on 2026-10-04. PVWatts uses NREL's NSRDB + // typical-year weather, not NASA POWER, so the two agree only within the data sets' differences. + const PVWATTS = { + Phoenix: { + lat: 33.45, lon: -112.07, tilt: 33, azimuth: 180, + monthly: [139.7, 139.1, 166.3, 175.4, 175.8, 165.2, 152.4, 155.2, 157.6, 157.8, 144.8, 140.3], + }, + Seattle: { + lat: 47.61, lon: -122.33, tilt: 48, azimuth: 180, + monthly: [51.9, 67.3, 89.6, 115.5, 127.1, 116.7, 130.3, 141.0, 111.9, 82.0, 48.6, 42.8], + }, + Anchorage: { + lat: 61.22, lon: -149.90, tilt: 61, azimuth: 180, + monthly: [19.6, 45.0, 119.2, 117.4, 124.2, 118.0, 109.4, 109.4, 87.8, 52.0, 23.3, 12.6], + }, + }; + + test.each(Object.entries(PVWATTS))('agrees with PVWatts in %s', (name, site) => { + const sunlight = lookupClimatology(grid, site.lat, site.lon); + const result = estimateProduction({ + latitude: site.lat, kwp: 1, tilt: site.tilt, azimuth: site.azimuth, ...sunlight, + }); + const annual = site.monthly.reduce((a, b) => a + b, 0); + expect(result.annual / annual).toBeGreaterThan(0.9); + expect(result.annual / annual).toBeLessThan(1.1); + // Dark months are small and the two data sets disagree most there; check the months that matter. + const meanMonth = annual / 12; + result.months.forEach((m, i) => { + if (site.monthly[i] >= 0.6 * meanMonth) { + expect(m.energy / site.monthly[i]).toBeGreaterThan(0.75); + expect(m.energy / site.monthly[i]).toBeLessThan(1.25); + } + }); + }); + + test('has polar night', () => { + expect(lookupClimatology(grid, 80, 15).ghi[11]).toBeLessThan(0.05); + }); +}); diff --git a/scripts/build_solar_climatology.py b/scripts/build_solar_climatology.py new file mode 100644 index 00000000..f4673c56 --- /dev/null +++ b/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('