Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_geom/src/tile.rs

15.7 KiB, 3 runs

created by r1870400018:60371, 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//! Web Mercator map tiles: the `z/x/y` grid street maps are cut into.
2//!
3//! Zoom `z` divides the square Web Mercator world into `2^z` by `2^z` tiles, `x` eastward from
4//! the antimeridian and `y` southward from the map's top edge at
5//! [`crate::proj::WEB_MERCATOR_MAX_LAT`]. This is the grid of every slippy map, of Mapbox
6//! Vector Tiles and of PMTiles archives ([`pmtiles`]).
7//!
8//! A [`Viewport`] says which tiles it needs ([`Viewport::tiles_covering`]) and how to paint a
9//! tile's bitmap onto the screen ([`Viewport::tile_affine`]): exactly on the flat map, and on
10//! the globe as the flat map's transform with the error that costs, so a caller can paint
11//! bitmaps where the error is below a pixel and reproject vectors where it is not.
12
13pub mod pmtiles;
14
15use crate::proj::{
16 Projection,
17 Viewport,
18 unit_vec,
19};
20
21use oxedyne_fe2o3_core::prelude::*;
22
23use std::f64::consts::{
24 FRAC_PI_2,
25 PI,
26 TAU,
27};
28
29pub const MAX_ZOOM: u8 = 31; // the deepest zoom whose PMTiles ids fit a u64
30
31/// One tile of the grid.
32#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
33pub struct TileId {
34 pub z: u8,
35 pub x: u32, // 0 at the antimeridian, eastward
36 pub y: u32, // 0 at the top of the map, southward
37}
38
39impl std::fmt::Display for TileId {
40 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
41 write!(f, "{}/{}/{}", self.z, self.x, self.y)
42 }
43}
44
45/// A longitude's position across the grid at zoom `z`, in tiles: 0 at 180 degrees west.
46pub fn lng_to_tx(lng: f64, z: u8) -> f64 {
47 (lng + 180.0) / 360.0 * (1u64 << z) as f64
48}
49
50/// A latitude's position down the grid at zoom `z`, in tiles: 0 at the top of the map.
51pub fn lat_to_ty(lat: f64, z: u8) -> f64 {
52 let phi = lat.to_radians();
53 (1.0 - phi.tan().asinh() / PI) / 2.0 * (1u64 << z) as f64
54}
55
56pub fn tx_to_lng(tx: f64, z: u8) -> f64 {
57 tx / (1u64 << z) as f64 * 360.0 - 180.0
58}
59
60pub fn ty_to_lat(ty: f64, z: u8) -> f64 {
61 (PI * (1.0 - 2.0 * ty / (1u64 << z) as f64)).sinh().atan().to_degrees()
62}
63
64impl TileId {
65 pub fn new(z: u8, x: u32, y: u32) -> Outcome<Self> {
66 if z > MAX_ZOOM {
67 return Err(err!("Zoom {} is past the deepest, {}.", z, MAX_ZOOM; Invalid, Input, Range));
68 }
69 let n = 1u64 << z;
70 if x as u64 >= n || y as u64 >= n {
71 return Err(err!("Tile {}/{}/{} is off a grid {} tiles across.", z, x, y, n;
72 Invalid, Input, Range));
73 }
74 Ok(Self { z, x, y })
75 }
76
77 /// The tile holding a position. Latitude is clamped to the map and longitude wrapped, so
78 /// a pole lands on the top or bottom row.
79 pub fn at(lat: f64, lng: f64, z: u8) -> Outcome<Self> {
80 if z > MAX_ZOOM {
81 return Err(err!("Zoom {} is past the deepest, {}.", z, MAX_ZOOM; Invalid, Input, Range));
82 }
83 if !(lat.is_finite() && lng.is_finite()) {
84 return Err(err!("{}, {} is not a position.", lat, lng; Invalid, Input));
85 }
86 let n = (1u64 << z) as f64;
87 let lat = lat.clamp(-crate::proj::WEB_MERCATOR_MAX_LAT, crate::proj::WEB_MERCATOR_MAX_LAT);
88 let tx = lng_to_tx(lng, z).rem_euclid(n);
89 let ty = lat_to_ty(lat, z);
90 let clamp = |t: f64| -> u32 { t.floor().clamp(0.0, n - 1.0) as u32 };
91 Ok(Self { z, x: clamp(tx), y: clamp(ty) })
92 }
93
94 /// The position at a fraction `u` across and `v` down the tile, as `(lat, lng)`.
95 pub fn point(&self, u: f64, v: f64) -> (f64, f64) {
96 (ty_to_lat(self.y as f64 + v, self.z), tx_to_lng(self.x as f64 + u, self.z))
97 }
98
99 /// The tile's edges: south, west, north and east, in degrees.
100 pub fn bounds(&self) -> (f64, f64, f64, f64) {
101 let (north, west) = self.point(0.0, 0.0);
102 let (south, east) = self.point(1.0, 1.0);
103 (south, west, north, east)
104 }
105
106 /// The tile containing this one at a zoom no deeper than its own.
107 pub fn ancestor(&self, z: u8) -> Outcome<Self> {
108 if z > self.z {
109 return Err(err!("A zoom-{} tile has no ancestor at zoom {}.", self.z, z; Invalid, Input, Range));
110 }
111 let d = self.z - z;
112 Ok(Self { z, x: self.x >> d, y: self.y >> d })
113 }
114
115 pub fn children(&self) -> Outcome<[Self; 4]> {
116 if self.z >= MAX_ZOOM {
117 return Err(err!("A zoom-{} tile has no children.", self.z; Invalid, Input, Range));
118 }
119 let (z, x, y) = (self.z + 1, self.x * 2, self.y * 2);
120 Ok([Self { z, x, y }, Self { z, x: x + 1, y }, Self { z, x, y: y + 1 }, Self { z, x: x + 1, y: y + 1 }])
121 }
122
123 /// Where this tile sits inside an ancestor, for drawing it from the ancestor's data past
124 /// an archive's deepest zoom: the fraction across and down the ancestor at which it begins,
125 /// and how many of it span the ancestor. `None` if `ancestor` does not contain it.
126 pub fn within(&self, ancestor: &TileId) -> Option<(f64, f64, f64)> {
127 if ancestor.z > self.z {
128 return None;
129 }
130 let d = self.z - ancestor.z;
131 if self.x >> d != ancestor.x || self.y >> d != ancestor.y {
132 return None;
133 }
134 let span = (1u64 << d) as f64;
135 let mask = (1u64 << d) - 1;
136 Some(((self.x as u64 & mask) as f64 / span, (self.y as u64 & mask) as f64 / span, span))
137 }
138}
139
140/// A 2D affine transform in canvas order: `x' = a x + c y + e`, `y' = b x + d y + f`.
141#[derive(Clone, Copy, Debug, PartialEq)]
142pub struct Affine {
143 pub a: f64,
144 pub b: f64,
145 pub c: f64,
146 pub d: f64,
147 pub e: f64,
148 pub f: f64,
149}
150
151impl Affine {
152 pub fn apply(&self, x: f64, y: f64) -> (f64, f64) {
153 (self.a * x + self.c * y + self.e, self.b * x + self.d * y + self.f)
154 }
155}
156
157impl Viewport {
158 /// The tiles of zoom `z` the screen shows, nearest the centre first, refusing more than
159 /// `max`.
160 ///
161 /// On the flat map these are the tiles under the screen, turned by its heading if it is
162 /// turned, wrapping east and west, each listed once however many copies of the world the
163 /// screen shows. On the globe they are the tiles whose bounding caps meet the view's
164 /// [`Viewport::bounding_cap`].
165 pub fn tiles_covering(&self, z: u8, max: usize) -> Outcome<Vec<TileId>> {
166 res!(self.check());
167 if z > MAX_ZOOM {
168 return Err(err!("Zoom {} is past the deepest, {}.", z, MAX_ZOOM; Invalid, Input, Range));
169 }
170 let n = 1u64 << z;
171 let nf = n as f64;
172 let mut out: Vec<(f64, TileId)> = Vec::new();
173 let mut seen = std::collections::HashSet::new();
174 match self.kind {
175 Projection::WebMercator => {
176 let f = self.frame();
177 // The screen's corners in fractional tiles, longitude unwrapped.
178 let corners = [(0.0, 0.0), (self.w, 0.0), (self.w, self.h), (0.0, self.h)];
179 let mut quad = [(0.0f64, 0.0f64); 4];
180 for (k, (sx, sy)) in corners.iter().enumerate() {
181 let (x, yd) = f.from_screen(*sx, *sy);
182 let lam = f.lam_0 + x / f.r_px;
183 let my = (-yd + f.y_0) / f.r_px;
184 quad[k] = ((lam + PI) / TAU * nf, (PI - my) / TAU * nf);
185 }
186 let (cx, cy) = {
187 let my = f.y_0 / f.r_px;
188 ((f.lam_0 + PI) / TAU * nf, (PI - my) / TAU * nf)
189 };
190 let xl = quad.iter().map(|p| p.0).fold(f64::MAX, f64::min).floor();
191 // A screen wider than the world needs each column once.
192 let xh = quad.iter().map(|p| p.0).fold(f64::MIN, f64::max).ceil().min(xl + nf);
193 let yl = quad.iter().map(|p| p.1).fold(f64::MAX, f64::min).floor().max(0.0);
194 let yh = quad.iter().map(|p| p.1).fold(f64::MIN, f64::max).ceil().min(nf);
195 let span = (xh - xl) * (yh - yl);
196 if span > (4 * max.max(1)) as f64 {
197 return Err(err!("The screen spans about {} zoom-{} tiles, more than {}.",
198 span as u64, z, max; Excessive, Size));
199 }
200 let turned = self.heading.rem_euclid(90.0) != 0.0;
201 let mut ty = yl;
202 while ty < yh {
203 let mut tx = xl;
204 while tx < xh {
205 let square = [(tx, ty), (tx + 1.0, ty), (tx + 1.0, ty + 1.0), (tx, ty + 1.0)];
206 if !turned || overlaps(&square, &quad) {
207 let t = TileId { z, x: (tx as i64).rem_euclid(n as i64) as u32, y: ty as u32 };
208 if seen.insert(t) {
209 let d = (tx + 0.5 - cx).hypot(ty + 0.5 - cy);
210 out.push((d, t));
211 if out.len() > max {
212 return Err(err!("The screen shows more than {} zoom-{} tiles.",
213 max, z; Excessive, Size));
214 }
215 }
216 }
217 tx += 1.0;
218 }
219 ty += 1.0;
220 }
221 },
222 Projection::Orthographic => {
223 let (c, r) = self.bounding_cap();
224 let lat_c = self.lat_0.clamp(-90.0, 90.0);
225 let lat_hi = (lat_c + r.to_degrees()).min(crate::proj::WEB_MERCATOR_MAX_LAT);
226 let lat_lo = (lat_c - r.to_degrees()).max(-crate::proj::WEB_MERCATOR_MAX_LAT);
227 let polar = lat_c + r.to_degrees() >= 90.0 || lat_c - r.to_degrees() <= -90.0;
228 let (tx_lo, tx_hi) = if polar || r >= FRAC_PI_2 {
229 (0.0, nf)
230 } else {
231 let half = (r.sin() / lat_c.to_radians().cos()).min(1.0).asin().to_degrees();
232 (lng_to_tx(self.lon_0 - half, z).floor(), lng_to_tx(self.lon_0 + half, z).ceil())
233 };
234 let (ty_lo, ty_hi) = (lat_to_ty(lat_hi, z).floor().max(0.0), lat_to_ty(lat_lo, z).ceil().min(nf));
235 let (cx, cy) = (lng_to_tx(self.lon_0, z), lat_to_ty(lat_c.clamp(-85.0, 85.0), z));
236 let span = (tx_hi - tx_lo).min(nf) * (ty_hi - ty_lo);
237 if span > (16 * max.max(1)) as f64 {
238 return Err(err!("The globe spans about {} zoom-{} tiles, more than {}.",
239 span as u64, z, max; Excessive, Size));
240 }
241 let mut ty = ty_lo;
242 while ty < ty_hi {
243 let mut tx = tx_lo;
244 while tx < tx_hi.min(tx_lo + nf) {
245 let t = TileId { z, x: (tx as i64).rem_euclid(n as i64) as u32, y: ty as u32 };
246 let (tc, tr) = tile_cap(&t);
247 let apart = (tc[0] * c[0] + tc[1] * c[1] + tc[2] * c[2]).clamp(-1.0, 1.0).acos();
248 if apart <= tr + r && seen.insert(t) {
249 let mut dx = (tx + 0.5 - cx).abs() % nf;
250 dx = dx.min(nf - dx);
251 out.push((dx.hypot(ty + 0.5 - cy), t));
252 if out.len() > max {
253 return Err(err!("The globe shows more than {} zoom-{} tiles.", max, z;
254 Excessive, Size));
255 }
256 }
257 tx += 1.0;
258 }
259 ty += 1.0;
260 }
261 },
262 }
263 out.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
264 Ok(out.into_iter().map(|(_, t)| t).collect())
265 }
266
267 /// The zoom, fractional, at which a tile `tile_px` pixels across is drawn one tile pixel
268 /// to one screen pixel at the centre. Both projections share their scale at the centre,
269 /// so the globe asks for the same tiles as the flat map with the same camera.
270 pub fn tile_zoom(&self, tile_px: f64) -> f64 {
271 let flat = Viewport { kind: Projection::WebMercator, ..*self };
272 if flat.check().is_err() || !(tile_px > 0.0) {
273 return 0.0;
274 }
275 let f = flat.frame();
276 (TAU * f.r_px / tile_px).log2()
277 }
278
279 /// The transform that paints a tile's bitmap, from its own coordinates -- `u` across and
280 /// `v` down, 0 to 1 -- onto the screen, with how far it strays from the true projection in
281 /// pixels.
282 ///
283 /// The transform is the flat map's for the same camera, the copy of the tile nearest the
284 /// centre. On the flat map it is exact. On the globe it is exact at the centre, where the
285 /// two projections share their scale, and the error is measured at the tile's corners,
286 /// edge midpoints and middle; `None` if any of those is behind the globe. At 32 degrees a
287 /// view 12 km high keeps the error under half a pixel on a phone, which is where a globe
288 /// can paint the same bitmaps as the map with nothing to hand off.
289 pub fn tile_affine(&self, t: &TileId) -> Option<(Affine, f64)> {
290 if self.check().is_err() {
291 return None;
292 }
293 let flat = Viewport { kind: Projection::WebMercator, ..*self };
294 let f = flat.frame();
295 let n = (1u64 << t.z) as f64;
296 let s = TAU * f.r_px / n;
297 // The tile's north-west corner on the unturned plane, on the copy nearest the centre.
298 let mut lam = t.x as f64 / n * TAU - PI - f.lam_0;
299 let mid = lam + PI / n;
300 lam -= TAU * ((mid + PI).div_euclid(TAU));
301 let x_nw = lam * f.r_px;
302 let y_nw = f.r_px * (PI - TAU * t.y as f64 / n) - f.y_0;
303 let aff = Affine {
304 a: s * f.ch,
305 b: -s * f.sh,
306 c: s * f.sh,
307 d: s * f.ch,
308 e: f.cx + x_nw * f.ch - y_nw * f.sh,
309 f: f.cy - x_nw * f.sh - y_nw * f.ch,
310 };
311 if self.kind == Projection::WebMercator {
312 return Some((aff, 0.0));
313 }
314 let mut err: f64 = 0.0;
315 for (u, v) in [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (0.0, 0.5), (0.5, 0.5), (1.0, 0.5),
316 (0.0, 1.0), (0.5, 1.0), (1.0, 1.0)]
317 {
318 let (lat, lng) = t.point(u, v);
319 let p = match self.forward(lat, lng) {
320 Some(p) => p,
321 None => return None,
322 };
323 let (x, y) = aff.apply(u, v);
324 err = err.max((p.x - x).hypot(p.y - y));
325 }
326 Some((aff, err))
327 }
328}
329
330/// A cap about a tile's middle holding the whole tile. A tile's north and south edges are
331/// parallels, not great circles, so the edge midpoints are measured as well as the corners,
332/// and a little is added for the bulge between them.
333fn tile_cap(t: &TileId) -> ([f64; 3], f64) {
334 let (lat, lng) = t.point(0.5, 0.5);
335 let c = unit_vec(lat, lng);
336 let mut far: f64 = 0.0;
337 for (u, v) in [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (0.0, 0.5), (1.0, 0.5), (0.0, 1.0),
338 (0.5, 1.0), (1.0, 1.0)]
339 {
340 let (a, b) = t.point(u, v);
341 let p = unit_vec(a, b);
342 far = far.max((p[0] * c[0] + p[1] * c[1] + p[2] * c[2]).clamp(-1.0, 1.0).acos());
343 }
344 (c, far * 1.02 + 1.0e-9)
345}
346
347/// Do two convex quadrilaterals overlap? Separating axes: the edge normals of both.
348fn overlaps(a: &[(f64, f64); 4], b: &[(f64, f64); 4]) -> bool {
349 for poly in [a, b] {
350 for k in 0..4 {
351 let (p, q) = (poly[k], poly[(k + 1) % 4]);
352 let (nx, ny) = (q.1 - p.1, p.0 - q.0);
353 let proj = |pts: &[(f64, f64); 4]| -> (f64, f64) {
354 let mut lo = f64::MAX;
355 let mut hi = f64::MIN;
356 for r in pts.iter() {
357 let d = r.0 * nx + r.1 * ny;
358 lo = lo.min(d);
359 hi = hi.max(d);
360 }
361 (lo, hi)
362 };
363 let (alo, ahi) = proj(a);
364 let (blo, bhi) = proj(b);
365 if ahi < blo || bhi < alo {
366 return false;
367 }
368 }
369 }
370 true
371}