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
123
124
125
126
127
128
129

//! Albers Equal-Area Conic projection.

use crate::error::{ProjectionError, Result};
use crate::{to_degrees, to_radians};
use super::{ProjectionImpl, ProjectionParams};

pub(super) struct AlbersProj {
    lon0: f64,
    a: f64,
    e2: f64,
    n: f64,
    c: f64,
    rho0: f64,
    fe: f64,
    fn_: f64,
}

fn alpha(e2: f64, lat: f64) -> f64 {
    let e = e2.sqrt();
    let sin_lat = lat.sin();
    let esin = e * sin_lat;
    (1.0 - e2) * (
        sin_lat / (1.0 - e2 * sin_lat * sin_lat)
        - (1.0 / (2.0 * e)) * ((1.0 - esin) / (1.0 + esin)).ln()
    )
}

fn m(e2: f64, lat: f64) -> f64 {
    let sin_lat = lat.sin();
    lat.cos() / (1.0 - e2 * sin_lat * sin_lat).sqrt()
}

impl AlbersProj {
    pub fn new(p: &ProjectionParams, lat1_deg: f64, lat2_deg: f64) -> Result<Self> {
        let lat1 = to_radians(lat1_deg);
        let lat2 = to_radians(lat2_deg);
        let lat0 = to_radians(p.lat0);
        let a = p.ellipsoid.a;
        let e2 = p.ellipsoid.e2;

        let m1 = m(e2, lat1);
        let m2 = m(e2, lat2);
        let a1 = alpha(e2, lat1);
        let a2 = alpha(e2, lat2);
        let a0 = alpha(e2, lat0);

        let n = (m1 * m1 - m2 * m2) / (a2 - a1);
        let c = m1 * m1 + n * a1;
        let rho0 = a * (c - n * a0).sqrt() / n;

        Ok(AlbersProj {
            lon0: to_radians(p.lon0),
            a,
            e2,
            n,
            c,
            rho0,
            fe: p.false_easting,
            fn_: p.false_northing,
        })
    }
}

impl ProjectionImpl for AlbersProj {
    fn forward(&self, lon_deg: f64, lat_deg: f64) -> Result<(f64, f64)> {
        let lat = to_radians(lat_deg);
        let lon = to_radians(lon_deg);

        let a_val = alpha(self.e2, lat);
        let rho = self.a * (self.c - self.n * a_val).sqrt() / self.n;
        let theta = self.n * (lon - self.lon0);

        let x = rho * theta.sin() + self.fe;
        let y = self.rho0 - rho * theta.cos() + self.fn_;
        Ok((x, y))
    }

    fn inverse(&self, x: f64, y: f64) -> Result<(f64, f64)> {
        let mut x = x - self.fe;
        let mut y = self.rho0 - (y - self.fn_);

        let mut rho = (x * x + y * y).sqrt();
        if rho != 0.0 && self.n < 0.0 {
            rho = -rho;
            x = -x;
            y = -y;
        }

        let theta = x.atan2(y);
        let lon = theta / self.n + self.lon0;

        let q = (self.c - (rho * self.n / self.a).powi(2)) / self.n;
        let e2 = self.e2;
        let e = e2.sqrt();

        let phi = if e < 1e-15 {
            (q / 2.0).clamp(-1.0, 1.0).asin()
        } else {
            let mut phi = (q / 2.0).clamp(-1.0, 1.0).asin();
            let mut converged = false;
            for _ in 0..25 {
                let sin_phi = phi.sin();
                let cos_phi = phi.cos();
                if cos_phi.abs() < 1e-15 {
                    break;
                }
                let esin = e * sin_phi;
                let one_minus = 1.0 - e2 * sin_phi * sin_phi;
                let phi_new = phi + one_minus * one_minus / (2.0 * cos_phi) * (
                    q / (1.0 - e2)
                    - sin_phi / one_minus
                    + (1.0 / (2.0 * e)) * ((1.0 - esin) / (1.0 + esin)).ln()
                );
                if (phi_new - phi).abs() < 1e-12 {
                    phi = phi_new;
                    converged = true;
                    break;
                }
                phi = phi_new;
            }
            if !converged {
                return Err(ProjectionError::ConvergenceFailure { iterations: 25 });
            }
            phi
        };

        Ok((to_degrees(lon), to_degrees(phi)))
    }
}