cascadia.mjs

1	/**
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	*/
23	import { readFileSync, writeFileSync, mkdirSync, existsSync } from 'node:fs'
24	import { 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	*/
35	const ROT_DEG = Number(process.env.ROT ?? 0)
36	/** A real place pinned to a grid pixel. ROT_DEG rotates the frame about this. */
37	const 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	*/
44	const 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
57	const GRID_W = 2048 // x is z, the axis. Sun cap at x = 0.
58	const GRID_H = 1192 // y is theta, and it wraps.
59	const LENGTH_M = 540_000
60	const CIRCUMFERENCE_M = 2 * Math.PI * 50_000
61	const MPP_X = LENGTH_M / GRID_W
62	const 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	*/
78	const HEIGHT_MIN_M = -560
79	const HEIGHT_MAX_M = 4400
80	const SNOWLINE_M = 1900
81
82	const ZOOM = 9 // ~206 m/px at 47.6°N, against the grid's 264 m/px
83	const TILE = 256
84	const CACHE = '.cache/terrain-tiles'
85	const TILE_URL = (z, x, y) =>
86	`https://s3.amazonaws.com/elevation-tiles-prod/terrarium/${z}/${x}/${y}.png`
87
88	// ------------------------------------------------------------------- progress
89	const T0 = Date.now()
90	const secs = (from) => `${((Date.now() - from) / 1000).toFixed(1)}s`
91	const mb = (n) => `${(n / 1048576).toFixed(1)} MB`
92	function 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	}
96	function 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.
105	const A = 6378137, F = 1 / 298.257223563, K0 = 0.9996, LON0 = (-123 * Math.PI) / 180
106	const E2 = F * (2 - F), EP2 = E2 / (1 - E2)
107	const M0 = 1 - E2 / 4 - (3 * E2 ** 2) / 64 - (5 * E2 ** 3) / 256
108
109	function 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
126	function 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
147	const BETA = (ROT_DEG * Math.PI) / 180
148	// +z (deepward) points BETA west of north; +theta (spinward) is 90° clockwise of it.
149	const ZE = -Math.sin(BETA), ZN = Math.cos(BETA)
150	const YE = Math.cos(BETA), YN = Math.sin(BETA)
151	const [AE, AN] = toUTM(ANCHOR.lat, ANCHOR.lon)
152	// Grid origin: back off from the anchor along both axes to reach pixel (0,0).
153	const OE = AE - (ANCHOR.x + 0.5) * MPP_X * ZE - ((ANCHOR.y + 0.5) * MPP_Y - CUT_SHIFT_M) * YE
154	const 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. */
157	const 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	*/
168	const lonToGX = (lon) => ((lon + 180) / 360) * TILE * 2 ** ZOOM
169	const 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
174	async 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. */
192	async 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. */
246	function 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	*/
266	const 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	]
271	function 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
283	function 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
328	mkdirSync(CACHE, { recursive: true })
329	console.log(`\nTower — Cascadia at 1:1, ${GRID_W}x${GRID_H} at ${MPP_X.toFixed(1)} m/px`)
330	line(1, 'frame', `anchor ${ANCHOR.name}, rotation ${ROT_DEG}°, cut ${(CUT_SHIFT_M / 1000).toFixed(1)} km east`)
331
332	const src = await loadElevation()
333	const sample = sampler(src)
334
335	const 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
346	let lo = Infinity, hi = -Infinity, wet = 0
347	for (const v of height) {
348	if (v < lo) lo = v
349	if (v > hi) hi = v
350	if (v <= 0) wet++
351	}
352	line(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	*/
361	const seam = []
362	for (let x = 0; x < GRID_W; x++) {
363	seam.push(Math.abs(height[x] - height[(GRID_H - 1) * GRID_W + x]))
364	}
365	const sorted = [...seam].sort((a, b) => a - b)
366	line(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`)
368	if (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
382	const rgb = step(1, 'colourise', () => colorise(height))
383
384	const span = HEIGHT_MAX_M - HEIGHT_MIN_M
385	const 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	*/
393	let drowned = 0, overflowed = 0
394	for (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	}
400	line(1, 'height range',
401	`${HEIGHT_MIN_M}..${HEIGHT_MAX_M} m, ${(span / 255).toFixed(1)} m per 8-bit step`)
402	if (drowned) line(1, 'seafloor clamped', `${drowned} px below ${HEIGHT_MIN_M} m (intended)`)
403	if (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
409	line(1, 'encode')
410	const 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)}`)
416	write('cascadia-height.png', () => encodeGray16(GRID_W, GRID_H, norm))
417	write('cascadia-color.png', () => encodeRGB(GRID_W, GRID_H, rgb))
418
419	console.log(`\ndone in ${secs(T0)}\n`)