wbprojection 0.1.0

Whitebox Projections is a map projection library for Rust, inspired by PROJ
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

//! Winkel Tripel projection.

use super::{ProjectionImpl, ProjectionParams};
use crate::error::{ProjectionError, Result};
use crate::{to_degrees, to_radians};
use std::f64::consts::{FRAC_PI_2, PI};

pub(super) struct WinkelTripelProj {
    lon0: f64,
    a: f64,
    fe: f64,
    fn_: f64,
}

impl WinkelTripelProj {
    pub fn new(p: &ProjectionParams) -> Result<Self> {
        Ok(WinkelTripelProj {
            lon0: to_radians(p.lon0),
            a: p.ellipsoid.a,
            fe: p.false_easting,
            fn_: p.false_northing,
        })
    }

    fn wrap_lon(mut lon: f64) -> f64 {
        while lon > PI {
            lon -= 2.0 * PI;
        }
        while lon < -PI {
            lon += 2.0 * PI;
        }
        lon
    }

    fn sinc(alpha: f64) -> f64 {
        if alpha.abs() < 1e-12 {
            1.0 - alpha * alpha / 6.0
        } else {
            alpha.sin() / alpha
        }
    }

    fn normalized_forward(lon_rel: f64, lat: f64) -> (f64, f64) {
        // Winkel Tripel uses standard parallel φ1 = arccos(2/π)
        let phi1 = (2.0 / PI).acos();
        let cos_phi1 = phi1.cos();

        let cos_lat = lat.cos();
        let sin_lat = lat.sin();
        let half_lon = lon_rel * 0.5;
        let alpha = (cos_lat * half_lon.cos()).acos();
        let s = Self::sinc(alpha);

        let x_aitoff = 2.0 * cos_lat * half_lon.sin() / s;
        let y_aitoff = sin_lat / s;

        let x_equi = lon_rel * cos_phi1;
        let y_equi = lat;

        (0.5 * (x_aitoff + x_equi), 0.5 * (y_aitoff + y_equi))
    }
}

impl ProjectionImpl for WinkelTripelProj {
    fn forward(&self, lon_deg: f64, lat_deg: f64) -> Result<(f64, f64)> {
        let lon = to_radians(lon_deg);
        let lat = to_radians(lat_deg);
        if !(-FRAC_PI_2..=FRAC_PI_2).contains(&lat) {
            return Err(ProjectionError::out_of_bounds(
                "latitude outside valid range [-90, 90]",
            ));
        }

        let lon_rel = Self::wrap_lon(lon - self.lon0);
        let (x_norm, y_norm) = Self::normalized_forward(lon_rel, lat);
        Ok((self.a * x_norm + self.fe, self.a * y_norm + self.fn_))
    }

    fn inverse(&self, x: f64, y: f64) -> Result<(f64, f64)> {
        let x_target = (x - self.fe) / self.a;
        let y_target = (y - self.fn_) / self.a;

        let mut lat = y_target.clamp(-FRAC_PI_2, FRAC_PI_2);
        let mut lon_rel = (x_target * 0.8).clamp(-PI, PI);

        for _ in 0..60 {
            let (fx, fy) = Self::normalized_forward(lon_rel, lat);
            let r1 = fx - x_target;
            let r2 = fy - y_target;

            if r1.abs() < 1e-13 && r2.abs() < 1e-13 {
                let lon = Self::wrap_lon(self.lon0 + lon_rel);
                return Ok((to_degrees(lon), to_degrees(lat)));
            }

            let h = 1e-7;
            let lat_h = (lat + h).clamp(-FRAC_PI_2, FRAC_PI_2);
            let lon_h = Self::wrap_lon(lon_rel + h);

            let (fx_lon_h, fy_lon_h) = Self::normalized_forward(lon_h, lat);
            let (fx_lat_h, fy_lat_h) = Self::normalized_forward(lon_rel, lat_h);

            let j11 = (fx_lon_h - fx) / h;
            let j12 = (fx_lat_h - fx) / h;
            let j21 = (fy_lon_h - fy) / h;
            let j22 = (fy_lat_h - fy) / h;

            let det = j11 * j22 - j12 * j21;
            if det.abs() < 1e-16 {
                break;
            }

            let d_lon = (-r1 * j22 + r2 * j12) / det;
            let d_lat = (-j11 * r2 + j21 * r1) / det;

            lon_rel = Self::wrap_lon(lon_rel + d_lon.clamp(-0.5, 0.5));
            lat = (lat + d_lat.clamp(-0.5, 0.5)).clamp(-FRAC_PI_2, FRAC_PI_2);
        }

        Err(ProjectionError::ConvergenceFailure { iterations: 60 })
    }
}