gistools/readers/nadgrid/
mod.rs

1use crate::parsers::{FeatureReader, Reader};
2use alloc::{string::String, vec, vec::Vec};
3use libm::round;
4use s2json::{
5    MValue, Point, Properties, VectorFeature, VectorGeometry, VectorMultiPoint, VectorPoint,
6};
7
8/// Seconds to degrees (S / 3_600)
9const SEC2DEG: f64 = 0.00000484813681109536;
10
11/// A Subgrid contained inside a NadGrid
12#[derive(Default, Debug, Clone, PartialEq)]
13pub struct NadSubGrid {
14    /// CVS => lonLat coords
15    pub cvs: VectorMultiPoint,
16    /// ll => lower_lon_lat
17    pub ll: VectorPoint,
18    /// del => lon_lat_interval
19    pub del: VectorPoint,
20    /// lim => lon_lat_column_count
21    pub lim: VectorPoint,
22    /// count
23    pub count: i32,
24}
25
26/// A Subgrid contained inside a NadGrid
27#[derive(Default, Debug, Clone, PartialEq)]
28pub struct NadMetadata {
29    /// ll => lower_lon_lat
30    pub ll: Point,
31    /// del => lon_lat_interval
32    pub del: Point,
33    /// lim => lon_lat_column_count
34    pub lim: Point,
35    /// count
36    pub count: i32,
37}
38
39/// An LAS Shaped Vector Feature
40pub type NadVectorFeature = VectorFeature<NadMetadata, Properties, MValue>;
41
42/// A grid wrapper around a parsed .gsb file
43#[derive(Default, Debug, Clone)]
44pub struct NadGridDefinition<'a, T: Reader> {
45    /// Grid name
46    pub name: String,
47    /// If the grid is mandatory
48    pub mandatory: bool,
49    /// Grid data
50    pub grid: Option<&'a NadGridReader<T>>,
51    /// If the grid is null
52    pub is_null: bool,
53}
54
55/// The header of a NTv2 file
56#[derive(Default, Debug, Clone, PartialEq)]
57pub struct NadGridHeader {
58    /// Grid fields count
59    pub n_fields: i32,
60    /// Subgrid fields count
61    pub n_subgrid_fields: i32,
62    /// Subgrids count
63    pub n_subgrids: i32,
64    /// Shift type
65    pub shift_type: String,
66    /// from semi major axis
67    pub from_semi_major_axis: f64,
68    /// from semi minor axis
69    pub from_semi_minor_axis: f64,
70    /// to semi major axis
71    pub to_semi_major_axis: f64,
72    /// to semi minor axis
73    pub to_semi_minor_axis: f64,
74}
75
76/// Each subgrid has it's own data on how to decode the points inside the subgrid
77#[derive(Default, Debug, Clone, PartialEq)]
78pub struct NadSubGridHeader {
79    /// The name of the subgrid
80    pub name: String,
81    /// The name of the parent grid
82    pub parent: String,
83    /// The lower latitude of the subgrid
84    pub lower_latitude: f64,
85    /// The upper latitude of the subgrid
86    pub upper_latitude: f64,
87    /// The lower longitude of the subgrid
88    pub lower_longitude: f64,
89    /// The upper longitude of the subgrid
90    pub upper_longitude: f64,
91    /// The latitude interval
92    pub latitude_interval: f64,
93    /// The longitude interval
94    pub longitude_interval: f64,
95    /// The number of points in the subgrid
96    pub grid_node_count: i32,
97}
98
99/// A single Node describing how to decode the point
100#[derive(Default, Debug, Clone, PartialEq)]
101pub struct NadGridNode {
102    /// The latitude shift
103    pub latitude_shift: f64,
104    /// The longitude shift
105    pub longitude_shift: f64,
106    /// The latitude accuracy
107    pub latitude_accuracy: f64,
108    /// The longitude accuracy
109    pub longitude_accuracy: f64,
110}
111
112/// The metadata inside each vector feature
113#[derive(Default, Debug, Clone, PartialEq)]
114pub struct NadGridMetadata {
115    /// The lower longitude and latitude
116    pub lower_lon_lat: VectorPoint,
117    /// The longitude and latitude interval
118    pub lon_lat_interval: VectorPoint,
119    /// The number of longitude and latitude columns
120    pub lon_lat_column_count: VectorPoint,
121    /// The number of points
122    pub count: u64,
123}
124
125/// # NAD Grid Reader
126///
127/// ## Description
128/// Loads/reads a binary NTv2 file (.gsb)
129///
130/// Implements the [`FeatureReader`] trait
131///
132/// It should be noted that a proj4 Transformer usually uses this class internally. But if you want
133/// to manually parse a .gsb file, you can use this class directly.
134///
135/// ## Usage
136///
137/// The methods you have access to:
138/// - [`NadGridReader::new`]: Create a new NadGridReader
139/// - [`NadGridReader::len`]: Get the length of the feature count
140/// - [`NadGridReader::is_empty`]: Check if the reader is empty
141/// - [`NadGridReader::header`]: Read the header
142/// - [`NadGridReader::get_points`]: Read a subgrid into a Point
143/// - [`NadGridReader::get_feature`]: Read a subgrid into a Vector Feature
144/// - [`NadGridReader::iter`]: Create an iterator to collect the features
145/// - [`NadGridReader::par_iter`]: Create a parallel iterator to collect the features
146///
147/// ```rust
148/// use gistools::{parsers::{FeatureReader, FileReader}, readers::NadGridReader};
149/// use std::path::PathBuf;
150///
151/// let mut path = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
152/// path.push("tests/proj4/fixtures/BETA2007.gsb");
153///
154/// let nadgrid_reader = NadGridReader::new("test".into(), FileReader::from(path.clone()));
155/// let features = nadgrid_reader.iter().collect::<Vec<_>>();
156/// assert_eq!(features.len(), 1);
157/// ```
158///
159/// ## Links
160/// - <https://web.archive.org/web/20140127204822if_/http://www.mgs.gov.on.ca:80/stdprodconsume/groups/content/@mgs/@iandit/documents/resourcelist/stel02_047447.pdf>
161/// - <http://mimaka.com/help/gs/html/004_NTV2%20Data%20Format.htm>
162#[derive(Debug, Clone)]
163pub struct NadGridReader<T: Reader> {
164    /// The name of the grid
165    pub key: String,
166    reader: T,
167    is_little_endian: bool,
168    header: NadGridHeader,
169    subgrids: Vec<NadSubGrid>,
170}
171impl<T: Reader> NadGridReader<T> {
172    /// Create a new NadGridReader
173    pub fn new(key: String, reader: T) -> Self {
174        let mut nad_grid = Self {
175            key,
176            reader,
177            is_little_endian: false,
178            header: NadGridHeader::default(),
179            subgrids: vec![],
180        };
181        nad_grid.detect_little_endian();
182        nad_grid.read_header();
183        nad_grid.read_sub_grids();
184
185        nad_grid
186    }
187
188    /// Get the length of the feature count
189    pub fn len(&self) -> u64 {
190        self.subgrids.len() as u64
191    }
192
193    /// Check if the reader is empty
194    pub fn is_empty(&self) -> bool {
195        self.len() == 0
196    }
197
198    /// read the header
199    pub fn header(&self) -> &NadGridHeader {
200        &self.header
201    }
202
203    /// Read a subgrid into a Point
204    pub fn get_points(&self, index: u64) -> Option<VectorMultiPoint> {
205        if let Some(grid) = self.subgrids.get(index as usize) {
206            return Some(grid.cvs.clone());
207        }
208        None
209    }
210
211    /// Read a subgrid into a Vector Feature
212    pub fn get_feature(&self, index: u64) -> Option<NadVectorFeature> {
213        if let Some(grid) = self.subgrids.get(index as usize) {
214            let NadSubGrid { cvs, ll, del, lim, count, .. } = grid;
215            return Some(VectorFeature::new_wm(
216                None,
217                Properties::default(),
218                VectorGeometry::new_multipoint(cvs.clone(), None),
219                Some(NadMetadata {
220                    ll: Point(ll.x, ll.y),
221                    del: Point(del.x, del.y),
222                    lim: Point(lim.x, lim.y),
223                    count: *count,
224                }),
225            ));
226        }
227        None
228    }
229
230    // INTERNAL
231
232    /// Set the little endian flag
233    fn detect_little_endian(&mut self) {
234        let NadGridReader { reader, .. } = self;
235        let mut n_fields = reader.int32_be(Some(8));
236        if n_fields == 11 {
237            return;
238        }
239        n_fields = reader.int32_le(Some(8));
240        if n_fields != 11 {
241            panic!("Failed to detect nadgrid endian-ness, defaulting to little-endian");
242        }
243        self.is_little_endian = true;
244    }
245
246    /// Parse the main header data. Describes the subgrids to decode lon-lat
247    fn read_header(&mut self) {
248        let NadGridReader { reader, is_little_endian, .. } = self;
249        let le = *is_little_endian;
250        self.header.n_fields = reader.int32(Some(8), Some(le));
251        self.header.n_subgrid_fields = reader.int32(Some(24), Some(le));
252        self.header.n_subgrids = reader.int32(Some(40), Some(le));
253        self.header.shift_type = reader.parse_string(Some(56), Some(8));
254        self.header.from_semi_major_axis = reader.f64(Some(120), Some(le));
255        self.header.from_semi_minor_axis = reader.f64(Some(136), Some(le));
256        self.header.to_semi_major_axis = reader.f64(Some(152), Some(le));
257        self.header.to_semi_minor_axis = reader.f64(Some(168), Some(le));
258    }
259
260    /// Build all grid data
261    fn read_sub_grids(&mut self) {
262        let mut grid_offset = 176;
263        let mut i = 0;
264        while i < self.header.n_subgrids {
265            let sub_header = self.read_sub_grid_header(grid_offset);
266            let nodes = self.read_grid_nodes(grid_offset, &sub_header);
267            let lon_column_count = round(
268                1.0 + (sub_header.upper_longitude - sub_header.lower_longitude)
269                    / sub_header.longitude_interval,
270            );
271            let lat_column_count = round(
272                1.0 + (sub_header.upper_latitude - sub_header.lower_latitude)
273                    / sub_header.latitude_interval,
274            );
275
276            self.subgrids.push(NadSubGrid {
277                cvs: nodes
278                    .iter()
279                    .map(|node| {
280                        let NadGridNode { longitude_shift, latitude_shift, .. } = node;
281                        VectorPoint::new_xy(
282                            longitude_shift * SEC2DEG,
283                            latitude_shift * SEC2DEG,
284                            None,
285                        )
286                    })
287                    .collect(),
288                ll: VectorPoint::new_xy(
289                    sub_header.lower_longitude * SEC2DEG,
290                    sub_header.lower_latitude * SEC2DEG,
291                    None,
292                ),
293                del: VectorPoint::new_xy(
294                    sub_header.longitude_interval * SEC2DEG,
295                    sub_header.latitude_interval * SEC2DEG,
296                    None,
297                ),
298                lim: VectorPoint::new_xy(lon_column_count, lat_column_count, None),
299                count: sub_header.grid_node_count,
300            });
301            grid_offset += 176 + sub_header.grid_node_count as u64 * 16;
302
303            i += 1;
304        }
305    }
306
307    /// Read a subgrid header
308    ///
309    /// ## Parameters
310    /// - `offset`: offset to read in the subgrid header
311    ///
312    /// ## Returns
313    /// The subgrid header
314    fn read_sub_grid_header(&self, offset: u64) -> NadSubGridHeader {
315        let NadGridReader { reader, is_little_endian, .. } = self;
316        let le = *is_little_endian;
317        NadSubGridHeader {
318            name: reader.parse_string(Some(offset + 8), Some(8)),
319            parent: reader.parse_string(Some(offset + 24), Some(8)),
320            lower_latitude: reader.f64(Some(offset + 72), Some(le)),
321            upper_latitude: reader.f64(Some(offset + 88), Some(le)),
322            lower_longitude: reader.f64(Some(offset + 104), Some(le)),
323            upper_longitude: reader.f64(Some(offset + 120), Some(le)),
324            latitude_interval: reader.f64(Some(offset + 136), Some(le)),
325            longitude_interval: reader.f64(Some(offset + 152), Some(le)),
326            grid_node_count: reader.int32(Some(offset + 168), Some(le)),
327        }
328    }
329
330    /// Read the grid nodes
331    ///
332    /// ## Parameters
333    /// - `offset`: offset of the grid
334    /// - `grid_header`: header of the grid
335    ///
336    /// ## Returns
337    /// An array of grid nodes
338    fn read_grid_nodes(&self, offset: u64, grid_header: &NadSubGridHeader) -> Vec<NadGridNode> {
339        let NadGridReader { reader, is_little_endian, .. } = self;
340        let le = *is_little_endian;
341        let node_count = grid_header.grid_node_count as u64;
342        let nodes_offset = offset + 176;
343        let grl: u64 = 16; // grid_record_length
344        let mut grid_shift_records = vec![];
345        let mut i: u64 = 0;
346        while i < node_count {
347            grid_shift_records.push(NadGridNode {
348                latitude_shift: reader.f32(Some(nodes_offset + i * grl), Some(le)) as f64,
349                longitude_shift: reader.f32(Some(nodes_offset + i * grl + 4), Some(le)) as f64,
350                latitude_accuracy: reader.f32(Some(nodes_offset + i * grl + 8), Some(le)) as f64,
351                longitude_accuracy: reader.f32(Some(nodes_offset + i * grl + 12), Some(le)) as f64,
352            });
353            i += 1;
354        }
355        grid_shift_records
356    }
357}
358
359/// The NadGrid Iterator tool
360#[derive(Debug)]
361pub struct NadGridIterator<'a, T: Reader> {
362    reader: &'a NadGridReader<T>,
363    index: u64,
364    len: u64,
365}
366impl<T: Reader> Iterator for NadGridIterator<'_, T> {
367    type Item = NadVectorFeature;
368
369    fn next(&mut self) -> Option<Self::Item> {
370        let cdf_reader = &self.reader;
371        if self.index >= self.len {
372            None
373        } else if let Some(point) = cdf_reader.get_feature(self.index) {
374            self.index += 1;
375            Some(point)
376        } else {
377            None
378        }
379    }
380}
381/// A feature reader trait with a callback-based approach
382impl<T: Reader> FeatureReader<NadMetadata, Properties, MValue> for NadGridReader<T> {
383    type FeatureIterator<'a>
384        = NadGridIterator<'a, T>
385    where
386        T: 'a;
387
388    fn iter(&self) -> Self::FeatureIterator<'_> {
389        NadGridIterator { reader: self, index: 0, len: self.len() }
390    }
391
392    fn par_iter(&self, pool_size: usize, thread_id: usize) -> Self::FeatureIterator<'_> {
393        let pool_size = pool_size as u64;
394        let thread_id = thread_id as u64;
395        let start = self.len() * thread_id / pool_size;
396        let end = self.len() * (thread_id + 1) / pool_size;
397        NadGridIterator { reader: self, index: start, len: end }
398    }
399}
gistools/readers/nadgrid/mod.rs

gistools/readers/nadgrid/
mod.rs
