diff --git a/.changeset/pv-array-draw-on-map.md b/.changeset/pv-array-draw-on-map.md new file mode 100644 index 00000000..c01ea4fa --- /dev/null +++ b/.changeset/pv-array-draw-on-map.md @@ -0,0 +1,11 @@ +--- +"ftw": minor +--- + +PV arrays can be drawn on the map instead of measured by hand. The Weather tab gains a rectangle tool over the existing Leaflet map: draw over your panels and it becomes a `weather.pv_arrays` entry, with area from the shape and azimuth from the way the shape is turned. + +The rectangle is drawn at an angle rather than square to north, because the angle is the point — the long edge follows the ridge, and the face is perpendicular to it. Two directions are perpendicular to a ridge and a flat outline genuinely does not say which, so the equatorward one is offered as a default with a one-click flip, never as a measurement. + +Tilt is the one number an overhead outline cannot contain, so it is typed once before drawing — and it is also what turns the outline into panel area. What you trace on a map is the horizontal projection of a sloped rectangle, so a 35° roof carries about 22 % more panel than its outline suggests; ignoring that would quietly under-size every drawn array. Capacity uses the same packing factor and module density as the Lantmäteriet roof model, so a drawn array and a derived one are comparable. + +Drawing is progressive enhancement: [Terra Draw](https://github.com/JamesLMilner/terra-draw) and its Leaflet adapter (both MIT) are lazy-loaded from CDN with SRI hashes only when the tool is first used, and if either fails to load the numeric editor is untouched and says so. The editor stays the final word for drawn and typed arrays alike. diff --git a/web/components/pv-array-geometry.js b/web/components/pv-array-geometry.js new file mode 100644 index 00000000..aa190d5f --- /dev/null +++ b/web/components/pv-array-geometry.js @@ -0,0 +1,208 @@ +// Turning a rectangle drawn on the map into a PV array. +// +// Drawing supplies two of the three numbers an array needs. Area comes from +// the shape; azimuth from how the shape is turned. Tilt cannot be seen from +// directly above at all, so it stays typed — and it is also what converts the +// drawn outline into real panel area, because what you trace on a map is the +// *horizontal projection* of a sloped rectangle, not the rectangle itself. +// +// Everything here is pure so it can be tested without a browser or a map. + +export const DEFAULT_PACKING_FACTOR = 0.7; +export const DEFAULT_MODULE_W_PER_M2 = 200; +export const DEFAULT_TILT_DEG = 35; + +// IUGG mean Earth radius. At the scale of one roof the radius matters far +// less than the flat-Earth approximation below, which is exact enough for a +// 20 m rectangle and meaningless across a county. +const EARTH_RADIUS_M = 6371008.8; +const DEG = Math.PI / 180; + +// Beyond this a roof is a wall: cos(tilt) approaches zero and the plan-area +// division runs away. A wall has no horizontal projection to trace anyway, so +// clamping here bounds the arithmetic instead of returning Infinity. +const MAX_TILT_FOR_PROJECTION_DEG = 89; + +function stripClosingVertex(ring) { + if (ring.length > 1) { + const first = ring[0]; + const last = ring[ring.length - 1]; + if (first[0] === last[0] && first[1] === last[1]) return ring.slice(0, -1); + } + return ring.slice(); +} + +/** + * Project a WGS84 ring ([[lon, lat], …]) to metres about its own centroid. + * + * A local tangent plane, not a real projection: over one building the error + * is well under the precision anyone draws with, and it avoids carrying a + * projection library into the settings page. + */ +export function toLocalMetres(ring) { + const pts = stripClosingVertex(ring || []); + if (pts.length === 0) return []; + let lon0 = 0; + let lat0 = 0; + for (const [lon, lat] of pts) { + lon0 += lon; + lat0 += lat; + } + lon0 /= pts.length; + lat0 /= pts.length; + const mPerDegLat = EARTH_RADIUS_M * DEG; + const mPerDegLon = mPerDegLat * Math.cos(lat0 * DEG); + return pts.map(([lon, lat]) => [(lon - lon0) * mPerDegLon, (lat - lat0) * mPerDegLat]); +} + +/** Area of the drawn outline in m², as seen from above. */ +export function planAreaM2(ring) { + const p = toLocalMetres(ring); + if (p.length < 3) return 0; + let twiceArea = 0; + for (let i = 0; i < p.length; i++) { + const [x1, y1] = p[i]; + const [x2, y2] = p[(i + 1) % p.length]; + twiceArea += x1 * y2 - x2 * y1; + } + return Math.abs(twiceArea) / 2; +} + +/** Compass bearing of a local vector, 0 = north, 90 = east. */ +function bearingDeg(dx, dy) { + return (((Math.atan2(dx, dy) / DEG) % 360) + 360) % 360; +} + +/** + * Direction of the ring's longest edge, as a line in [0, 180). + * + * For a panel rectangle that edge runs along the ridge, which is a line and + * not an arrow — calling it "north" rather than "south" would be a + * distinction the drawing does not contain. + */ +export function ridgeAzimuthDeg(ring) { + const p = toLocalMetres(ring); + if (p.length < 2) return null; + let longest = 0; + let bx = 0; + let by = 0; + for (let i = 0; i < p.length; i++) { + const [x1, y1] = p[i]; + const [x2, y2] = p[(i + 1) % p.length]; + const dx = x2 - x1; + const dy = y2 - y1; + const len = Math.hypot(dx, dy); + if (len > longest) { + longest = len; + bx = dx; + by = dy; + } + } + if (longest <= 0) return null; + return bearingDeg(bx, by) % 180; +} + +/** Shortest angle between two compass bearings, in degrees. */ +export function angularDistanceDeg(a, b) { + const d = Math.abs((((a - b) % 360) + 360) % 360); + return d > 180 ? 360 - d : d; +} + +/** Turn an azimuth to face the opposite way. */ +export function flipAzimuthDeg(azimuthDeg) { + return ((((azimuthDeg + 180) % 360) + 360) % 360); +} + +/** + * The two directions the face could point: perpendicular to the ridge, either + * side of it. A flat outline genuinely does not say which. + */ +export function faceAzimuthCandidates(ring) { + const ridge = ridgeAzimuthDeg(ring); + if (ridge === null) return []; + return [(ridge + 90) % 360, (ridge + 270) % 360]; +} + +/** + * The candidate a panel is more likely to use: the equatorward one. + * + * This is a default, not a measurement. Both perpendiculars fit the drawing + * equally well, so the UI offers a flip rather than pretending to know. + */ +export function preferredAzimuthDeg(ring, latitudeDeg) { + const candidates = faceAzimuthCandidates(ring); + if (candidates.length === 0) return null; + const target = (latitudeDeg || 0) >= 0 ? 180 : 0; + const [a, b] = candidates; + return angularDistanceDeg(a, target) <= angularDistanceDeg(b, target) ? a : b; +} + +/** + * Real panel area from the traced outline. + * + * A sloped rectangle of area A casts a shadow of A·cos(tilt) on the map, so + * recovering it divides that back out. A 35° roof carries about 22 % more + * panel than its outline suggests, which is the difference between a + * believable kWp and a quietly low one. + */ +export function slopeAreaM2(planArea, tiltDeg) { + const tilt = Math.min(Math.max(tiltDeg || 0, 0), MAX_TILT_FOR_PROJECTION_DEG); + return planArea / Math.cos(tilt * DEG); +} + +/** Installable DC capacity for a roof area, matching the roof model's basis. */ +export function kwpFromSlopeArea(areaM2, packingFactor, moduleWPerM2) { + const packing = packingFactor == null ? DEFAULT_PACKING_FACTOR : packingFactor; + const wPerM2 = moduleWPerM2 == null ? DEFAULT_MODULE_W_PER_M2 : moduleWPerM2; + return (areaM2 * packing * wPerM2) / 1000; +} + +/** Human-readable face name, mirroring the roof model's naming. */ +export function compassName(azimuthDeg, tiltDeg) { + if (tiltDeg < 5) return "Roof flat"; + const points = [ + [0, "north"], [45, "north-east"], [90, "east"], [135, "south-east"], + [180, "south"], [225, "south-west"], [270, "west"], [315, "north-west"], + [360, "north"], + ]; + let best = points[0]; + for (const p of points) { + if (Math.abs(p[0] - azimuthDeg) < Math.abs(best[0] - azimuthDeg)) best = p; + } + return `Roof ${best[1]}`; +} + +function round(value, places) { + const factor = 10 ** places; + return Math.round(value * factor) / factor; +} + +/** + * Everything a drawn rectangle says about one array. + * + * Returns the config-shaped entry separately from the measurements, so only + * the four fields weather.pv_arrays actually defines are ever written back. + */ +export function arrayFromRing(ring, options) { + const opts = options || {}; + const plan = planAreaM2(ring); + if (!(plan > 0)) return null; + const tiltDeg = opts.tiltDeg == null ? DEFAULT_TILT_DEG : opts.tiltDeg; + const candidates = faceAzimuthCandidates(ring); + const azimuth = opts.azimuthDeg == null + ? preferredAzimuthDeg(ring, opts.latitude) + : opts.azimuthDeg; + const azimuthDeg = azimuth == null ? 180 : Math.round(azimuth); + const slope = slopeAreaM2(plan, tiltDeg); + return { + array: { + name: opts.name || compassName(azimuthDeg, tiltDeg), + kwp: round(kwpFromSlopeArea(slope, opts.packingFactor, opts.moduleWPerM2), 2), + tilt_deg: tiltDeg, + azimuth_deg: azimuthDeg, + }, + planAreaM2: round(plan, 1), + slopeAreaM2: round(slope, 1), + azimuthCandidates: candidates.map((c) => Math.round(c)), + }; +} diff --git a/web/components/pv-array-geometry.test.mjs b/web/components/pv-array-geometry.test.mjs new file mode 100644 index 00000000..0dd99679 --- /dev/null +++ b/web/components/pv-array-geometry.test.mjs @@ -0,0 +1,203 @@ +// node --test web/components/pv-array-geometry.test.mjs + +import assert from "node:assert/strict"; +import { describe, it } from "node:test"; + +import { + DEFAULT_MODULE_W_PER_M2, + DEFAULT_PACKING_FACTOR, + angularDistanceDeg, + arrayFromRing, + compassName, + faceAzimuthCandidates, + flipAzimuthDeg, + kwpFromSlopeArea, + planAreaM2, + preferredAzimuthDeg, + ridgeAzimuthDeg, + slopeAreaM2, +} from "./pv-array-geometry.js"; + +// Fixtures are built with the standard ellipsoidal metres-per-degree series, +// deliberately *not* with the module's own spherical projection — a shape +// round-tripped through the code under test would agree with itself and prove +// nothing. The two disagree by roughly half a percent in area at this +// latitude, which is the known bias of a spherical Earth against WGS84 and is +// far below the precision anyone draws a roof with. +const STOCKHOLM = { lat: 59.3293, lon: 18.0686 }; + +function metresPerDegree(latDeg) { + const p = (latDeg * Math.PI) / 180; + return { + lat: 111132.92 - 559.82 * Math.cos(2 * p) + 1.175 * Math.cos(4 * p) + - 0.0023 * Math.cos(6 * p), + lon: 111412.84 * Math.cos(p) - 93.5 * Math.cos(3 * p) + 0.118 * Math.cos(5 * p), + }; +} + +/** Build a WGS84 ring from local east/north offsets in metres. */ +function ringFromMetres(offsets, origin) { + const m = metresPerDegree(origin.lat); + return offsets.map(([x, y]) => [origin.lon + x / m.lon, origin.lat + y / m.lat]); +} + +/** Rotate local offsets counter-clockwise in the east/north plane. */ +function rotate(offsets, degrees) { + const r = (degrees * Math.PI) / 180; + const c = Math.cos(r); + const s = Math.sin(r); + return offsets.map(([x, y]) => [x * c - y * s, x * s + y * c]); +} + +// 10 m along the ridge (east-west) by 6 m down the slope. +const RECT_10x6 = [[-5, -3], [5, -3], [5, 3], [-5, 3]]; + +describe("plan area", () => { + it("recovers the drawn size in square metres", () => { + const ring = ringFromMetres(RECT_10x6, STOCKHOLM); + const area = planAreaM2(ring); + assert.ok(Math.abs(area - 60) < 0.6, `area ${area} should be ~60 m²`); + }); + + it("does not care whether the ring repeats its first point", () => { + const ring = ringFromMetres(RECT_10x6, STOCKHOLM); + const closed = [...ring, ring[0]]; + assert.ok(Math.abs(planAreaM2(ring) - planAreaM2(closed)) < 1e-9); + }); + + it("is unsigned, so winding order cannot produce a negative roof", () => { + const ring = ringFromMetres(RECT_10x6, STOCKHOLM); + assert.ok(Math.abs(planAreaM2(ring) - planAreaM2([...ring].reverse())) < 1e-9); + }); + + it("treats a shape with no area as no array", () => { + assert.equal(planAreaM2([]), 0); + assert.equal(planAreaM2([[18, 59], [18.001, 59]]), 0); + assert.equal(arrayFromRing([[18, 59], [18.001, 59]], {}), null); + }); +}); + +describe("orientation", () => { + it("reads the ridge from the longest edge", () => { + const ring = ringFromMetres(RECT_10x6, STOCKHOLM); + // The 10 m edges run east-west: a ridge bearing of 90°. + assert.ok(Math.abs(ridgeAzimuthDeg(ring) - 90) < 0.5); + }); + + it("offers both faces the outline permits, and no others", () => { + const ring = ringFromMetres(RECT_10x6, STOCKHOLM); + const [a, b] = faceAzimuthCandidates(ring).map(Math.round); + assert.deepEqual([a, b].sort((x, y) => x - y), [0, 180]); + }); + + it("defaults to the equatorward face, per hemisphere", () => { + const north = ringFromMetres(RECT_10x6, STOCKHOLM); + assert.ok(Math.abs(preferredAzimuthDeg(north, STOCKHOLM.lat) - 180) < 0.5); + + const south = ringFromMetres(RECT_10x6, { lat: -33.87, lon: 151.21 }); + const picked = preferredAzimuthDeg(south, -33.87); + assert.ok(angularDistanceDeg(picked, 0) < 0.5, `expected ~0°, got ${picked}`); + }); + + it("follows the rectangle round as it turns", () => { + // Turning the shape 30° counter-clockwise swings the ridge from 90° to + // 60°, so the faces move with it: 150° and 330°, and south-ish wins. + const ring = ringFromMetres(rotate(RECT_10x6, 30), STOCKHOLM); + assert.ok(Math.abs(ridgeAzimuthDeg(ring) - 60) < 0.5); + assert.ok(Math.abs(preferredAzimuthDeg(ring, STOCKHOLM.lat) - 150) < 0.5); + }); + + it("keeps the ridge a line rather than an arrow", () => { + // Drawing the same rectangle the other way round is the same roof. + const ring = ringFromMetres(RECT_10x6, STOCKHOLM); + const reversed = ringFromMetres([...RECT_10x6].reverse(), STOCKHOLM); + assert.ok(Math.abs(ridgeAzimuthDeg(ring) - ridgeAzimuthDeg(reversed)) < 0.5); + }); + + it("flips to the opposite face", () => { + assert.equal(flipAzimuthDeg(180), 0); + assert.equal(flipAzimuthDeg(0), 180); + assert.equal(flipAzimuthDeg(270), 90); + assert.equal(flipAzimuthDeg(350), 170); + }); + + it("measures the shorter way round the compass", () => { + assert.equal(angularDistanceDeg(350, 10), 20); + assert.equal(angularDistanceDeg(10, 350), 20); + assert.equal(angularDistanceDeg(0, 180), 180); + }); +}); + +describe("tilt turns an outline into panel area", () => { + it("leaves a flat roof alone", () => { + assert.ok(Math.abs(slopeAreaM2(60, 0) - 60) < 1e-9); + }); + + it("recovers the area hidden by the slope", () => { + // cos 60° = 0.5, so a 60 m² shadow is cast by 120 m² of roof. + assert.ok(Math.abs(slopeAreaM2(60, 60) - 120) < 1e-9); + // A 35° roof carries ~22 % more panel than its outline suggests. + assert.ok(Math.abs(slopeAreaM2(60, 35) - 73.24) < 0.05); + }); + + it("stays finite at a wall, where there is no outline to trace", () => { + assert.ok(Number.isFinite(slopeAreaM2(60, 90))); + }); + + it("means a steeper roof is a bigger array for the same drawing", () => { + const ring = ringFromMetres(RECT_10x6, STOCKHOLM); + const flat = arrayFromRing(ring, { latitude: STOCKHOLM.lat, tiltDeg: 0 }); + const steep = arrayFromRing(ring, { latitude: STOCKHOLM.lat, tiltDeg: 45 }); + assert.ok(steep.array.kwp > flat.array.kwp, + `${steep.array.kwp} should exceed ${flat.array.kwp}`); + }); +}); + +describe("capacity", () => { + it("uses the same basis as the roof model", () => { + // 60 m² × 0.70 packing × 200 W/m² = 8.4 kWp. + assert.ok(Math.abs(kwpFromSlopeArea(60) - 8.4) < 1e-9); + assert.equal(DEFAULT_PACKING_FACTOR, 0.7); + assert.equal(DEFAULT_MODULE_W_PER_M2, 200); + }); + + it("honours an overridden packing factor", () => { + assert.ok(Math.abs(kwpFromSlopeArea(60, 0.5, 200) - 6) < 1e-9); + }); +}); + +describe("the entry written back to config", () => { + it("carries only the four fields weather.pv_arrays defines", () => { + const ring = ringFromMetres(RECT_10x6, STOCKHOLM); + const out = arrayFromRing(ring, { latitude: STOCKHOLM.lat }); + assert.deepEqual( + Object.keys(out.array).sort(), + ["azimuth_deg", "kwp", "name", "tilt_deg"], + ); + }); + + it("describes a south-facing 35° roof from the drawing alone", () => { + const ring = ringFromMetres(RECT_10x6, STOCKHOLM); + const out = arrayFromRing(ring, { latitude: STOCKHOLM.lat }); + assert.equal(out.array.azimuth_deg, 180); + assert.equal(out.array.tilt_deg, 35); + assert.equal(out.array.name, "Roof south"); + assert.ok(Math.abs(out.planAreaM2 - 60) < 0.6); + assert.ok(Math.abs(out.slopeAreaM2 - 73.2) < 0.6); + assert.ok(Math.abs(out.array.kwp - 10.25) < 0.1); + assert.deepEqual(out.azimuthCandidates.slice().sort((a, b) => a - b), [0, 180]); + }); + + it("lets an explicit azimuth override the guess", () => { + const ring = ringFromMetres(RECT_10x6, STOCKHOLM); + const out = arrayFromRing(ring, { latitude: STOCKHOLM.lat, azimuthDeg: 0 }); + assert.equal(out.array.azimuth_deg, 0); + assert.equal(out.array.name, "Roof north"); + }); + + it("names a flat roof for what it is", () => { + assert.equal(compassName(180, 0), "Roof flat"); + assert.equal(compassName(90, 35), "Roof east"); + assert.equal(compassName(225, 35), "Roof south-west"); + }); +}); diff --git a/web/index.html b/web/index.html index ec2a315e..f051a44f 100644 --- a/web/index.html +++ b/web/index.html @@ -937,7 +937,7 @@
' + + 'Drawing gives the size and the direction; tilt is the one thing an overhead ' + + 'outline cannot show, so set it above before you draw — a 35° roof holds about ' + + '22 % more panel than its outline suggests. ' + 'Tilt: 0° = flat roof, 35° = typical pitched roof, 90° = wall. Azimuth: 0 = N, 90 = E, 180 = S, 270 = W.' + '
' + ''; @@ -198,6 +350,21 @@ ctx.config.weather.pv_arrays.push({ name: "", kwp: 0, tilt_deg: 35, azimuth_deg: 180 }); renderPVArrays(ctx); }); + var drawBtn = document.getElementById("pv-array-draw"); + if (drawBtn) drawBtn.addEventListener("click", function () { startArrayDrawing(ctx); }); + var doneBtn = document.getElementById("pv-array-draw-done"); + if (doneBtn) doneBtn.addEventListener("click", stopArrayDrawing); + var status = document.getElementById("pv-draw-status"); + if (status) status.addEventListener("click", function (e) { + if (!e.target || e.target.id !== "pv-draw-flip" || !lastDrawnArray) return; + // Flip the object, not an index: removing another row above it would + // otherwise silently turn a different roof around. + lastDrawnArray.azimuth_deg = drawGeometry.flipAzimuthDeg(lastDrawnArray.azimuth_deg); + lastDrawnArray.name = drawGeometry.compassName( + lastDrawnArray.azimuth_deg, lastDrawnArray.tilt_deg); + renderPVArrays(ctx); + drawStatus("Now facing " + lastDrawnArray.azimuth_deg + "°."); + }); }, }; })();