fuzzyworld

tower / cascadia / files / tools / cascadia.mjs

Branch cascadia · default branch · all branches

cascadia.mjs

Blame: the commit that last changed each line · file · history

1affced2 2026-09-08 Pavel Safronov1/**
2 * Cascadia ingest — the "real world, real size" mode. `npm run cascadia`
3 *
4 * Tower's inner surface is 540 km along the axis by 314.159 km around. Measured in
5 * UTM 10N, the ground from the Pacific surf at Kalaloch east to the Columbia River
6 * at Entiat is 312.8 km — 1.4 km short of the circumference. So western Washington
7 * closes the ring on real terrain: walk east from the beach, over the Olympics,
8 * across Puget Sound, over the Cascades, down the east flank to the Columbia, and
9 * you arrive back at the beach. Nothing is invented and nothing is rescaled. This
10 * is Tower at 1:1, and the point directly over Seattle is a place you can drive to.
11 *
12 * That is the whole reason this mode exists: it is a RULER, not canon. It contains
13 * no Tower geography at all. See the decision record before wiring it to anything.
14 *
15 * Writes, alongside the authored rasters that `npm run terrain` produces:
16 * static/cascadia-height.png 16-bit greyscale, HEIGHT_MIN_M..HEIGHT_MAX_M
17 * static/cascadia-color.png 8-bit RGB, hypsometric tint + hillshade
18 *
19 * Elevation comes from the AWS open terrain tiles (terrarium encoding), which are
20 * SRTM/3DEP over land and GEBCO below water, and which — unlike 3DEP alone — carry
21 * both bathymetry and the Canadian ground north of 49°N in one consistent grid.
22 */
23import { readFileSync, writeFileSync, mkdirSync, existsSync } from 'node:fs'
24import { encodeGray16, encodeRGB, decodeRGB } from './png.mjs'
25
26// ------------------------------------------------------------------ the frame
27/**
28 * Tower has no north. The compass is sunward/deepward/spinward/antispinward
29 * (content/geography.json), so rotating a real region to fit the cylinder is not a
30 * distortion — there is nothing to distort against. ROT_DEG turns the sampling
31 * frame that many degrees west of north, which is how the strip is made to follow
32 * the Cascadia trend rather than a meridian. Kept at 0 until measured against real
33 * elevations; see the seam report this tool prints.
34 */
35const ROT_DEG = Number(process.env.ROT ?? 0)
36/** A real place pinned to a grid pixel. ROT_DEG rotates the frame about this. */
37const ANCHOR = { name: 'Seattle', lat: 47.6062, lon: -122.3321, x: 570, y: 582 }
38/**
39 * Where the ring is cut, as metres spinward (east) of the Kalaloch shoreline. The
40 * cut is the one place the world cannot be continuous, so it is worth choosing
41 * rather than inheriting: sliding it east lifts BOTH edges off the water and onto
42 * ground, which is what actually shrinks the step. See the seam report below.
43 */
44const CUT_SHIFT_M = Number(process.env.CUT ?? 0)
45/**
46 * Zero, after measuring. Sliding the cut 36 km east more than halves the step
47 * through the settled half (195-378 m against 369-844 m), but it selects a
48 * different 314 km of ground: the outer Pacific coast leaves the world, and with
49 * it the fact this mode exists for — that 124.42W and 120.24W are ONE point,
50 * 100 km straight over Seattle, reachable by walking either way round. Rotating
51 * the frame was tried too and does not help at any angle from 0 to 20 degrees:
52 * the mismatch is vertical, not alignment. So the seam stays where it is, an
53 * 800 m scarp running the length of the world exactly where the wet side meets
54 * the dry one, which is the truest thing about this coast.
55 */
56
57const GRID_W = 2048 // x is z, the axis. Sun cap at x = 0.
58const GRID_H = 1192 // y is theta, and it wraps.
59const LENGTH_M = 540_000
60const CIRCUMFERENCE_M = 2 * Math.PI * 50_000
61const MPP_X = LENGTH_M / GRID_W
62const MPP_Y = CIRCUMFERENCE_M / GRID_H
63
64/**
65 * The working range of cascadia-height.png.
66 *
67 * These are not chosen, they are BORROWED: they must equal RASTER_MIN_HEIGHT_M and
68 * RASTER_MAX_HEIGHT_M in src/views/raster.ts exactly. The interior view bakes that
69 * range into its GLSL as literals, so a mode with its own window would need a
70 * second shader; sharing canon's makes this mode a swap of two PNG filenames and
71 * nothing else. It costs nothing to share, because Cascadia happens to fit: the
72 * highest land here is Rainier at 4,352 m, under the 4,400 m ceiling.
73 *
74 * Deep ocean is clamped rather than carried. The abyssal plain is 2.6 km down and
75 * under water, so flattening it is invisible, while spending a third of the range
76 * on it would coarsen every hillside in the world to pay for it.
77 */
78const HEIGHT_MIN_M = -560
79const HEIGHT_MAX_M = 4400
80const SNOWLINE_M = 1900
81
82const ZOOM = 9 // ~206 m/px at 47.6°N, against the grid's 264 m/px
83const TILE = 256
84const CACHE = '.cache/terrain-tiles'
85const TILE_URL = (z, x, y) =>
86 `https://s3.amazonaws.com/elevation-tiles-prod/terrarium/${z}/${x}/${y}.png`
87
88// ------------------------------------------------------------------- progress
89const T0 = Date.now()
90const secs = (from) => `${((Date.now() - from) / 1000).toFixed(1)}s`
91const mb = (n) => `${(n / 1048576).toFixed(1)} MB`
92function line(indent, label, right = '') {
93 const pad = Math.max(1, 58 - indent * 2 - label.length - right.length)
94 console.log(' '.repeat(indent * 2) + label + ' '.repeat(pad) + right)
95}
96function step(indent, label, fn, extra = () => '') {
97 const t = Date.now()
98 const out = fn()
99 line(indent, label, `${extra(out)}${secs(t)}`.trim())
100 return out
101}
102
103// --------------------------------------------------------- WGS84 / UTM zone 10N
104// Snyder's series, good to a centimetre — far past what a 264 m grid can notice.
105const A = 6378137, F = 1 / 298.257223563, K0 = 0.9996, LON0 = (-123 * Math.PI) / 180
106const E2 = F * (2 - F), EP2 = E2 / (1 - E2)
107const M0 = 1 - E2 / 4 - (3 * E2 ** 2) / 64 - (5 * E2 ** 3) / 256
108
109function toUTM(latDeg, lonDeg) {
110 const p = (latDeg * Math.PI) / 180, l = (lonDeg * Math.PI) / 180
111 const sp = Math.sin(p), cp = Math.cos(p), tp = Math.tan(p)
112 const N = A / Math.sqrt(1 - E2 * sp * sp)
113 const T = tp * tp, C = EP2 * cp * cp, a = cp * (l - LON0)
114 const M = A * (M0 * p
115 - (3 * E2 / 8 + 3 * E2 ** 2 / 32 + 45 * E2 ** 3 / 1024) * Math.sin(2 * p)
116 + (15 * E2 ** 2 / 256 + 45 * E2 ** 3 / 1024) * Math.sin(4 * p)
117 - (35 * E2 ** 3 / 3072) * Math.sin(6 * p))
118 return [
119 K0 * N * (a + (1 - T + C) * a ** 3 / 6
120 + (5 - 18 * T + T * T + 72 * C - 58 * EP2) * a ** 5 / 120) + 500000,
121 K0 * (M + N * tp * (a * a / 2 + (5 - T + 9 * C + 4 * C * C) * a ** 4 / 24
122 + (61 - 58 * T + T * T + 600 * C - 330 * EP2) * a ** 6 / 720)),
123 ]
124}
125
126function fromUTM(e, n) {
127 const mu = n / K0 / (A * M0)
128 const e1 = (1 - Math.sqrt(1 - E2)) / (1 + Math.sqrt(1 - E2))
129 const p1 = mu + (3 * e1 / 2 - 27 * e1 ** 3 / 32) * Math.sin(2 * mu)
130 + (21 * e1 ** 2 / 16 - 55 * e1 ** 4 / 32) * Math.sin(4 * mu)
131 + (151 * e1 ** 3 / 96) * Math.sin(6 * mu)
132 + (1097 * e1 ** 4 / 512) * Math.sin(8 * mu)
133 const sp = Math.sin(p1), cp = Math.cos(p1), tp = Math.tan(p1)
134 const C1 = EP2 * cp * cp, T1 = tp * tp
135 const N1 = A / Math.sqrt(1 - E2 * sp * sp)
136 const R1 = (A * (1 - E2)) / (1 - E2 * sp * sp) ** 1.5
137 const D = (e - 500000) / (N1 * K0)
138 const lat = p1 - (N1 * tp / R1) * (D * D / 2
139 - (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * EP2) * D ** 4 / 24
140 + (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * EP2 - 3 * C1 * C1) * D ** 6 / 720)
141 const lon = LON0 + (D - (1 + 2 * T1 + C1) * D ** 3 / 6
142 + (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * EP2 + 24 * T1 * T1) * D ** 5 / 120) / cp
143 return [(lat * 180) / Math.PI, (lon * 180) / Math.PI]
144}
145
146// ------------------------------------------------------------- grid <-> ground
147const BETA = (ROT_DEG * Math.PI) / 180
148// +z (deepward) points BETA west of north; +theta (spinward) is 90° clockwise of it.
149const ZE = -Math.sin(BETA), ZN = Math.cos(BETA)
150const YE = Math.cos(BETA), YN = Math.sin(BETA)
151const [AE, AN] = toUTM(ANCHOR.lat, ANCHOR.lon)
152// Grid origin: back off from the anchor along both axes to reach pixel (0,0).
153const OE = AE - (ANCHOR.x + 0.5) * MPP_X * ZE - ((ANCHOR.y + 0.5) * MPP_Y - CUT_SHIFT_M) * YE
154const ON = AN - (ANCHOR.x + 0.5) * MPP_X * ZN - ((ANCHOR.y + 0.5) * MPP_Y - CUT_SHIFT_M) * YN
155
156/** Grid pixel centre -> UTM 10N easting/northing. */
157const gridToUTM = (x, y) => {
158 const dz = (x + 0.5) * MPP_X, dy = (y + 0.5) * MPP_Y
159 return [OE + dz * ZE + dy * YE, ON + dz * ZN + dy * YN]
160}
161
162// ----------------------------------------------------------------- the tiles
163/**
164 * Web Mercator global pixel coordinates at ZOOM. Terrain tiles are square in
165 * Mercator, so the ground sample distance varies with latitude; at ZOOM 9 it runs
166 * 206 m at Seattle to 192 m at the far cap, comfortably finer than the 264 m grid.
167 */
168const lonToGX = (lon) => ((lon + 180) / 360) * TILE * 2 ** ZOOM
169const latToGY = (lat) => {
170 const s = Math.sin((lat * Math.PI) / 180)
171 return (0.5 - Math.log((1 + s) / (1 - s)) / (4 * Math.PI)) * TILE * 2 ** ZOOM
172}
173
174async function fetchTile(tx, ty) {
175 const path = `${CACHE}/${ZOOM}-${tx}-${ty}.png`
176 if (existsSync(path)) return readFileSync(path)
177 for (let attempt = 1; ; attempt++) {
178 try {
179 const res = await fetch(TILE_URL(ZOOM, tx, ty))
180 if (!res.ok) throw new Error(`HTTP ${res.status}`)
181 const buf = Buffer.from(await res.arrayBuffer())
182 writeFileSync(path, buf)
183 return buf
184 } catch (err) {
185 if (attempt >= 4) throw new Error(`tile ${ZOOM}/${tx}/${ty}: ${err.message}`)
186 await new Promise((r) => setTimeout(r, 400 * attempt))
187 }
188 }
189}
190
191/** Fetch the tiles covering the strip and assemble one elevation field over them. */
192async function loadElevation() {
193 // Bounding box of the (possibly rotated) strip, walked around its edge.
194 let latMin = 90, latMax = -90, lonMin = 180, lonMax = -180
195 for (let i = 0; i <= 64; i++) {
196 for (const [x, y] of [
197 [(i / 64) * GRID_W, 0], [(i / 64) * GRID_W, GRID_H],
198 [0, (i / 64) * GRID_H], [GRID_W, (i / 64) * GRID_H],
199 ]) {
200 const [lat, lon] = fromUTM(...gridToUTM(x, y))
201 latMin = Math.min(latMin, lat); latMax = Math.max(latMax, lat)
202 lonMin = Math.min(lonMin, lon); lonMax = Math.max(lonMax, lon)
203 }
204 }
205 const tx0 = Math.floor(lonToGX(lonMin) / TILE), tx1 = Math.floor(lonToGX(lonMax) / TILE)
206 const ty0 = Math.floor(latToGY(latMax) / TILE), ty1 = Math.floor(latToGY(latMin) / TILE)
207 const cols = tx1 - tx0 + 1, rows = ty1 - ty0 + 1
208 line(1, 'extent',
209 `${latMin.toFixed(2)}..${latMax.toFixed(2)}N ` +
210 `${(-lonMax).toFixed(2)}..${(-lonMin).toFixed(2)}W`)
211
212 const W = cols * TILE, H = rows * TILE
213 const field = new Float32Array(W * H)
214 const jobs = []
215 for (let ty = ty0; ty <= ty1; ty++) for (let tx = tx0; tx <= tx1; tx++) jobs.push([tx, ty])
216
217 let done = 0, fetched = 0
218 const t = Date.now()
219 // Eight at a time: enough to saturate the link, few enough to stay polite.
220 const workers = Array.from({ length: 8 }, async () => {
221 for (;;) {
222 const job = jobs.shift()
223 if (!job) return
224 const [tx, ty] = job
225 if (!existsSync(`${CACHE}/${ZOOM}-${tx}-${ty}.png`)) fetched++
226 const { data, channels } = decodeRGB(await fetchTile(tx, ty))
227 const ox = (tx - tx0) * TILE, oy = (ty - ty0) * TILE
228 for (let py = 0; py < TILE; py++) {
229 for (let px = 0; px < TILE; px++) {
230 const s = (py * TILE + px) * channels
231 // terrarium: (R * 256 + G + B / 256) - 32768 metres
232 field[(oy + py) * W + ox + px] =
233 data[s] * 256 + data[s + 1] + data[s + 2] / 256 - 32768
234 }
235 }
236 if (++done % 20 === 0) process.stdout.write(`\r tiles ${done}/${cols * rows}`)
237 }
238 })
239 await Promise.all(workers)
240 process.stdout.write('\r' + ' '.repeat(30) + '\r')
241 line(1, 'tiles', `${cols}x${rows} at z${ZOOM}, ${fetched} fetched, ${secs(t)}`)
242 return { field, W, H, gx0: tx0 * TILE, gy0: ty0 * TILE }
243}
244
245/** Bilinear sample of the assembled field, in degrees. */
246function sampler(src) {
247 const { field, W, H, gx0, gy0 } = src
248 return (lat, lon) => {
249 const fx = lonToGX(lon) - gx0, fy = latToGY(lat) - gy0
250 const x0 = Math.max(0, Math.min(W - 2, Math.floor(fx)))
251 const y0 = Math.max(0, Math.min(H - 2, Math.floor(fy)))
252 const tx = Math.max(0, Math.min(1, fx - x0)), ty = Math.max(0, Math.min(1, fy - y0))
253 const i = y0 * W + x0
254 return (field[i] * (1 - tx) + field[i + 1] * tx) * (1 - ty)
255 + (field[i + W] * (1 - tx) + field[i + W + 1] * tx) * ty
256 }
257}
258
259// -------------------------------------------------------------------- colour
260/**
261 * Hypsometric tint for the Pacific Northwest. Sea is depth-shaded; land runs from
262 * the wet coastal green through the dry east-flank olive to rock and snow. The
263 * ramp is steeper below 800 m than a global one would be, because almost everyone
264 * in this world lives under 200 m and the Sound country needs the contrast.
265 */
266const LAND = [
267 [0, [72, 104, 66]], [120, [86, 118, 70]], [400, [104, 130, 76]],
268 [900, [126, 134, 84]], [1400, [140, 128, 92]], [1900, [138, 126, 114]],
269 [2600, [150, 148, 148]], [3400, [186, 190, 196]],
270]
271function tint(h) {
272 if (h <= LAND[0][0]) return LAND[0][1]
273 for (let i = 1; i < LAND.length; i++) {
274 if (h <= LAND[i][0]) {
275 const [a, ca] = LAND[i - 1], [b, cb] = LAND[i]
276 const t = (h - a) / (b - a)
277 return [0, 1, 2].map((k) => ca[k] + (cb[k] - ca[k]) * t)
278 }
279 }
280 return LAND[LAND.length - 1][1]
281}
282
283function colorise(height) {
284 const az = (315 * Math.PI) / 180, alt = (42 * Math.PI) / 180, strength = 0.55
285 const lx = Math.cos(alt) * Math.sin(az), ly = Math.cos(alt) * Math.cos(az)
286 const lz = Math.sin(alt)
287 const rgb = new Uint8Array(GRID_W * GRID_H * 3)
288 for (let y = 0; y < GRID_H; y++) {
289 for (let x = 0; x < GRID_W; x++) {
290 const i = y * GRID_W + x
291 const nx = (height[y * GRID_W + Math.max(0, x - 1)]
292 - height[y * GRID_W + Math.min(GRID_W - 1, x + 1)]) / (2 * MPP_X)
293 // theta wraps, so the shading does too — no seam in the hillshade.
294 const ny = (height[((y - 1 + GRID_H) % GRID_H) * GRID_W + x]
295 - height[((y + 1) % GRID_H) * GRID_W + x]) / (2 * MPP_Y)
296 const len = Math.hypot(nx, ny, 1)
297 let shade = (nx * lx + ny * ly + lz) / len
298 shade = 1 - strength + strength * Math.max(0, shade) * 1.35
299
300 const h = height[i]
301 let c
302 if (h <= 0) {
303 // Depth ramp tuned to the Sound and the Strait, not to the abyssal plain:
304 // most of this world's water is under 300 m and wants the whole gradient.
305 const d = Math.min(1, -h / 320)
306 c = [56 - 26 * d, 94 - 34 * d, 130 - 30 * d]
307 shade = 0.93 + 0.07 * shade
308 } else {
309 c = tint(h)
310 if (h > SNOWLINE_M) {
311 // Snow does not lie on cliffs; blending it out with slope is what turns a
312 // white blob into a mountain with rock showing through.
313 const slope = Math.hypot(nx, ny)
314 const holds = 1 - Math.min(1, Math.max(0, (slope - 0.5) / 0.5))
315 const t = Math.min(1, (h - SNOWLINE_M) / 700) * holds
316 c = c.map((v, k) => v + ([247, 249, 253][k] - v) * t)
317 }
318 }
319 rgb[i * 3] = Math.max(0, Math.min(255, c[0] * shade))
320 rgb[i * 3 + 1] = Math.max(0, Math.min(255, c[1] * shade))
321 rgb[i * 3 + 2] = Math.max(0, Math.min(255, c[2] * shade))
322 }
323 }
324 return rgb
325}
326
327// -------------------------------------------------------------------- process
328mkdirSync(CACHE, { recursive: true })
329console.log(`\nTower — Cascadia at 1:1, ${GRID_W}x${GRID_H} at ${MPP_X.toFixed(1)} m/px`)
330line(1, 'frame', `anchor ${ANCHOR.name}, rotation ${ROT_DEG}°, cut ${(CUT_SHIFT_M / 1000).toFixed(1)} km east`)
331
332const src = await loadElevation()
333const sample = sampler(src)
334
335const height = step(1, 'sample', () => {
336 const out = new Float32Array(GRID_W * GRID_H)
337 for (let y = 0; y < GRID_H; y++) {
338 for (let x = 0; x < GRID_W; x++) {
339 const [lat, lon] = fromUTM(...gridToUTM(x, y))
340 out[y * GRID_W + x] = sample(lat, lon)
341 }
342 }
343 return out
344})
345
346let lo = Infinity, hi = -Infinity, wet = 0
347for (const v of height) {
348 if (v < lo) lo = v
349 if (v > hi) hi = v
350 if (v <= 0) wet++
351}
352line(1, 'terrain', `${lo.toFixed(0)}..${hi.toFixed(0)} m, ` +
353 `${((1 - wet / height.length) * 100).toFixed(0)}% land`)
354
355/**
356 * The seam. theta wraps, so column 0 and column GRID_H-1 are neighbours on the
357 * cylinder even though they are 314 km apart on Earth. This mode's whole claim is
358 * that the ring closes on real ground, so the size of that step is the number that
359 * says whether the claim holds — print it rather than trusting the geometry.
360 */
361const seam = []
362for (let x = 0; x < GRID_W; x++) {
363 seam.push(Math.abs(height[x] - height[(GRID_H - 1) * GRID_W + x]))
364}
365const sorted = [...seam].sort((a, b) => a - b)
366line(1, 'seam step', `median ${sorted[GRID_W >> 1].toFixed(0)} m, ` +
367 `p90 ${sorted[Math.floor(GRID_W * 0.9)].toFixed(0)} m, max ${sorted[GRID_W - 1].toFixed(0)} m`)
368if (process.env.SEAM) {
369 // Where the cliff is, not just how big. A step that is small through the
370 // inhabited middle and large only in the far north is a different proposition
371 // from one that is uniformly bad, and the median alone cannot tell them apart.
372 for (let i = 0; i < 9; i++) {
373 const x = Math.min(GRID_W - 1, Math.round((i / 8) * GRID_W))
374 const [lat] = fromUTM(...gridToUTM(x, 0))
375 const band = seam.slice(Math.max(0, x - 64), x + 64)
376 const mean = band.reduce((a, b) => a + b, 0) / band.length
377 line(2, `${lat.toFixed(2)}N z=${(x * MPP_X / 1000).toFixed(0)} km`,
378 `${'#'.repeat(Math.round(mean / 60))} ${mean.toFixed(0)} m`)
379 }
380}
381
382const rgb = step(1, 'colourise', () => colorise(height))
383
384const span = HEIGHT_MAX_M - HEIGHT_MIN_M
385const norm = new Float64Array(height.length)
386/**
387 * Clipping the seafloor is deliberate and clipping the land is a bug, so they are
388 * counted apart. terrain.mjs fails on any clipping at all because every metre of
389 * its terrain was authored on purpose (D-024); here the abyssal plain is 2.6 km
390 * down and under water, and flattening it costs nothing anybody can see, while
391 * carrying it would coarsen every hillside in the world to pay for it.
392 */
393let drowned = 0, overflowed = 0
394for (let i = 0; i < norm.length; i++) {
395 const v = (height[i] - HEIGHT_MIN_M) / span
396 if (v < 0) drowned++
397 if (v > 1) overflowed++
398 norm[i] = Math.min(1, Math.max(0, v))
399}
400line(1, 'height range',
401 `${HEIGHT_MIN_M}..${HEIGHT_MAX_M} m, ${(span / 255).toFixed(1)} m per 8-bit step`)
402if (drowned) line(1, 'seafloor clamped', `${drowned} px below ${HEIGHT_MIN_M} m (intended)`)
403if (overflowed) {
404 console.error(`\n ! ${overflowed} pixel(s) of LAND above ${HEIGHT_MAX_M} m were ` +
405 `flattened; terrain reaches ${hi.toFixed(0)} m. Raise HEIGHT_MAX_M.\n`)
406 process.exitCode = 1
407}
408
409line(1, 'encode')
410const write = (name, make) =>
411 step(2, name, () => {
412 const buf = make()
413 writeFileSync(`static/${name}`, buf)
414 return buf
415 }, (b) => `${mb(b.length).padEnd(10)}`)
416write('cascadia-height.png', () => encodeGray16(GRID_W, GRID_H, norm))
417write('cascadia-color.png', () => encodeRGB(GRID_W, GRID_H, rgb))
418
419console.log(`\ndone in ${secs(T0)}\n`)