1use super::{Event, SurveyRegion, fingerprint, provider_client, read_bounded, write_atomic};
2use anyhow::{Context as _, Result, bail, ensure};
3use serde::{Deserialize, Serialize};
4use std::{
5 collections::{BTreeMap, BTreeSet},
6 env, fs,
7 io::Cursor,
8 path::{Path, PathBuf},
9 sync::Arc,
10 time::Duration,
11};
12use trailgen_core::{
13 Coord, ElevationSample, ElevationSampler, Provenance, source::SourceFingerprint,
14};
15
16const DEFAULT_ENDPOINT: &str = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium";
17const PREFERRED_ZOOM: u8 = 12;
18const MINIMUM_ZOOM: u8 = 8;
19const MAX_TILES: usize = 256;
20const MAX_TILE_BYTES: u64 = 4 * 1024 * 1024;
21const PLAUSIBLE_ELEVATION_M: std::ops::RangeInclusive<f64> = -150.0..=9_000.0;
22
23#[derive(Clone, Copy, Debug, Deserialize, Eq, Ord, PartialEq, PartialOrd, Serialize)]
24pub struct TerrainTileId {
25 pub z: u8,
26 pub x: u32,
27 pub y: u32,
28}
29
30impl TerrainTileId {
31 fn relative_path(self) -> PathBuf {
32 PathBuf::from("sources/mapzen-terrain")
33 .join(self.z.to_string())
34 .join(self.x.to_string())
35 .join(format!("{}.png", self.y))
36 }
37
38 fn url(self, endpoint: &str) -> String {
39 format!(
40 "{}/{}/{}/{}.png",
41 endpoint.trim_end_matches('/'),
42 self.z,
43 self.x,
44 self.y
45 )
46 }
47}
48
49#[derive(Clone, Debug, Deserialize, Serialize)]
50pub struct TerrainReceipt {
51 pub tile: TerrainTileId,
52 pub raw_path: PathBuf,
53 pub raw: SourceFingerprint,
54}
55
56pub struct TerrainSource {
57 pub receipt: TerrainReceipt,
58 bytes: Vec<u8>,
59}
60
61#[derive(Clone, Debug)]
62pub struct TopographicTile {
63 pub id: TerrainTileId,
64 pub width: u32,
65 pub height: u32,
66 elevations_m: Arc<[f32]>,
67}
68
69impl TopographicTile {
70 #[must_use]
71 pub fn elevation(&self, x: u32, y: u32) -> Option<f32> {
72 (x < self.width && y < self.height)
73 .then(|| self.elevations_m[(y as usize * self.width as usize) + x as usize])
74 }
75}
76
77#[allow(
78 clippy::cast_possible_truncation,
79 reason = "sub-centimeter DEM precision is immaterial to ten-meter isohypses"
80)]
81pub fn topographic_tile(project: &Path, receipt: &TerrainReceipt) -> Result<TopographicTile> {
82 let path = project.join(&receipt.raw_path);
83 let bytes = fs::read(&path).with_context(|| format!("read {}", path.display()))?;
84 ensure!(
85 fingerprint(&bytes) == receipt.raw,
86 "topographic tile {} does not match its index receipt",
87 path.display()
88 );
89 let tile = Tile::decode(&bytes).with_context(|| format!("decode {}", path.display()))?;
90 let mut elevations_m = Vec::with_capacity(tile.width as usize * tile.height as usize);
91 for y in 0..tile.height {
92 for x in 0..tile.width {
93 elevations_m.push(tile.elevation(x, y) as f32);
94 }
95 }
96 Ok(TopographicTile {
97 id: receipt.tile,
98 width: tile.width,
99 height: tile.height,
100 elevations_m: elevations_m.into(),
101 })
102}
103
104pub fn acquire(
105 project: &Path,
106 regions: &[SurveyRegion],
107 fetch_missing: bool,
108 emit: &mut impl FnMut(Event),
109) -> Result<Vec<TerrainSource>> {
110 let tiles = desired_tiles(regions);
111 let total = tiles.len();
112 let endpoint =
113 env::var("TRAILGEN_TERRAIN_ENDPOINT").unwrap_or_else(|_| DEFAULT_ENDPOINT.to_owned());
114 let mut client = None;
115 let mut sources = Vec::with_capacity(tiles.len());
116 emit(Event::Elevating { complete: 0, total });
117 for (slot, tile) in tiles.into_iter().enumerate() {
118 let relative = tile.relative_path();
119 let path = project.join(&relative);
120 let cached = fs::read(&path)
121 .ok()
122 .filter(|bytes| Tile::decode(bytes).is_ok());
123 let bytes = if let Some(bytes) = cached {
124 bytes
125 } else {
126 ensure!(
127 fetch_missing,
128 "terrain tile z{}/x{}/y{} is absent or corrupt",
129 tile.z,
130 tile.x,
131 tile.y
132 );
133 let client = client.get_or_insert_with(|| {
134 provider_client("terrain", Duration::from_mins(1))
135 .expect("static terrain client configuration is valid")
136 });
137 let response = client
138 .get(tile.url(&endpoint))
139 .send()
140 .with_context(|| format!("fetch terrain tile z{}/x{}/y{}", tile.z, tile.x, tile.y))?
141 .error_for_status()
142 .with_context(|| {
143 format!(
144 "terrain provider rejected z{}/x{}/y{}",
145 tile.z, tile.x, tile.y
146 )
147 })?;
148 let bytes = read_bounded(response, MAX_TILE_BYTES, "terrain tile")?;
149 Tile::decode(&bytes).context("decode terrain provider response")?;
150 write_atomic(&path, &bytes)?;
151 bytes
152 };
153 sources.push(TerrainSource {
154 receipt: TerrainReceipt {
155 tile,
156 raw_path: relative,
157 raw: fingerprint(&bytes),
158 },
159 bytes,
160 });
161 emit(Event::Elevating {
162 complete: slot + 1,
163 total,
164 });
165 }
166 Ok(sources)
167}
168
169pub fn desired_tiles(regions: &[SurveyRegion]) -> Vec<TerrainTileId> {
170 for zoom in (MINIMUM_ZOOM..=PREFERRED_ZOOM).rev() {
171 let tiles = regions
172 .iter()
173 .flat_map(|region| tiles_for_bounds(region, zoom))
174 .collect::<BTreeSet<_>>();
175 if tiles.len() <= MAX_TILES || zoom == MINIMUM_ZOOM {
176 return tiles.into_iter().collect();
177 }
178 }
179 unreachable!("terrain zoom interval is nonempty")
180}
181
182fn tiles_for_bounds(region: &SurveyRegion, z: u8) -> Vec<TerrainTileId> {
183 let north_west = tile_at(Coord::new(region.bounds.west, region.bounds.north), z);
184 let south_east = tile_at(Coord::new(region.bounds.east, region.bounds.south), z);
185 (north_west.x..=south_east.x)
186 .flat_map(move |x| (north_west.y..=south_east.y).map(move |y| TerrainTileId { z, x, y }))
187 .collect()
188}
189
190#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
191fn tile_at(coord: Coord, z: u8) -> TerrainTileId {
192 let n = f64::from(1_u32 << z);
193 let x = ((coord.lon + 180.0) / 360.0 * n)
194 .floor()
195 .clamp(0.0, n - 1.0) as u32;
196 let latitude = coord.lat.clamp(-85.051_128_78, 85.051_128_78).to_radians();
197 let y = ((1.0 - latitude.tan().asinh() / std::f64::consts::PI) * 0.5 * n)
198 .floor()
199 .clamp(0.0, n - 1.0) as u32;
200 TerrainTileId { z, x, y }
201}
202
203pub struct TerrainAtlas {
204 zoom: u8,
205 tiles: BTreeMap<(u32, u32), Tile>,
206}
207
208impl TerrainAtlas {
209 pub fn decode(sources: &[TerrainSource]) -> Result<Option<Self>> {
210 let Some(zoom) = sources.first().map(|source| source.receipt.tile.z) else {
211 return Ok(None);
212 };
213 ensure!(
214 sources.iter().all(|source| source.receipt.tile.z == zoom),
215 "terrain corpus mixes zoom levels"
216 );
217 let tiles = sources
218 .iter()
219 .map(|source| {
220 let id = source.receipt.tile;
221 Ok(((id.x, id.y), Tile::decode(&source.bytes)?))
222 })
223 .collect::<Result<_>>()?;
224 Ok(Some(Self { zoom, tiles }))
225 }
226}
227
228impl ElevationSampler for TerrainAtlas {
229 #[allow(clippy::cast_possible_truncation)]
230 fn sample(&self, coord: Coord) -> Option<ElevationSample> {
231 let id = tile_at(coord, self.zoom);
232 let tile = self.tiles.get(&(id.x, id.y))?;
233 let n = f64::from(1_u32 << self.zoom);
234 let world_x = (coord.lon + 180.0) / 360.0 * n;
235 let latitude = coord.lat.clamp(-85.051_128_78, 85.051_128_78).to_radians();
236 let world_y = (1.0 - latitude.tan().asinh() / std::f64::consts::PI) * 0.5 * n;
237 let pixel_x = world_x.mul_add(f64::from(tile.width), -0.5);
238 let pixel_y = world_y.mul_add(f64::from(tile.height), -0.5);
239 let x_floor = pixel_x.floor();
240 let y_floor = pixel_y.floor();
241 let x0 = x_floor as i64;
242 let y0 = y_floor as i64;
243 let dx = pixel_x - x_floor;
244 let dy = pixel_y - y_floor;
245 let mut elevation = 0.0;
246 let mut weight = 0.0;
247 for (x, x_weight) in [(x0, 1.0 - dx), (x0 + 1, dx)] {
248 for (y, y_weight) in [(y0, 1.0 - dy), (y0 + 1, dy)] {
249 let pixel_weight = x_weight * y_weight;
250 if let Some(value) = self.pixel(x, y, tile.width, tile.height)
251 && PLAUSIBLE_ELEVATION_M.contains(&value)
252 {
253 elevation = value.mul_add(pixel_weight, elevation);
254 weight += pixel_weight;
255 }
256 }
257 }
258 (weight > f64::EPSILON).then(|| ElevationSample {
259 ele_m: elevation / weight,
260 confidence: 0.82,
261 provenance: Provenance {
262 source: "mapzen-terrain-tiles".to_owned(),
263 layer: Some(format!("terrarium-z{}", self.zoom)),
264 source_id: Some(format!("{}/{}/{}", id.z, id.x, id.y)),
265 license: Some(
266 "source-specific; https://github.com/tilezen/joerd/blob/master/docs/attribution.md"
267 .to_owned(),
268 ),
269 },
270 })
271 }
272}
273
274impl TerrainAtlas {
275 fn pixel(&self, x: i64, y: i64, width: u32, height: u32) -> Option<f64> {
276 let side = i64::from(1_u32 << self.zoom);
277 let world_width = side * i64::from(width);
278 let world_height = side * i64::from(height);
279 let x = x.rem_euclid(world_width);
280 let y = y.clamp(0, world_height - 1);
281 let tile_x = u32::try_from(x / i64::from(width)).ok()?;
282 let tile_y = u32::try_from(y / i64::from(height)).ok()?;
283 let tile = self.tiles.get(&(tile_x, tile_y))?;
284 if tile.width != width || tile.height != height {
285 return None;
286 }
287 Some(tile.elevation(
288 u32::try_from(x % i64::from(width)).ok()?,
289 u32::try_from(y % i64::from(height)).ok()?,
290 ))
291 }
292}
293
294struct Tile {
295 width: u32,
296 height: u32,
297 rgb: Vec<u8>,
298}
299
300impl Tile {
301 fn decode(bytes: &[u8]) -> Result<Self> {
302 let mut decoder = png::Decoder::new(Cursor::new(bytes));
303 decoder.set_transformations(png::Transformations::EXPAND | png::Transformations::STRIP_16);
304 let mut reader = decoder.read_info().context("read terrain PNG header")?;
305 let mut buffer = vec![
306 0;
307 reader
308 .output_buffer_size()
309 .context("terrain PNG is too large")?
310 ];
311 let info = reader
312 .next_frame(&mut buffer)
313 .context("decode terrain PNG")?;
314 let pixels = &buffer[..info.buffer_size()];
315 let rgb = match info.color_type {
316 png::ColorType::Rgb => pixels.to_vec(),
317 png::ColorType::Rgba => pixels
318 .chunks_exact(4)
319 .flat_map(|pixel| pixel[..3].iter().copied())
320 .collect(),
321 other => bail!("terrain PNG has unsupported color type {other:?}"),
322 };
323 ensure!(
324 rgb.len() == info.width as usize * info.height as usize * 3,
325 "terrain PNG pixel count is inconsistent"
326 );
327 Ok(Self {
328 width: info.width,
329 height: info.height,
330 rgb,
331 })
332 }
333
334 fn elevation(&self, x: u32, y: u32) -> f64 {
335 let offset = (y as usize * self.width as usize + x as usize) * 3;
336 f64::from(self.rgb[offset]).mul_add(
337 256.0,
338 f64::from(self.rgb[offset + 1]) + f64::from(self.rgb[offset + 2]) / 256.0,
339 ) - 32_768.0
340 }
341}
342
343#[cfg(test)]
344mod tests {
345 use super::*;
346 use trailgen_core::source::GeoBounds;
347
348 #[test]
349 fn bounded_regions_choose_one_common_zoom_and_tile_set() {
350 let regions = [
351 SurveyRegion::new(GeoBounds::new(-74.15, 41.15, -74.0, 41.30))
352 .expect("valid Harriman bounds"),
353 ];
354 let tiles = desired_tiles(®ions);
355 assert!(!tiles.is_empty() && tiles.len() <= MAX_TILES);
356 assert!(tiles.iter().all(|tile| tile.z == tiles[0].z));
357 }
358
359 #[test]
360 fn void_and_bathymetric_pixels_do_not_poison_hiking_profiles() {
361 let atlas = TerrainAtlas {
362 zoom: 0,
363 tiles: BTreeMap::from([(
364 (0, 0),
365 Tile {
366 width: 1,
367 height: 1,
368 rgb: vec![107, 0, 0],
369 },
370 )]),
371 };
372 assert!(atlas.sample(Coord::new(-74.0, 41.0)).is_none());
373 }
374}