Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_geom/examples/world_gen.rs

18.6 KiB, 1 run

created by r1870400018:59716, which is this file's identity for as long as the history lasts, whatever it is later renamed to

download · who wrote it · its history

1//! Draws a world file from Natural Earth, and optionally a picture of it.
2//!
3//! ```text
4//! cargo run --release -p oxedyne_fe2o3_geom --example world_gen -- \
5//! [--ne DIR] [--spec ochre] [--level TOL,ISLAND,LAKE,BORDER]... [--detail N] \
6//! [--places FILE --rank-max R] [--out FILE] [--compare FILE] [--svg FILE]
7//! ```
8//!
9//! Natural Earth is public domain. The GeoJSON is read from `DIR` (default
10//! `~/.cache/natural-earth`), and a missing layer is named with the address to fetch it from;
11//! nothing here reaches the network.
12//!
13//! `--spec ochre` (the default) is the three levels Ochre ships: thirty kilometres, six and
14//! one, each with its least island in square kilometres, least lake and least border in
15//! kilometres. `--level` replaces them, once per level, coarsest first. `--detail N` keeps
16//! only level `N`. `--compare FILE` checks the result against another world file byte for
17//! byte after the eight-byte magic, which is how this is held to Ochre's `world.bin`.
18//!
19//! The numbers are Ochre's `dev/gen_world.py`'s, and so is the arithmetic, down to the order
20//! of the additions: a ring's area decides whether it is drawn, and a rounding difference at a
21//! threshold would drop an island. Python's `sum` of floats is compensated since 3.12 and its
22//! `%` takes the sign of the divisor, and both are reproduced below.
23
24use oxedyne_fe2o3_geom::{
25 cell::Cell,
26 proj::{
27 EARTH_RADIUS_M,
28 Projection,
29 RingMode,
30 ScreenPaths,
31 Viewport,
32 unit_vec,
33 },
34 world::{
35 self,
36 Label,
37 Layer,
38 LayerKind,
39 World,
40 },
41};
42
43use oxedyne_fe2o3_core::prelude::*;
44use oxedyne_fe2o3_jdat::{
45 prelude::*,
46 string::dec::DecoderConfig,
47 usr::{
48 UsrKind,
49 UsrKindCode,
50 UsrKindId,
51 },
52};
53
54use std::{
55 collections::BTreeMap,
56 fmt::Write as _,
57};
58
59const NE_URL: &str = "https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/";
60
61// Ochre's levels, as (tolerance m, least island km2, least lake km2, least border km).
62const OCHRE_LEVELS: [[f64; 4]; 3] = [
63 [30_000.0, 2_000.0, 20_000.0, 200.0],
64 [ 6_000.0, 100.0, 2_000.0, 40.0],
65 [ 1_000.0, 2.0, 550.0, 5.0],
66];
67
68// The layers, in file order: name, kind, sources, which of a level's thresholds applies.
69const LAYERS: [(&str, LayerKind, &[&str], usize); 3] = [
70 ("land", LayerKind::Fill, &["ne_10m_land", "ne_10m_minor_islands"], 1),
71 ("lakes", LayerKind::Fill, &["ne_10m_lakes"], 2),
72 ("borders", LayerKind::Stroke, &["ne_10m_admin_0_boundary_lines_land"], 3),
73];
74
75/// A GeoJSON coordinate as Python would hold it: the value, and whether it was written as a
76/// whole number, which Python's `sum` treats differently.
77#[derive(Clone, Copy)]
78struct Num {
79 v: f64,
80 int: bool,
81}
82
83type Ring = Vec<(Num, Num)>; // (lng, lat), as GeoJSON orders them
84
85struct Args {
86 ne: String,
87 levels: Vec<[f64; 4]>,
88 detail: Option<u8>,
89 places: Option<String>,
90 rank_max: u8,
91 out: Option<String>,
92 compare: Option<String>,
93 svg: Option<String>,
94}
95
96fn args() -> Outcome<Args> {
97 let home = std::env::var("HOME").unwrap_or_default();
98 let mut a = Args {
99 ne: fmt!("{}/.cache/natural-earth", home), levels: Vec::new(), detail: None,
100 places: None, rank_max: 3, out: None, compare: None, svg: None,
101 };
102 let mut it = std::env::args().skip(1);
103 while let Some(flag) = it.next() {
104 let mut val = || -> Outcome<String> {
105 it.next().ok_or_else(|| err!("{} wants a value.", flag; Input, Missing))
106 };
107 match flag.as_str() {
108 "--ne" => a.ne = res!(val()),
109 "--spec" => {
110 let s = res!(val());
111 if s != "ochre" {
112 return Err(err!("No spec called {:?}; there is `ochre`.", s; Input, Invalid));
113 }
114 },
115 "--level" => {
116 let s = res!(val());
117 let parts: Vec<f64> = res!(s.split(',').map(|p| p.trim().parse::<f64>())
118 .collect::<Result<Vec<f64>, _>>(), Input, Invalid);
119 if parts.len() != 4 {
120 return Err(err!("--level {:?} wants four numbers.", s; Input, Invalid));
121 }
122 a.levels.push([parts[0], parts[1], parts[2], parts[3]]);
123 },
124 "--detail" => a.detail = Some(res!(res!(val()).parse::<u8>(), Input, Invalid)),
125 "--places" => a.places = Some(res!(val())),
126 "--rank-max" => a.rank_max = res!(res!(val()).parse::<u8>(), Input, Invalid),
127 "--out" => a.out = Some(res!(val())),
128 "--compare" => a.compare = Some(res!(val())),
129 "--svg" => a.svg = Some(res!(val())),
130 other => return Err(err!("Unknown argument {:?}.", other; Input, Invalid)),
131 }
132 }
133 if a.levels.is_empty() {
134 a.levels = OCHRE_LEVELS.to_vec();
135 }
136 Ok(a)
137}
138
139fn main() -> Outcome<()> {
140 let a = res!(args());
141 let mut sources: BTreeMap<String, Vec<Ring>> = BTreeMap::new();
142 for (_, _, names, _) in LAYERS.iter() {
143 for name in names.iter() {
144 if !sources.contains_key(*name) {
145 let path = fmt!("{}/{}.geojson", a.ne, name);
146 let dat = res!(load(&path, name));
147 sources.insert(name.to_string(), res!(rings_of_file(&dat, name)));
148 }
149 }
150 }
151
152 let mut w = World::default();
153 for (name, kind, names, which) in LAYERS.iter() {
154 for (detail, level) in a.levels.iter().enumerate() {
155 if a.detail.map_or(false, |d| d as usize != detail) {
156 continue;
157 }
158 let mut rings: Vec<Vec<(f64, f64)>> = Vec::new();
159 for src in names.iter() {
160 if let Some(list) = sources.get(*src) {
161 prepare(list, *kind, level[0], level[*which], &mut rings);
162 }
163 }
164 w.layers.push(Layer {
165 name: name.to_string(), kind: *kind, detail: detail as u8, tol_m: level[0],
166 rings, labels: Vec::new(),
167 });
168 }
169 }
170 if let Some(path) = &a.places {
171 let dat = res!(load(path, "ne_10m_populated_places"));
172 let labels = res!(places(&dat, a.rank_max));
173 w.layers.push(Layer {
174 name: "places".to_string(), kind: LayerKind::Label, detail: 0, tol_m: 0.0,
175 rings: Vec::new(), labels,
176 });
177 }
178
179 let bytes = res!(world::write(&w));
180 report(&w, bytes.len());
181 // What was written must read back as what was meant.
182 let back = res!(world::read(&bytes));
183 let again = res!(world::write(&back));
184 if again != bytes {
185 return Err(err!("The world did not survive a round trip through its own file."; Bug));
186 }
187 if let Some(path) = &a.out {
188 res!(std::fs::write(path, &bytes), File, Write);
189 println!("wrote {} ({} bytes)", path, bytes.len());
190 }
191 if let Some(path) = &a.compare {
192 let other = res!(std::fs::read(path), File, Read);
193 let same = other.len() == bytes.len() && other.get(8..) == bytes.get(8..);
194 if same {
195 println!("byte-identical to {} after the magic ({} bytes)", path, bytes.len());
196 } else {
197 let at = other.iter().skip(8).zip(bytes.iter().skip(8)).position(|(x, y)| x != y);
198 return Err(err!("Differs from {}: {} bytes against {}, first difference at {:?}.",
199 path, bytes.len(), other.len(), at.map(|p| p + 8); Mismatch));
200 }
201 }
202 if let Some(path) = &a.svg {
203 res!(svg(&back, path));
204 println!("wrote {}", path);
205 }
206 Ok(())
207}
208
209fn load(path: &str, name: &str) -> Outcome<Dat> {
210 let text = match std::fs::read_to_string(path) {
211 Ok(t) => t,
212 Err(e) => return Err(err!(e,
213 "Cannot read {}. Fetch it once with\n curl -o {} {}{}.geojson",
214 path, path, NE_URL, name; File, Read)),
215 };
216 let cfg = DecoderConfig::<BTreeMap<UsrKindCode, UsrKind>, BTreeMap<String, UsrKindId>>::json(None);
217 Dat::decode_string_with_config(text, &cfg)
218}
219
220fn num(d: &Dat) -> Outcome<Num> {
221 let int = match d {
222 Dat::U8(_) | Dat::U16(_) | Dat::U32(_) | Dat::U64(_) |
223 Dat::I8(_) | Dat::I16(_) | Dat::I32(_) | Dat::I64(_) => true,
224 _ => false,
225 };
226 match d.get_float64() {
227 Some(f) => Ok(Num { v: f.0, int }),
228 None => Err(err!("A coordinate {:?} is not a number.", d; Input, Mismatch)),
229 }
230}
231
232fn position(d: &Dat) -> Outcome<(Num, Num)> {
233 match d {
234 Dat::List(v) if v.len() >= 2 => Ok((res!(num(&v[0])), res!(num(&v[1])))),
235 other => Err(err!("A position {:?} is not a pair.", other.kind(); Input, Mismatch)),
236 }
237}
238
239fn line(d: &Dat) -> Outcome<Ring> {
240 match d {
241 Dat::List(v) => {
242 let mut out = Vec::with_capacity(v.len());
243 for p in v {
244 out.push(res!(position(p)));
245 }
246 Ok(out)
247 },
248 other => Err(err!("A line {:?} is not a list.", other.kind(); Input, Mismatch)),
249 }
250}
251
252fn list(d: &Dat) -> Outcome<&Vec<Dat>> {
253 match d {
254 Dat::List(v) => Ok(v),
255 other => Err(err!("Expected a list, found {:?}.", other.kind(); Input, Mismatch)),
256 }
257}
258
259/// Every ring or line of every feature, whichever geometry it carries, in file order.
260fn rings_of_file(dat: &Dat, name: &str) -> Outcome<Vec<Ring>> {
261 let features = res!(dat.map_get_list(&dat!("features")));
262 let mut out = Vec::new();
263 for f in features {
264 let geom = match res!(f.map_get(&dat!("geometry"))) {
265 Some(g @ Dat::Map(_)) | Some(g @ Dat::OrdMap(_)) => g,
266 _ => continue,
267 };
268 let kind = res!(geom.map_get_string(&dat!("type")));
269 let coords = res!(geom.map_get_must(&dat!("coordinates")));
270 match kind.as_str() {
271 "MultiPolygon" => for poly in res!(list(coords)) {
272 for r in res!(list(poly)) {
273 out.push(res!(line(r)));
274 }
275 },
276 "Polygon" | "MultiLineString" => for r in res!(list(coords)) {
277 out.push(res!(line(r)));
278 },
279 "LineString" => out.push(res!(line(coords))),
280 _ => (),
281 }
282 }
283 println!("{}: {} rings", name, out.len());
284 Ok(out)
285}
286
287/// Python's `sum` over floats and small integers: Neumaier's compensated sum for the floats,
288/// integers added plainly, and an integer prefix summed exactly before the first float.
289fn py_sum(vals: impl Iterator<Item = Num>) -> f64 {
290 let mut i_sum: i64 = 0;
291 let mut floating = false;
292 let mut f = 0.0f64;
293 let mut c = 0.0f64;
294 for n in vals {
295 if !floating {
296 if n.int {
297 i_sum += n.v as i64;
298 continue;
299 }
300 floating = true;
301 f = i_sum as f64 + n.v;
302 continue;
303 }
304 if n.int {
305 f += n.v;
306 continue;
307 }
308 let t = f + n.v;
309 if f.abs() >= n.v.abs() {
310 c += (f - t) + n.v;
311 } else {
312 c += (n.v - t) + f;
313 }
314 f = t;
315 }
316 if !floating {
317 return i_sum as f64;
318 }
319 if c != 0.0 && c.is_finite() {
320 f += c;
321 }
322 f
323}
324
325/// Python's `%` on floats, which takes the sign of the divisor.
326fn py_mod(a: f64, b: f64) -> f64 {
327 let m = a % b;
328 if m != 0.0 {
329 if (b < 0.0) != (m < 0.0) { m + b } else { m }
330 } else {
331 0.0f64.copysign(b)
332 }
333}
334
335/// Roughly how much ground a ring encloses, in square kilometres: the shoelace formula in a
336/// plane laid on the ring's own middle latitude.
337fn area_km2(ring: &Ring) -> f64 {
338 let n = ring.len();
339 if n < 3 {
340 return 0.0;
341 }
342 let mid = py_sum(ring.iter().map(|c| c.1)) / n as f64;
343 let k = mid.to_radians().cos();
344 let mut a = 0.0;
345 for i in 0..n {
346 let (x1, y1) = (ring[i].0.v * k, ring[i].1.v);
347 let (x2, y2) = (ring[(i + 1) % n].0.v * k, ring[(i + 1) % n].1.v);
348 a += x1 * y2 - x2 * y1;
349 }
350 a.abs() / 2.0 * 111.32f64.powf(2.0)
351}
352
353/// Roughly how long a line is, in kilometres.
354fn length_km(ring: &Ring) -> f64 {
355 let mut total = 0.0;
356 for i in 0..ring.len().saturating_sub(1) {
357 let k = ((ring[i].1.v + ring[i + 1].1.v) / 2.0).to_radians().cos();
358 let dx = (py_mod(ring[i + 1].0.v - ring[i].0.v + 180.0, 360.0) - 180.0) * k;
359 let dy = ring[i + 1].1.v - ring[i].1.v;
360 total += dx.hypot(dy);
361 }
362 total * 111.32
363}
364
365/// Every ring of a source, pruned by size and simplified, appended as `(lat, lng)`.
366fn prepare(src: &[Ring], kind: LayerKind, eps_m: f64, least: f64, out: &mut Vec<Vec<(f64, f64)>>) {
367 for ring in src {
368 if ring.len() < 2 {
369 continue;
370 }
371 let small = match kind {
372 LayerKind::Fill => area_km2(ring) < least,
373 _ => length_km(ring) < least,
374 };
375 if small {
376 continue;
377 }
378 let vecs: Vec<[f64; 3]> = ring.iter().map(|(lng, lat)| unit_vec(lat.v, lng.v)).collect();
379 let keep = world::simplify_sphere(&vecs, eps_m);
380 let least_n = if kind == LayerKind::Fill { 4 } else { 2 };
381 if keep.len() < least_n {
382 continue;
383 }
384 out.push(keep.iter().map(|i| (ring[*i].1.v, ring[*i].0.v)).collect());
385 }
386}
387
388/// Populated places up to a rank, most important first.
389fn places(dat: &Dat, rank_max: u8) -> Outcome<Vec<Label>> {
390 let features = res!(dat.map_get_list(&dat!("features")));
391 let mut out = Vec::new();
392 for f in features {
393 let props = res!(f.map_get_must(&dat!("properties")));
394 let rank = match res!(props.map_get(&dat!("SCALERANK"))) {
395 Some(r) => r,
396 None => res!(props.map_get_must(&dat!("scalerank"))),
397 };
398 let rank = match rank.get_float64() {
399 Some(r) => r.0,
400 None => continue,
401 };
402 if !(rank >= 0.0 && rank <= rank_max as f64) {
403 continue;
404 }
405 let name = match res!(props.map_get(&dat!("NAME"))) {
406 Some(Dat::Str(s)) => s.clone(),
407 _ => res!(props.map_get_string(&dat!("name"))),
408 };
409 let geom = res!(f.map_get_must(&dat!("geometry")));
410 let (lng, lat) = res!(position(res!(geom.map_get_must(&dat!("coordinates")))));
411 out.push(Label { lat: lat.v, lng: lng.v, rank: rank as u8, name });
412 }
413 out.sort_by(|a, b| a.rank.cmp(&b.rank).then(a.name.cmp(&b.name)));
414 Ok(out)
415}
416
417fn report(w: &World, size: usize) {
418 println!("{:<9} {:<6} {:>9} {:>7} {:>9}", "layer", "drawn", "tolerance", "rings", "vertices");
419 for l in &w.layers {
420 let verts: usize = l.rings.iter().map(|r| r.len()).sum();
421 println!("{:<9} {:<6} {:>7.0} m {:>7} {:>9}",
422 l.name, fmt!("{:?}", l.kind).to_lowercase(), l.tol_m,
423 if l.kind == LayerKind::Label { l.labels.len() } else { l.rings.len() }, verts);
424 }
425 println!("total {:.1} KiB", size as f64 / 1024.0);
426}
427
428// ---------------------------------------------------------------------------------------------
429// The picture
430// ---------------------------------------------------------------------------------------------
431
432fn path_d(paths: &ScreenPaths, dx: f64, dy: f64) -> String {
433 let mut d = String::new();
434 for i in 0..paths.len() {
435 if let Some((pts, closed)) = paths.path(i) {
436 for (k, p) in pts.chunks(2).enumerate() {
437 let _ = write!(d, "{}{:.1} {:.1}", if k == 0 { "M" } else { "L" },
438 p[0] as f64 + dx, p[1] as f64 + dy);
439 }
440 if closed {
441 d.push('Z');
442 }
443 }
444 }
445 d
446}
447
448/// The world on a globe and on a flat map, with level-3 cells over both.
449fn svg(w: &World, path: &str) -> Outcome<()> {
450 let detail = if w.details() > 1 { 1 } else { 0 };
451 let (gw, mw, h) = (620.0, 1000.0, 620.0);
452 let globe = res!(Viewport::new(Projection::Orthographic, 10.0, 110.0, 0.0,
453 EARTH_RADIUS_M / 290.0, gw, h));
454 let flat = res!(Viewport::new(Projection::WebMercator, 20.0, 150.0, 0.0,
455 std::f64::consts::TAU * EARTH_RADIUS_M / mw * (20.0f64).to_radians().cos(), mw, h));
456 let mut cells: Vec<Vec<[f64; 3]>> = Vec::new();
457 for face in 0..6u8 {
458 for i in 0..8u32 {
459 for j in 0..8u32 {
460 cells.push(res!(Cell::from_face_ij(face, 3, i, j)).outline(8));
461 }
462 }
463 }
464 let mut s = String::new();
465 let _ = write!(s, "<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"{}\" height=\"{}\" \
466 viewBox=\"0 0 {} {}\">\n<rect width=\"100%\" height=\"100%\" fill=\"#f4f1ea\"/>\n",
467 gw + mw + 20.0, h + 40.0, gw + mw + 20.0, h + 40.0);
468 let _ = write!(s, "<text x=\"10\" y=\"24\" font-family=\"sans-serif\" font-size=\"15\">\
469 Natural Earth, level of detail {} ({:.0} km), with level-3 cells: orthographic globe and \
470 Web Mercator map</text>\n", detail,
471 w.layer("land", detail).map_or(0.0, |l| l.tol_m / 1000.0));
472 for (view, dx) in [(globe, 0.0), (flat, gw + 20.0)] {
473 let dy = 40.0;
474 // Each panel is clipped to itself, since the flat map carries a margin off its edge.
475 let _ = write!(s, "<clipPath id=\"p{}\"><rect x=\"{}\" y=\"{}\" width=\"{}\" \
476 height=\"{}\"/></clipPath>\n<g clip-path=\"url(#p{})\">\n",
477 dx as u32, dx, dy, view.w, view.h, dx as u32);
478 if view.kind == Projection::Orthographic {
479 let _ = write!(s, "<circle cx=\"{}\" cy=\"{}\" r=\"290\" fill=\"#bcd9ea\"/>\n",
480 dx + gw / 2.0, dy + h / 2.0);
481 } else {
482 let _ = write!(s, "<rect x=\"{}\" y=\"{}\" width=\"{}\" height=\"{}\" fill=\"#bcd9ea\"/>\n",
483 dx, dy, mw, h);
484 }
485 for (name, fill, stroke, mode) in [
486 ("land", "#e8e2cf", "#8a8270", RingMode::Fill),
487 ("lakes", "#bcd9ea", "none", RingMode::Fill),
488 ("borders", "none", "#a39a86", RingMode::Line),
489 ] {
490 if let Some(layer) = w.layer(name, detail) {
491 let mut out = ScreenPaths::new();
492 res!(view.project_rings(&layer.unit_rings(), mode, layer.tol_m, 0.3, &mut out));
493 let _ = write!(s, "<path d=\"{}\" fill=\"{}\" stroke=\"{}\" stroke-width=\"0.6\" \
494 fill-rule=\"nonzero\"/>\n", path_d(&out, dx, dy), fill, stroke);
495 }
496 }
497 let mut out = ScreenPaths::new();
498 res!(view.project_rings(&cells, RingMode::Outline, 0.0, 0.3, &mut out));
499 let _ = write!(s, "<path d=\"{}\" fill=\"none\" stroke=\"#0096b4\" stroke-width=\"0.7\" \
500 stroke-opacity=\"0.8\"/>\n", path_d(&out, dx, dy));
501 s.push_str("</g>\n");
502 }
503 s.push_str("</svg>\n");
504 res!(std::fs::write(path, s), File, Write);
505 Ok(())
506}