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    fn par_iter(&self, pool_size: usize, thread_id: usize) -> Self::FeatureIterator<'_> {
216        let pool_size = pool_size as u64;
217        let thread_id = thread_id as u64;
218        let start = self.len() * thread_id / pool_size;
219        let end = self.len() * (thread_id + 1) / pool_size;
220        LASIterator { reader: self, index: start, len: end }
221    }
222}
223
224/// The Public Header Block is followed by one or more Variable Length Records (There is one
225/// mandatory Variable Length Record, GeoKeyDirectoryTag). The number of Variable Length
226/// Records is specified in the "Number of Variable Length Records" field in the Public Header Block.
227/// The Variable Length Records must be accessed sequentially since the size of each variable length
228/// record is contained in the Variable Length Record Header. Each Variable Length Record Header
229/// is 54 bytes in length.
230///
231/// Each record is as follows:
232/// - Reserved `unsigned short 2 bytes`
233/// - User ID `char[16] 16 bytes`
234/// - Record ID `unsigned short 2 bytes`
235/// - Record Length After Header `unsigned short 2 bytes`
236/// - Description `char[32] 32 bytes`
237/// - optional data: variable size
238pub fn las_parse_variable_length_records<T: Reader>(
239    header: &LASHeader,
240    reader: &T,
241) -> BTreeMap<u32, LASExtendedVariableLengthRecord> {
242    let LASHeader { header_size, num_variable_length_records, .. } = header;
243    let mut res = BTreeMap::new();
244    let mut offset = *header_size as u64;
245    let mut i = 0;
246    while i < *num_variable_length_records {
247        let record = LASExtendedVariableLengthRecord::from_reader(reader, offset);
248        offset += 54 + record.record_length;
249        res.insert(record.record_id as u32, record);
250        i += 1;
251    }
252
253    res
254}
255
256/// WKT Parsing
257///
258/// For definition of WKT, we refer to Open Geospatial Consortium (OGC) specification “OpenGIS
259/// coordinate transformation service implementation specification” revision 1.00 released 12
260/// January 2001, section 7 (coordinate transformation services spec). This specification may be
261/// found at www.opengeospatial.org/standards/ct. As there are a few dialects of WKT, please note
262/// that LAS is not using the “ESRI WKT” dialect, which does not include TOWGS84 and authority
263/// nodes.
264/// - OGC MATH TRANSFORM WKT RECORD (2111)
265/// - OGC COORDINATE SYSTEM WKT (2112)
266///
267/// NOTE: It is required to use WKT if the point type is 6-10
268///
269/// ## Returns
270/// The WKT string if it exists
271pub fn build_wkt(
272    header: &LASHeader,
273    variable_length_records: &BTreeMap<u32, LASExtendedVariableLengthRecord>,
274    transformer: &mut Transformer,
275) -> Option<String> {
276    // 4th bit of global encoding must be set
277    if (header.encoding & (1 << 3)) != 0 {
278        return None;
279    }
280    // OGC MATH TRANSFORM WKT RECORD:
281    let wkt_math_ogc = variable_length_records.get(&2111).and_then(|v| v.data.clone());
282    // OGC COORDINATE SYSTEM WKT:
283    let wkt_coord_system_data = variable_length_records.get(&2112).and_then(|v| v.data.clone());
284    if wkt_math_ogc.is_none() && wkt_coord_system_data.is_none() {
285        return None;
286    }
287    let wkt_coord_system = wkt_math_ogc.or(wkt_coord_system_data).unwrap();
288    let wkt_str: String = String::from_utf8_lossy(&wkt_coord_system).into();
289    transformer.set_source(wkt_str.clone());
290
291    Some(wkt_str)
292}
293
294/// userID of "LASF_Projection" will contain at least 3 records:
295/// - GeoKeyDirectoryTag (34735)
296/// - GeoDoubleParamsTag (34736)
297/// - GeoASCIIParamsTag (34737)
298///
299/// Only the `GeoKeyDirectoryTag` record is required. This parses the `GeoKeyDirectoryTag`.
300/// This record contains the key values that define the coordinate system. A complete description
301/// can be found in the GeoTIFF format specification. Here is a summary from a programmatic point
302/// of view for someone interested in implementation.
303///
304/// The `GeoKeyDirectoryTag` is defined as just an array of unsigned short values. But,
305/// programmatically, the data can be seen as something like this:
306///
307/// ## Returns
308/// The parsed GeoKeyDirectory
309pub fn build_geo_key_directory(
310    variable_length_records: &BTreeMap<u32, LASExtendedVariableLengthRecord>,
311    transformer: &mut Transformer,
312    wkt_is_none: bool,
313) -> GeoStore {
314    let mut file_dir = GeoStore::default();
315    // GeoKeyDirectoryTag (34735)
316    let geokey_record = variable_length_records
317        .get(&(FieldTagNames::GeoKeyDirectory as u32))
318        .and_then(|v| v.data.clone());
319    if geokey_record.is_none() {
320        return GeoStore::default();
321    }
322    let raw_geo_keys: Vec<u16> = geokey_record
323        .unwrap()
324        .chunks_exact(2)
325        .map(|chunk| u16::from_le_bytes([chunk[0], chunk[1]]))
326        .collect();
327    // GeoDoubleParamsTag (34736)
328    let double_record = variable_length_records
329        .get(&(FieldTagNames::GeoDoubleParams as u32))
330        .and_then(|v| v.data.clone());
331    if let Some(double_record) = double_record {
332        file_dir.set(
333            FieldTagNames::GeoDoubleParams as u16,
334            double_record.to_vec(),
335            GeoTIFFTypes::DOUBLE,
336        );
337    }
338    // GeoAsciiParamsTag (34737)
339    let ascii_record = variable_length_records
340        .get(&(FieldTagNames::GeoAsciiParams as u32))
341        .and_then(|v| v.data.clone());
342    if let Some(ascii_record) = ascii_record {
343        file_dir.set(
344            FieldTagNames::GeoAsciiParams as u16,
345            ascii_record.to_vec(),
346            GeoTIFFTypes::ASCII,
347        );
348    }
349    let gkd = parse_geotiff_raw_geokeys(&raw_geo_keys, &file_dir);
350    if wkt_is_none {
351        build_transform_from_geo_keys(transformer, &gkd);
352    }
353
354    gkd
355}
gistools/readers/las/reader.rs

gistools/readers/las/
reader.rs
