gistools/readers/las/
reader.rs

1use super::{LASExtendedVariableLengthRecord, LASHeader, LASPoint};
2use crate::{
3    parsers::{FeatureReader, Reader},
4    proj::Transformer,
5    readers::{
6        FieldTagNames, GeoStore, GeoTIFFTypes, build_transform_from_geo_keys,
7        parse_geotiff_raw_geokeys,
8    },
9};
10use alloc::{collections::BTreeMap, string::String, vec::Vec};
11use s2json::{Properties, VectorFeature, VectorGeometry, VectorPoint};
12
13/// Options for the LAS Reader
14#[derive(Debug, Default, Clone)]
15pub struct LASReaderOptions {
16    /// Whether to transform the data to WGS84 if it's not already in WGS84
17    pub dont_transform: bool,
18    /// List of EPSG codes to utilize e.g. `{ "4326": "WKT_STRING" }``
19    pub epsg_codes: BTreeMap<String, String>,
20}
21
22/// An LAS Shaped Vector Feature
23pub type LASVectorFeature = VectorFeature<(), Properties, LASPoint>;
24
25/// # LAS Reader
26///
27/// ## Description
28/// Reads LAS data. Supports up to the LAS 1.4 specification.
29/// [See specification](https://www.asprs.org/wp-content/uploads/2010/12/LAS_1_4_r13.pdf)
30///
31/// Implements the [`FeatureReader`] trait
32///
33/// Data is stored like so:
34/// ```txt
35/// |           PUBLIC HEADER BLOCK           |
36/// |         VARIABLE LENGTH RECORDS         |
37/// |            POINT DATA RECORDS           |
38/// ```
39///
40/// ## Usage
41///
42/// The methods you have access to:
43/// - [`LASReader::new`]: Create a new LASReader
44/// - [`LASReader::len`]: Get the number of points stored
45/// - [`LASReader::is_empty`]: Check if the reader is empty
46/// - [`LASReader::get_feature`]: Reads a point in at index as a feature
47/// - [`LASReader::get_point`]: Reads a point in at index
48/// - [`LASReader::iter`]: Create an iterator to collect the features
49/// - [`LASReader::par_iter`]: Create a parallel iterator to collect the features
50///
51/// ```rust
52/// use gistools::{parsers::{FeatureReader, FileReader}, readers::LASReader};
53/// use std::path::PathBuf;
54///
55/// let mut path = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
56/// path.push("tests/readers/las/fixtures/simple1_1.las");
57///
58/// let las_reader = LASReader::new(FileReader::from(path.clone()), None);
59/// let features = las_reader.iter().collect::<Vec<_>>();
60/// assert_eq!(features.len(), 1_065);
61/// ```
62///
63/// ## Links
64/// - <https://www.usgs.gov/ngp-standards-and-specifications/lidar-base-specification-online>
65/// - <https://www.asprs.org/wp-content/uploads/2010/12/LAS_1_4_r13.pdf>
66/// - <https://liblas.org/development/index.html>
67/// - <https://downloads.rapidlasso.de/doc/LAZ_Specification_1.4_R1.pdf>
68/// - <https://github.com/PDAL/PDAL>
69/// - <https://github.com/libLAS/libLAS> (deprecated for PDAL)
70/// - <https://github.com/LASzip>
71#[derive(Debug, Clone)]
72pub struct LASReader<T: Reader> {
73    reader: T,
74    /// Public Header Block
75    pub header: LASHeader,
76    /// Extended VARIABLE LENGTH RECORDS
77    pub variable_length_records: BTreeMap<u32, LASExtendedVariableLengthRecord>,
78    /// WKT projection string
79    pub wkt: Option<String>,
80    /// GeoKeyDirectory
81    pub geo_key_directory: GeoStore,
82    transformer: Transformer,
83    dont_transform: bool,
84}
85impl<T: Reader> LASReader<T> {
86    /// Create a new LASReader
87    pub fn new(reader: T, options: Option<LASReaderOptions>) -> Self {
88        let options = options.unwrap_or_default();
89        let header = LASHeader::from_reader(&reader);
90        let variable_length_records = las_parse_variable_length_records(&header, &reader);
91        let mut transformer = Transformer::new();
92        for (epsg_code, wkt) in options.epsg_codes.iter() {
93            transformer.insert_epsg_code(epsg_code.clone(), wkt.clone());
94        }
95        let wkt = build_wkt(&header, &variable_length_records, &mut transformer);
96        let geo_key_directory =
97            build_geo_key_directory(&variable_length_records, &mut transformer, wkt.is_none());
98        Self {
99            reader,
100            header,
101            variable_length_records,
102            wkt,
103            geo_key_directory,
104            transformer,
105            dont_transform: options.dont_transform,
106        }
107    }
108
109    /// Get the number of points stored
110    pub fn len(&self) -> u64 {
111        self.header.num_points as u64
112    }
113
114    /// Check if the reader is empty
115    pub fn is_empty(&self) -> bool {
116        self.len() == 0
117    }
118
119    /// Reads a point in at index as a feature
120    pub fn get_feature(&self, index: u64) -> Option<LASVectorFeature> {
121        self.get_point(index).map(|point| {
122            VectorFeature::new_wm(
123                None,
124                Properties::default(),
125                VectorGeometry::new_point(point, None),
126                None,
127            )
128        })
129    }
130
131    /// Reads a point in at index
132    pub fn get_point(&self, index: u64) -> Option<VectorPoint<LASPoint>> {
133        let Self { reader, header, dont_transform, .. } = self;
134        let LASHeader {
135            num_points,
136            offset_to_points,
137            point_data_format_id: format,
138            point_data_record_length,
139            ..
140        } = header;
141        if index + 1 > *num_points as u64 {
142            return None;
143        }
144        let offset_to_points = *offset_to_points as u64;
145        let point_data_record_length = *point_data_record_length as u64;
146        let offset = offset_to_points + index * point_data_record_length;
147        let point = if *format == 0 {
148            LASPoint::format0(reader, offset)
149        } else if *format == 1 {
150            LASPoint::format1(reader, offset)
151        } else if *format == 2 {
152            LASPoint::format2(reader, offset)
153        } else if *format == 3 {
154            LASPoint::format3(reader, offset)
155        } else if *format == 4 {
156            LASPoint::format4(reader, offset)
157        } else if *format == 5 {
158            LASPoint::format5(reader, offset)
159        } else if *format == 6 {
160            LASPoint::format6(reader, offset)
161        } else if *format == 7 {
162            LASPoint::format7(reader, offset)
163        } else if *format == 8 {
164            LASPoint::format8(reader, offset)
165        } else if *format == 9 {
166            LASPoint::format9(reader, offset)
167        } else if *format == 10 {
168            LASPoint::format10(reader, offset)
169        } else {
170            panic!("Unknown Point Data Format ID: {}", format);
171        };
172        let mut vp = point.to_vector_point(header);
173
174        if !*dont_transform {
175            self.transformer.forward_mut(&mut vp);
176        }
177
178        Some(vp)
179    }
180}
181
182/// The LAS Iterator tool
183#[derive(Debug)]
184pub struct LASIterator<'a, T: Reader> {
185    reader: &'a LASReader<T>,
186    index: u64,
187    len: u64,
188}
189impl<T: Reader> Iterator for LASIterator<'_, T> {
190    type Item = LASVectorFeature;
191
192    fn next(&mut self) -> Option<Self::Item> {
193        let las = &self.reader;
194        if self.index >= self.len {
195            None
196        } else if let Some(point) = las.get_feature(self.index) {
197            self.index += 1;
198            Some(point)
199        } else {
200            None
201        }
202    }
203}
204/// A feature reader trait with a callback-based approach
205impl<T: Reader> FeatureReader<(), Properties, LASPoint> for LASReader<T> {
206    type FeatureIterator<'a>
207        = LASIterator<'a, T>
208    where
209        T: 'a;
210
211    fn iter(&self) -> Self::FeatureIterator<'_> {
212        LASIterator { reader: self, index: 0, len: self.len() }
213    }
214
215    #[cfg(feature = "std")]
216    fn par_iter(&self, pool_size: usize, thread_id: usize) -> Self::FeatureIterator<'_> {
217        let pool_size = pool_size as u64;
218        let thread_id = thread_id as u64;
219        let start = self.len() * thread_id / pool_size;
220        let end = self.len() * (thread_id + 1) / pool_size;
221        LASIterator { reader: self, index: start, len: end }
222    }
223}
224
225/// The Public Header Block is followed by one or more Variable Length Records (There is one
226/// mandatory Variable Length Record, GeoKeyDirectoryTag). The number of Variable Length
227/// Records is specified in the "Number of Variable Length Records" field in the Public Header Block.
228/// The Variable Length Records must be accessed sequentially since the size of each variable length
229/// record is contained in the Variable Length Record Header. Each Variable Length Record Header
230/// is 54 bytes in length.
231///
232/// Each record is as follows:
233/// - Reserved `unsigned short 2 bytes`
234/// - User ID `char[16] 16 bytes`
235/// - Record ID `unsigned short 2 bytes`
236/// - Record Length After Header `unsigned short 2 bytes`
237/// - Description `char[32] 32 bytes`
238/// - optional data: variable size
239pub fn las_parse_variable_length_records<T: Reader>(
240    header: &LASHeader,
241    reader: &T,
242) -> BTreeMap<u32, LASExtendedVariableLengthRecord> {
243    let LASHeader { header_size, num_variable_length_records, .. } = header;
244    let mut res = BTreeMap::new();
245    let mut offset = *header_size as u64;
246    let mut i = 0;
247    while i < *num_variable_length_records {
248        let record = LASExtendedVariableLengthRecord::from_reader(reader, offset);
249        offset += 54 + record.record_length;
250        res.insert(record.record_id as u32, record);
251        i += 1;
252    }
253
254    res
255}
256
257/// WKT Parsing
258///
259/// For definition of WKT, we refer to Open Geospatial Consortium (OGC) specification “OpenGIS
260/// coordinate transformation service implementation specification” revision 1.00 released 12
261/// January 2001, section 7 (coordinate transformation services spec). This specification may be
262/// found at www.opengeospatial.org/standards/ct. As there are a few dialects of WKT, please note
263/// that LAS is not using the “ESRI WKT” dialect, which does not include TOWGS84 and authority
264/// nodes.
265/// - OGC MATH TRANSFORM WKT RECORD (2111)
266/// - OGC COORDINATE SYSTEM WKT (2112)
267///
268/// NOTE: It is required to use WKT if the point type is 6-10
269///
270/// ## Returns
271/// The WKT string if it exists
272pub fn build_wkt(
273    header: &LASHeader,
274    variable_length_records: &BTreeMap<u32, LASExtendedVariableLengthRecord>,
275    transformer: &mut Transformer,
276) -> Option<String> {
277    // 4th bit of global encoding must be set
278    if (header.encoding & (1 << 3)) != 0 {
279        return None;
280    }
281    // OGC MATH TRANSFORM WKT RECORD:
282    let wkt_math_ogc = variable_length_records.get(&2111).and_then(|v| v.data.clone());
283    // OGC COORDINATE SYSTEM WKT:
284    let wkt_coord_system_data = variable_length_records.get(&2112).and_then(|v| v.data.clone());
285    if wkt_math_ogc.is_none() && wkt_coord_system_data.is_none() {
286        return None;
287    }
288    let wkt_coord_system = wkt_math_ogc.or(wkt_coord_system_data).unwrap();
289    let wkt_str: String = String::from_utf8_lossy(&wkt_coord_system).into();
290    transformer.set_source(wkt_str.clone());
291
292    Some(wkt_str)
293}
294
295/// userID of "LASF_Projection" will contain at least 3 records:
296/// - GeoKeyDirectoryTag (34735)
297/// - GeoDoubleParamsTag (34736)
298/// - GeoASCIIParamsTag (34737)
299///
300/// Only the `GeoKeyDirectoryTag` record is required. This parses the `GeoKeyDirectoryTag`.
301/// This record contains the key values that define the coordinate system. A complete description
302/// can be found in the GeoTIFF format specification. Here is a summary from a programmatic point
303/// of view for someone interested in implementation.
304///
305/// The `GeoKeyDirectoryTag` is defined as just an array of unsigned short values. But,
306/// programmatically, the data can be seen as something like this:
307///
308/// ## Returns
309/// The parsed GeoKeyDirectory
310pub fn build_geo_key_directory(
311    variable_length_records: &BTreeMap<u32, LASExtendedVariableLengthRecord>,
312    transformer: &mut Transformer,
313    wkt_is_none: bool,
314) -> GeoStore {
315    let mut file_dir = GeoStore::default();
316    // GeoKeyDirectoryTag (34735)
317    let geokey_record = variable_length_records
318        .get(&(FieldTagNames::GeoKeyDirectory as u32))
319        .and_then(|v| v.data.clone());
320    if geokey_record.is_none() {
321        return GeoStore::default();
322    }
323    let raw_geo_keys: Vec<u16> = geokey_record
324        .unwrap()
325        .chunks_exact(2)
326        .map(|chunk| u16::from_le_bytes([chunk[0], chunk[1]]))
327        .collect();
328    // GeoDoubleParamsTag (34736)
329    let double_record = variable_length_records
330        .get(&(FieldTagNames::GeoDoubleParams as u32))
331        .and_then(|v| v.data.clone());
332    if let Some(double_record) = double_record {
333        file_dir.set(
334            FieldTagNames::GeoDoubleParams as u16,
335            double_record.to_vec(),
336            GeoTIFFTypes::DOUBLE,
337        );
338    }
339    // GeoAsciiParamsTag (34737)
340    let ascii_record = variable_length_records
341        .get(&(FieldTagNames::GeoAsciiParams as u32))
342        .and_then(|v| v.data.clone());
343    if let Some(ascii_record) = ascii_record {
344        file_dir.set(
345            FieldTagNames::GeoAsciiParams as u16,
346            ascii_record.to_vec(),
347            GeoTIFFTypes::ASCII,
348        );
349    }
350    let gkd = parse_geotiff_raw_geokeys(&raw_geo_keys, &file_dir);
351    if wkt_is_none {
352        build_transform_from_geo_keys(transformer, &gkd);
353    }
354
355    gkd
356}
gistools/readers/las/reader.rs

gistools/readers/las/
reader.rs
