diff options
| -rw-r--r-- | hittekaart/src/gpx.rs | 27 |
1 files changed, 22 insertions, 5 deletions
diff --git a/hittekaart/src/gpx.rs b/hittekaart/src/gpx.rs index 68bd6fb..86cb88b 100644 --- a/hittekaart/src/gpx.rs +++ b/hittekaart/src/gpx.rs @@ -38,15 +38,32 @@ impl Coordinates { /// projection](https://en.wikipedia.org/wiki/Web_Mercator_projection) of the coordinates. /// /// Returns the `(x, y)` projection, where both are in the range `[0, 256 * 2^zoom)`. + /// + /// Note that coordinates outside the range are silently truncated to the limits. pub fn web_mercator(self, zoom: u32) -> (u64, u64) { const WIDTH: f64 = super::layer::TILE_WIDTH as f64; const HEIGHT: f64 = super::layer::TILE_HEIGHT as f64; - let lambda = self.longitude.to_radians(); - let phi = self.latitude.to_radians(); - let x = 2u64.pow(zoom) as f64 / (2.0 * PI) * WIDTH * (lambda + PI); - let y = - 2u64.pow(zoom) as f64 / (2.0 * PI) * HEIGHT * (PI - (PI / 4.0 + phi / 2.0).tan().ln()); + // We support longitudes from 180° W/E, which is represented as PI in radians. + // For latitudes, we actually don't go all the way up to 90° N/S, see + // https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames + // To quote: + // This value will fall in the range (-π, π) for latitudes between 85.0511 °S and + // 85.0511 °N. + // For the curious, the number 85.0511 is the result of arctan(sinh(π)). By using + // this bound, the entire map becomes a (very large) square. + // We compute the exact bound here. + let max_lon: f64 = PI; + let max_lat: f64 = PI.sinh().atan(); + + let xmax = 2u64.pow(zoom) as f64 * WIDTH; + let ymax = 2u64.pow(zoom) as f64 * HEIGHT; + + let lambda = self.longitude.to_radians().clamp(-max_lon, max_lon); + let phi = self.latitude.to_radians().clamp(-max_lat, max_lat); + let phi = (phi.tan() + (1.0 / phi.cos())).ln(); + let x = xmax / (2.0 * PI) * (lambda + PI); + let y = ymax * (1.0 - phi / PI) / 2.0; (x.floor() as u64, y.floor() as u64) } } |