reproject_terrain/
reproject_terrain.rs

1//! Real-data web-mercator → geodetic-TMS reprojection smoke test.
2//!
3//! Targets a geodetic TMS tile, reprojects the covering web-mercator
4//! Terrarium DEM tiles (read from a directory of `z_x_y.png` files) onto its
5//! geodetic grid, encodes a `.terrain`, and — if given a real production
6//! tile for the same TMS coordinate — prints both for comparison.
7//!
8//! ```text
9//! cargo run -p terrain-codec --example reproject_terrain --features png -- \
10//!     <tms_z> <tms_x> <tms_y> <src_zoom> <tile_size> <png_dir> [real.terrain]
11//! ```
12
13use std::env;
14use std::path::PathBuf;
15
16use terrain_codec::heightmap::container::decode_image;
17use terrain_codec::heightmap::terrarium;
18use terrain_codec::mercator::MercatorDem;
19use terrain_codec::quantized_mesh::{DecodedMesh, TileBounds};
20use terrain_codec::terrain::{NormalMode, TerrainOptions, encode_terrain};
21use terrain_codec::tile_coords::geodetic_tms;
22
23fn main() {
24    let mut a = env::args().skip(1);
25    let tz: u8 = a.next().expect("tms_z").parse().unwrap();
26    let tx: u32 = a.next().expect("tms_x").parse().unwrap();
27    let ty: u32 = a.next().expect("tms_y").parse().unwrap();
28    let src_zoom: u8 = a.next().expect("src_zoom").parse().unwrap();
29    let tile_size: u32 = a.next().expect("tile_size").parse().unwrap();
30    let png_dir = PathBuf::from(a.next().expect("png_dir"));
31    let real = a.next();
32
33    // Target geodetic TMS tile.
34    let (w, s, e, n) = geodetic_tms::tile_to_bounds(tz, tx, ty);
35    let bounds = TileBounds::new(w, s, e, n);
36    let grid_size = tile_size + 1; // 2^n + 1
37    println!("target geodetic TMS z{tz}/{tx}/{ty}  bounds = [{w:.5}, {s:.5}, {e:.5}, {n:.5}]");
38
39    // Widen by ~one geodetic cell so the halo (buffer=1) stays in coverage.
40    let cell_lon = (e - w) / (grid_size - 1) as f64;
41    let cell_lat = (n - s) / (grid_size - 1) as f64;
42    let (x0, y0, ntx, nty) = MercatorDem::tiles_covering(
43        src_zoom,
44        w - cell_lon,
45        s - cell_lat,
46        e + cell_lon,
47        n + cell_lat,
48    );
49    println!(
50        "covering web-mercator z{src_zoom} tiles: x {x0}..{}  y {y0}..{}  ({} tiles, {tile_size}px)",
51        x0 + ntx - 1,
52        y0 + nty - 1,
53        ntx * nty
54    );
55
56    // Build the mercator DEM by reading each covering tile from disk.
57    let dem = MercatorDem::from_tiles(src_zoom, x0, y0, ntx, nty, tile_size, |z, x, y| {
58        let path = png_dir.join(format!("{z}_{x}_{y}.png"));
59        let bytes = std::fs::read(&path).unwrap_or_else(|_| panic!("missing tile {path:?}"));
60        let img = decode_image(&bytes).expect("decode png");
61        assert_eq!(img.width, tile_size);
62        assert_eq!(img.height, tile_size);
63        terrarium::decode(&img.rgb, img.width, img.height)
64    });
65
66    // Reproject onto the geodetic grid + a halo grid for seamless normals.
67    let grid = dem.geodetic_grid(&bounds, grid_size);
68    let buffered = dem.buffered_geodetic(&bounds, grid_size, 1);
69
70    let (gmin, gmax) = min_max(&grid);
71    println!("reprojected geodetic grid {grid_size}×{grid_size}: {gmin:.1} .. {gmax:.1} m");
72
73    let bytes = encode_terrain(
74        &grid,
75        grid_size,
76        &bounds,
77        &TerrainOptions {
78            max_error: 4.0,
79            compression_level: 6,
80            normals: NormalMode::BufferedGradient(buffered),
81            ..Default::default()
82        },
83    );
84    let mesh = DecodedMesh::decode(&bytes).expect("decode our output");
85    println!(
86        "ours: {} bytes (gzip), {} verts, {} tris, height {:.1}..{:.1} m, normals={}",
87        bytes.len(),
88        mesh.vertices.len(),
89        mesh.indices.len() / 3,
90        mesh.header.min_height,
91        mesh.header.max_height,
92        mesh.extensions.normals.is_some(),
93    );
94
95    if let Some(path) = real {
96        let raw = std::fs::read(&path).expect("read real .terrain");
97        let rm = DecodedMesh::decode(&raw).expect("decode real");
98        println!(
99            "real: {} bytes, {} verts, {} tris, height {:.1}..{:.1} m, normals={}  ({})",
100            raw.len(),
101            rm.vertices.len(),
102            rm.indices.len() / 3,
103            rm.header.min_height,
104            rm.header.max_height,
105            rm.extensions.normals.is_some(),
106            path,
107        );
108        println!(
109            "note: height offset vs real is expected — real is EGM2008 geoid-blended, ours is raw Terrarium (ellipsoidal)."
110        );
111    }
112
113    println!("\nOK ✅  web-mercator → geodetic-TMS reprojection produced valid .terrain");
114}
115
116fn min_max(v: &[f32]) -> (f32, f32) {
117    let mut lo = f32::INFINITY;
118    let mut hi = f32::NEG_INFINITY;
119    for &x in v {
120        lo = lo.min(x);
121        hi = hi.max(x);
122    }
123    (lo, hi)
124}
reproject_terrain/reproject_terrain.rs

reproject_terrain/
reproject_terrain.rs
