gistools/proj/project/
ocea.rs

1use crate::proj::{
2    AZIMUTH_PROJECTION_CENTRE, CoordinateStep, LATITUDE_OF_FIRST_POINT, LATITUDE_OF_SECOND_POINT,
3    LONGITUDE_OF_FIRST_POINT, LONGITUDE_OF_PROJECTION_CENTRE, LONGITUDE_OF_SECOND_POINT, Proj,
4    ProjValue, ProjectCoordinates, TransformCoordinates,
5};
6use alloc::rc::Rc;
7use core::{
8    cell::RefCell,
9    f64::consts::{FRAC_PI_2, PI},
10};
11use libm::{asin, atan, atan2, cos, sin, sqrt, tan};
12
13/// Oblique Cylindrical Equal Area Variables
14#[derive(Debug, Default, Clone, PartialEq)]
15pub struct OceaData {
16    rok: f64,
17    rtk: f64,
18    sinphi: f64,
19    cosphi: f64,
20}
21
22/// Oblique Cylindrical Equal Area Projection
23#[derive(Debug, Clone, PartialEq)]
24pub struct ObliqueCylindricalEqualAreaProjection {
25    proj: Rc<RefCell<Proj>>,
26    store: RefCell<OceaData>,
27}
28impl ProjectCoordinates for ObliqueCylindricalEqualAreaProjection {
29    fn code(&self) -> i64 {
30        -1
31    }
32    fn name(&self) -> &'static str {
33        "Oblique Cylindrical Equal Area"
34    }
35    fn names() -> &'static [&'static str] {
36        &["Oblique Cylindrical Equal Area", "ocea"]
37    }
38}
39impl CoordinateStep for ObliqueCylindricalEqualAreaProjection {
40    fn new(proj: Rc<RefCell<Proj>>) -> Self {
41        let mut store = OceaData::default();
42        {
43            let proj = &mut proj.borrow_mut();
44
45            store.rok = 1. / proj.k0;
46            store.rtk = proj.k0;
47            let mut lam_p;
48            let phi_p;
49            // If the keyword "alpha" is found in the sentence then use 1point+1azimuth
50            if let Some(alpha) = proj.params.get(&AZIMUTH_PROJECTION_CENTRE) {
51                // Define Pole of oblique transformation from 1 point & 1 azimuth
52                // ERO: I've added PI so that the alpha is the angle from point 1 to
53                // point 2 from the North in a clockwise direction (to be consistent
54                // with omerc behavior)
55                let alpha = PI + alpha.f64();
56                let lonz = proj
57                    .params
58                    .get(&LONGITUDE_OF_PROJECTION_CENTRE)
59                    .unwrap_or(&ProjValue::default())
60                    .f64();
61                // Equation 9-8 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)
62                // Actually slightliy modified to use atan2(), as it is suggested by
63                // Snyder for equation 9-1, but this is not mentioned here
64                lam_p = atan2(-cos(alpha), -sin(proj.phi0) * sin(alpha)) + lonz;
65                // Equation 9-7 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)
66                phi_p = asin(cos(proj.phi0) * sin(alpha));
67                // If the keyword "alpha" is NOT found in the sentence then use 2points
68            } else {
69                // Define Pole of oblique transformation from 2 points
70                let phi_1 = proj
71                    .params
72                    .get(&LATITUDE_OF_FIRST_POINT)
73                    .unwrap_or(&ProjValue::default())
74                    .f64()
75                    .to_radians();
76                let phi_2 = proj
77                    .params
78                    .get(&LATITUDE_OF_SECOND_POINT)
79                    .unwrap_or(&ProjValue::default())
80                    .f64()
81                    .to_radians();
82                let lam_1 = proj
83                    .params
84                    .get(&LONGITUDE_OF_FIRST_POINT)
85                    .unwrap_or(&ProjValue::default())
86                    .f64()
87                    .to_radians();
88                let lam_2 = proj
89                    .params
90                    .get(&LONGITUDE_OF_SECOND_POINT)
91                    .unwrap_or(&ProjValue::default())
92                    .f64()
93                    .to_radians();
94                // Equation 9-1 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)
95                lam_p = atan2(
96                    cos(phi_1) * sin(phi_2) * cos(lam_1) - sin(phi_1) * cos(phi_2) * cos(lam_2),
97                    sin(phi_1) * cos(phi_2) * sin(lam_2) - cos(phi_1) * sin(phi_2) * sin(lam_1),
98                );
99
100                //  take care of proj.lam0 wrap-around when +lam_1=-90
101                if lam_1 == -FRAC_PI_2 {
102                    lam_p = -lam_p;
103                }
104
105                // Equation 9-2 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)
106                let cos_lamp_m_minus_lam_1 = cos(lam_p - lam_1);
107                let tan_phi_1 = tan(phi_1);
108                if tan_phi_1 == 0.0 {
109                    // Not sure if we want to support this case, but at least this
110                    // avoids a division by zero, and gives the same result as the below
111                    // atan()
112                    phi_p = if cos_lamp_m_minus_lam_1 >= 0.0 { -FRAC_PI_2 } else { FRAC_PI_2 };
113                } else {
114                    phi_p = atan(-cos_lamp_m_minus_lam_1 / tan_phi_1);
115                }
116            }
117            proj.lam0 = lam_p + FRAC_PI_2;
118            store.cosphi = cos(phi_p);
119            store.sinphi = sin(phi_p);
120            proj.es = 0.;
121        }
122        ObliqueCylindricalEqualAreaProjection { proj, store: store.into() }
123    }
124    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
125        ocea_s_forward(&self.store.borrow(), p);
126    }
127    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
128        ocea_s_inverse(&self.store.borrow(), p);
129    }
130}
131
132/// Oblated Equal Area Spheroidal forward project
133pub fn ocea_s_forward<P: TransformCoordinates>(ocae: &OceaData, p: &mut P) {
134    let mut y = sin(p.lam());
135    let t = cos(p.lam());
136    let mut x = atan((tan(p.phi()) * ocae.cosphi + ocae.sinphi * y) / t);
137    if t < 0. {
138        x += PI;
139    }
140    x *= ocae.rtk;
141    y = ocae.rok * (ocae.sinphi * sin(p.phi()) - ocae.cosphi * cos(p.phi()) * y);
142
143    p.set_x(x);
144    p.set_y(y);
145}
146
147/// Oblated Equal Area Spheroidal inverse project
148pub fn ocea_s_inverse<P: TransformCoordinates>(ocae: &OceaData, p: &mut P) {
149    let y = p.y() / ocae.rok;
150    let x = p.x() / ocae.rtk;
151    let t = sqrt(1. - y * y);
152    let s = sin(x);
153    p.set_phi(asin(y * ocae.sinphi + t * ocae.cosphi * s));
154    p.set_lam(atan2(t * ocae.sinphi * s - y * ocae.cosphi, t * cos(x)));
155}
gistools/proj/project/ocea.rs

gistools/proj/project/
ocea.rs
