gistools/proj/project/
nzmg.rs

1use crate::proj::{
2    Complex, CoordinateStep, Proj, ProjectCoordinates, TransformCoordinates, zpoly1, zpolyd1,
3};
4use alloc::rc::Rc;
5use core::cell::RefCell;
6use libm::fabs;
7
8// /******************************************************************************
9//  * Project:  PROJ.4
10//  * Purpose:  Implementation of the nzmg (New Zealand Map Grid) projection.
11//  *           Very loosely based upon DMA code by Bradford W. Drew
12//  * Author:   Gerald Evenden
13//  *
14//  ******************************************************************************
15//  * Copyright (c) 1995, Gerald Evenden
16//  *
17//  * Permission is hereby granted, free of charge, to any person obtaining a
18//  * copy of this software and associated documentation files (the "Software"),
19//  * to deal in the Software without restriction, including without limitation
20//  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
21//  * and/or sell copies of the Software, and to permit persons to whom the
22//  * Software is furnished to do so, subject to the following conditions:
23//  *
24//  * The above copyright notice and this permission notice shall be included
25//  * in all copies or substantial portions of the Software.
26//  *
27//  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
28//  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
29//  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
30//  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
31//  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
32//  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
33//  * DEALINGS IN THE SOFTWARE.
34//  *****************************************************************************/
35const EPSLN: f64 = 1e-10;
36const SEC5_TO_RAD: f64 = 0.484_813_681_109_536;
37const RAD_TO_SEC5: f64 = 2.062_648_062_470_963_8;
38const N_BF: usize = 5;
39const N_TPSI: usize = 9;
40const N_TPHI: usize = 8;
41const BF: [Complex; 6] = [
42    Complex { r: 0.7557853228, i: 0.0 },
43    Complex { r: 0.249204646, i: 0.003371507 },
44    Complex { r: -0.001541739, i: 0.041058560 },
45    Complex { r: -0.10162907, i: 0.01727609 },
46    Complex { r: -0.26623489, i: -0.36249218 },
47    Complex { r: -0.6870983, i: -1.1651967 },
48];
49
50/// # New Zealand Map Grid (EPSG:27200)
51///
52/// **Classification**: Custom grid-based projection
53///
54/// **Available forms**: Forward and inverse, spherical and ellipsoidal
55///
56/// **Defined area**: New Zealand
57///
58/// **Alias**: nzmg
59///
60/// **Domain**: 2D
61///
62/// **Input type**: Geodetic coordinates
63///
64/// **Output type**: Projected coordinates
65///
66/// ## Projection String
67/// ```ini
68/// +proj=nzmg
69/// ```
70///
71/// ## Required Parameters
72/// - None (all standard projection parameters are hard-coded for this projection)
73///
74/// ## Reference:
75/// - Department of Land and Survey Technical Circular 1973/32
76///   <http://www.linz.govt.nz/docs/miscellaneous/nz-map-definition.pdf>
77/// - OSG Technical Report 4.1
78///   <http://www.linz.govt.nz/docs/miscellaneous/nzmg.pdf>
79///
80/// ![New Zealand Map Grid (EPSG:27200)](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/nzmg.png?raw=true)
81#[derive(Debug, Clone, PartialEq)]
82pub struct NewZealandMapGridProjection {
83    proj: Rc<RefCell<Proj>>,
84}
85impl ProjectCoordinates for NewZealandMapGridProjection {
86    fn code(&self) -> i64 {
87        -1
88    }
89    fn name(&self) -> &'static str {
90        "New Zealand Map Grid"
91    }
92    fn names() -> &'static [&'static str] {
93        &["New Zealand Map Grid", "NewZealandMapGrid", "New_Zealand_Map_Grid", "nzmg"]
94    }
95}
96impl CoordinateStep for NewZealandMapGridProjection {
97    fn new(proj: Rc<RefCell<Proj>>) -> Self {
98        {
99            let proj = &mut proj.borrow_mut();
100            // force to International major axis
101            proj.a = 6378388.0;
102            proj.ra = 1. / proj.a;
103            proj.lam0 = 173.0_f64.to_radians();
104            proj.phi0 = -41.0_f64.to_radians();
105            proj.x0 = 2510000.;
106            proj.y0 = 6023150.;
107        }
108        NewZealandMapGridProjection { proj }
109    }
110    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
111        nzmg_e_forward(&self.proj.borrow(), p);
112    }
113    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
114        nzmg_e_inverse(&self.proj.borrow(), p);
115    }
116}
117
118/// New Zealand Map Grid Ellipsoidal forward project
119pub fn nzmg_e_forward<P: TransformCoordinates>(proj: &Proj, p: &mut P) {
120    let mut c_p = Complex::default();
121    let tpsi: [f64; 10] = [
122        0.6399175073,
123        -0.1358797613,
124        0.063294409,
125        -0.02526853,
126        0.0117879,
127        -0.0055161,
128        0.0026906,
129        -0.001333,
130        0.00067,
131        -0.00034,
132    ];
133
134    p.set_phi((p.phi() - proj.phi0) * RAD_TO_SEC5);
135    let mut i = N_TPSI;
136    while i > 0 {
137        i -= 1;
138        c_p.r = tpsi[i] + p.phi() * c_p.r;
139    }
140    c_p.r *= p.phi();
141    c_p.i = p.lam();
142    c_p = zpoly1(c_p, &BF, N_BF);
143    p.set_x(c_p.i);
144    p.set_y(c_p.r);
145}
146
147/// New Zealand Map Grid Ellipsoidal inverse project
148pub fn nzmg_e_inverse<P: TransformCoordinates>(proj: &Proj, p: &mut P) {
149    let mut c_p = Complex::default();
150    let mut fp = Complex::default();
151    let mut dp = Complex::default();
152
153    let tphi: [f64; 9] = [
154        1.5627014243,
155        0.5185406398,
156        -0.03333098,
157        -0.1052906,
158        -0.0368594,
159        0.007317,
160        0.01220,
161        0.00394,
162        -0.0013,
163    ];
164
165    c_p.r = p.y();
166    c_p.i = p.x();
167    let mut nn: usize = 20;
168    while nn > 0 {
169        nn -= 1;
170        let mut f = zpolyd1(c_p, &BF, N_BF, &mut fp);
171        f.r -= p.y();
172        f.i -= p.x();
173        let den = fp.r * fp.r + fp.i * fp.i;
174        dp.r = -(f.r * fp.r + f.i * fp.i) / den;
175        dp.i = -(f.i * fp.r - f.r * fp.i) / den;
176        c_p.r += dp.r;
177        c_p.i += dp.i;
178        if (fabs(dp.r) + fabs(dp.i)) <= EPSLN {
179            break;
180        }
181    }
182    if nn != 0 {
183        p.set_lam(c_p.i);
184        let mut phi = 0.0;
185        let mut i = N_TPHI;
186        while i > 0 {
187            i -= 1;
188            phi = tphi[i] + c_p.r * phi;
189        }
190        p.set_phi(proj.phi0 + c_p.r * phi * SEC5_TO_RAD);
191    } else {
192        p.set_lam(f64::MAX);
193        p.set_phi(f64::MAX);
194    }
195}
gistools/proj/project/nzmg.rs

gistools/proj/project/
nzmg.rs
