proj4rs/nadgrids/files/
ntv2.rs

1//!
2//! Nadgrid parser
3//!
4use crate::errors::{Error, Result};
5use crate::log::trace;
6use crate::math::consts::SEC_TO_RAD;
7use crate::nadgrids::grid::{REL_TOLERANCE_HGRIDSHIFT, Grid, GridId, Lp};
8use crate::nadgrids::header::error_str::*;
9use crate::nadgrids::header::{Endianness, Header};
10use crate::nadgrids::Catalog;
11use std::io::Read;
12
13const NTV2_HEADER_SIZE: usize = 11 * 16;
14
15/// Ntv2 reader
16pub(super) fn read_ntv2<R: Read>(catalog: &Catalog, key: &str, read: &mut R) -> Result<()> {
17    let mut head = Header::<NTV2_HEADER_SIZE>::new();
18
19    trace!("Reading  ntv2 {}", key);
20
21    // Read overview header
22    head.read(read)?;
23    // Check endianness
24    head.endian = if head.get_u8(8) == 11 {
25        Endianness::native()
26    } else {
27        Endianness::other()
28    };
29
30    let nsubgrids = head.get_u32(40) as usize;
31
32    trace!("Reading ntv2 {} subgrids {}", key, nsubgrids);
33
34    // Read subsequent grids
35    (0..nsubgrids).try_for_each(|_| read_ntv2_grid(catalog, key, head.read(read)?, read))
36}
37
38/// Read ntv2 grid data
39fn read_ntv2_grid<R: Read>(
40    catalog: &Catalog,
41    key: &str,
42    head: &Header<NTV2_HEADER_SIZE>,
43    read: &mut R,
44) -> Result<()> {
45    match head.get_str(0, 8) {
46        Ok("SUB_NAME") => Ok(()),
47        _ => Err(Error::InvalidNtv2GridFormat(ERR_INVALID_HEADER)),
48    }?;
49
50    let id = head.get_id(8);
51    let mut lineage = head.get_id(24);
52    if lineage.as_str().trim() == "NONE" {
53        lineage = GridId::root();
54    }
55
56    let mut ll = Lp {
57        lam: -head.get_f64(120), // W_LONG
58        phi: head.get_f64(72),   // S_LAT
59    };
60
61    let mut ur = Lp {
62        lam: -head.get_f64(104), // E_LONG
63        phi: head.get_f64(88),   // N_LAT
64    };
65
66    let mut del = Lp {
67        lam: head.get_f64(152), // longitude interval
68        phi: head.get_f64(136), // latitude interval
69    };
70
71    let lim = Lp {
72        lam: (((ur.lam - ll.lam).abs() / del.lam + 0.5) + 1.).floor(),
73        phi: (((ur.phi - ll.phi).abs() / del.phi + 0.5) + 1.).floor(),
74    };
75
76    // units are in seconds of degree.
77    ll.lam *= SEC_TO_RAD;
78    ll.phi *= SEC_TO_RAD;
79    ur.lam *= SEC_TO_RAD;
80    ur.phi *= SEC_TO_RAD;
81    del.lam *= SEC_TO_RAD;
82    del.phi *= SEC_TO_RAD;
83
84    // Read matrix data
85    let nrows = lim.phi as usize;
86    let rowsize = lim.lam as usize;
87
88    let gs_count = head.get_u32(168) as usize;
89    if gs_count != nrows * rowsize {
90        return Err(Error::InvalidNtv2GridFormat(ERR_GSCOUNT_NOT_MATCHING));
91    }
92
93    // Read grid data
94    trace!(
95        "Reading  data for grid {}:{}:{}",
96        key,
97        id.as_str(),
98        lineage.as_str()
99    );
100
101    let mut buf = head.rebind::<16>();
102    let mut cvs: Vec<Lp> = (0..gs_count)
103        .map(|_| {
104            buf.read(read)?;
105            // NOTE: phi and lam are inverted
106            Ok(Lp {
107                phi: SEC_TO_RAD * (buf.get_f32(0) as f64),
108                lam: SEC_TO_RAD * (buf.get_f32(4) as f64) * -1.0, // NOTE: Compensate NT convention
109            })
110        })
111        .collect::<Result<Vec<_>>>()?;
112
113    // See https://geodesie.ign.fr/contenu/fichiers/documentation/algorithmes/notice/NT111_V1_HARMEL_TransfoNTF-RGF93_FormatGrilleNTV2.pdf
114
115    // In proj4, rows are stored in reverse order
116    for i in 0..nrows {
117        let offs = i * rowsize;
118        cvs[offs..(offs + rowsize)].reverse();
119    }
120
121    let epsilon = (del.lam.abs() + del.phi.abs()) * REL_TOLERANCE_HGRIDSHIFT;
122
123    catalog.add_grid(
124        key.into(),
125        Grid {
126            id,
127            lineage,
128            ll,
129            ur,
130            del,
131            lim,
132            epsilon,
133            cvs: cvs.into_boxed_slice(),
134        },
135    )
136}
137
138#[cfg(test)]
139mod tests {
140    use super::*;
141    use crate::nadgrids::Catalog;
142    use crate::tests::setup;
143    use std::env;
144    use std::fs::File;
145    use std::io::BufReader;
146    use std::path::Path;
147
148    macro_rules! fixture {
149        ($name:expr) => {
150            Path::new(&env::var("CARGO_MANIFEST_DIR").unwrap())
151                .join("fixtures")
152                .as_path()
153                .join($name)
154                .as_path()
155        };
156    }
157
158    macro_rules! load_ntv2 {
159        ($cat:expr, $name:expr) => {
160            // Use a BufReader or efficiency
161            let file = File::open(fixture!($name)).unwrap();
162            let mut read = BufReader::new(file);
163            read_ntv2($cat, $name, &mut read).unwrap();
164        };
165    }
166
167    #[test]
168    fn ntv2_100800401_gsb() {
169        setup();
170
171        let catalog = Catalog::default();
172        load_ntv2!(&catalog, "100800401.gsb");
173
174        let grids = catalog.find("100800401.gsb").unwrap().collect::<Vec<_>>();
175        assert_eq!(grids.len(), 1);
176
177        let grid = grids[0];
178        assert!(grid.is_root());
179        assert_eq!(grid.id.as_str(), "0INT2GRS");
180        assert_eq!(grid.cvs.len(), 1591);
181    }
182
183    #[test]
184    #[cfg(feature = "local_tests")]
185    fn ntv2_bwta2017_gsb() {
186        setup();
187
188        let catalog = Catalog::default();
189        load_ntv2!(&catalog, "BWTA2017.gsb");
190
191        let grids = catalog.find("BWTA2017.gsb").unwrap().collect::<Vec<_>>();
192        assert_eq!(grids.len(), 1);
193
194        let grid = grids[0];
195        assert!(grid.is_root());
196        assert_eq!(grid.id.as_str(), "DHDN90  ");
197        assert_eq!(grid.cvs.len(), 24514459);
198    }
199}
proj4rs/nadgrids/files/ntv2.rs

proj4rs/nadgrids/files/
ntv2.rs
