Sebastian Flachmeier (recipient of three LASmoons)
UniGIS Master of Science, University of Salzburg, AUSTRIA
Bavarian Forest National Park, administration, Grafenau, GERMANY
Background:
The Bavarian Forest National Park is located in South-Eastern Germany, along the border with the Czech Republic. It has a total area of 240 km² and its elevation ranges from 600 to 1453 m. In 2002 a project called “High-Tech-Offensive Bayern” was started and a few first/last return LiDAR transects were flown to compute some forest metrics. The results showed that LiDAR has an advantage over other methods, because the laser was able to get readings from below the canopy. New full waveform scanner were developed that produced many more returns in the lower canopy. The National Park experimented with this technology in several projects and improved their algorithms for single tree detection. In 2012 the whole park was flown with full waveform and strategies for LiDAR based forest inventory for the whole National Park were developed. This is the data that is used in the following workflow description.
Goal:
Several versions of the LiDAR existed on the server of the administration that didn’t have the attributes we needed to reconstruct the original flight lines. The number of returns per pulse, the flight line IDs, and the GPS time stamps were missing. The goal was a workflow to create a LAStools workflow to convert the LiDAR from the original ASCII text files provided by the flight company into LAS or compressed LAZ files with all fields properly populated.
Data:
+ ALS data flown in 2012 by Milan Geoservice GmbH 650 m above ground with overlap.
+ full waveform sensor (RIEGL 560 / Q680i S) with up to 7 returns per shot
+ total of 11.080.835.164 returns
+ in 1102 ASCII files with *.asc extension (changed to *.txt to avoid confusion with ASC raster)
+ covered area of 1.25 kilometers
+ last return density of 17.37 returns per square meter
This data is provided by the administration of Bavarian Forest National Park. The workflow was part of a Master’s thesis to get the academic degree UniGIS Master of Science at the University of Salzburg.
LAStools processing:
The LiDAR was provided as 1102 ASCII text files named ’spur000001.txt‘ to ’spur001102.txt‘ that looked like this:
more spur000001.txt 4589319.747 5436773.357 685.837 49 106 1 215248.851500 4589320.051 5436773.751 683.155 46 24 2 215248.851500 4589320.101 5436773.772 686.183 66 87 1 215248.851503 […]
Positions 1 to 3 store the x, y, and z coordinate in meter [m]. Position 4 stores the „echo width“ in 0.1 nanoseconds [ns], position 5 stores the intensity, position 6 stores the return number, position 7 stores the GPS time stamp in seconds [s] of the current GPS week. The „number of returns (of given pulse)“ information is not explicitly stored and will need to be reconstructed in order, for example, to identify which returns are last returns. The conversion from ASCII text to LAZ was done with the txt2las command line shown below that incorporates these rationals:
- Although the ASCII files list the three coordinates with millimeter resolution (three decimal digits), we store only centimeter resolution which is sufficient to capture all the precision in a typical airborne LiDAR survey.
- After computing histograms of the „return number“ and the „echo width“ for all points with lasinfo and determining their maximal ranges it was decided to use point type 1 which can store up to 7 returns per shot and store the „echo width“ as an additional attribute of type 3 („unsigned short“) using „extra bytes“.
- The conversion from GPS time stamp in GPS week time to Adjusted Standard time was done by finding out the exact week during which Milan Geoservice GmbH carried out the survey and looking up the corresponding GPS week 1698 using this online GPS time calculator.
- Information about the Coordinate Reference System „DHDN / 3-degree Gauss-Kruger zone 4“ as reported in the meta data is added in form of EPSG code 31468 to each LAS file.
txt2las -i ascii\spur*.txt ^
-parse xyz0irt ^
-set_scale 0.01 0.01 0.01 ^
-week_to_adjusted 1698 ^
-add_attribute 3 "echo width" "of returning waveform [ns]" 0.1 0 0.1 ^
-epsg 31468 ^
-odir spur_raw -olaz ^
-cores 4
The 1102 ASCII files are now 1102 LAZ files. Because we switched from GPS week time to Adjusted Standard GPS time stamps we also need to set the „global encoding“ flag in the LAS header from 0 to 1 (see ASPRS LAS specification). We can do this in-place (i.e. without creating another set of files) using the following lasinfo command:
lasinfo -i spur_raw\spur*.laz ^
-nh -nv -nc ^
-set_global_encoding 1
To reconstruct the missing flight line information we look for gaps in the sequence of GPS time stamps by computing GPS time histograms with lasinfo and bins of 10 seconds in size:
lasinfo -i spur_raw\spur*.laz -merged ^
-histo gps_time 10 ^
-o spur_raw_all.txt
The resulting histogram exhibits the expected gaps in the GPS time stamps that happen when the survey plane leaves the target area and turns around to approach the next flight line. The subsequent histogram entries marked in red show gaps of 120 and 90 seconds respectively.
more spur_raw_all.txt [...] bin [27165909.595196404,27165919.595196255) has 3878890 bin [27165919.595196255,27165929.595196106) has 4314401 bin [27165929.595196106,27165939.595195957) has 435788 bin [27166049.595194317,27166059.595194168) has 1317998 bin [27166059.595194168,27166069.595194019) has 4432534 bin [27166069.595194019,27166079.59519387) has 4261732 [...] bin [27166239.595191486,27166249.595191337) has 3289819 bin [27166249.595191337,27166259.595191188) has 3865892 bin [27166259.595191188,27166269.595191039) has 1989794 bin [27166349.595189847,27166359.595189698) has 2539936 bin [27166359.595189698,27166369.595189549) has 3948358 bin [27166369.595189549,27166379.5951894) has 3955071 [...]
Now that we validated their existence, we use these gaps in the GPS time stamps to split the LiDAR back into the original flightlines it was collected in. Using lassplit we produce one file per flightline as follows:
lassplit -i spur_raw\spur*.laz -merged ^
-recover_flightlines_interval 10 ^
-odir strips_raw -o strip.laz
In the next step we repair the missing „number of returns (per pulse)“ field that was not provided in the ASCII file. This can be done with lasreturn assuming that the point records in each file are sorted by increasing GPS time stamp. This happens to be true already in our case as the original ASCII files where storing the LiDAR returns in acquisition order and we have not changed this order. If the point records are not yet in this order it can be created with lassort as follows. As these strips can have many points per file it may be necessary to run the new 64 bit executables by adding ‚-cpu64‘ to the command line in order to avoid running out of memory.
lassort -i strips_raw\strips*.laz ^
-gpstime -return_number ^
-odir strips_sorted -olaz ^
-cores 4 -cpu64
An order sorted by GPS time stamp is necessary as lasreturn expects point records with the same GPS time stamp (i.e. returns generated by the same laser pulse) to be back to back in the input file. To ‚-repair_number_of_returns‘ the tool will load all returns with the same GPS time stamp and update the „number of returns (per pulse)“ attribute of each return to the highest „return number“ of the loaded set.
lasreturn -i strips_sorted\strips*.laz ^
-repair_number_of_returns ^
-odir strips_repaired -olaz ^
-cores 4
In a final step we use las2las with the ‚-files_are_flightlines‘ option (or short ‚-faf‘) to set the „file source ID“ field in the LAS header and the „point source ID“ attribute of every point record in the file to the same unique value per strip. The first file in the folder will have all its field set to 1, the next file will have all its field set to 2, the next file to 3 and so on. Please do not run this on multiple cores for the time being.
las2las -i strips_repaired\strips*.laz ^
-files_are_flightlines ^
-odir strips_final -olaz
It’s always useful to run a final validation of the files using lasvalidate to reassure yourself and the people you will be sharing the data with that nothing funky has happened during any of these conversion steps.
lasvalidate -i strips_final\strip*.laz ^
-o strips_final\report.xml
And it can also be useful to add an overview in SHP or KML format to the delivery that can be created with lasboundary as follows:
lasboundary -i strips_final\strip*.laz ^
-overview -labels ^
-o strips_final\overview.kml
The result was 89 LAZ files (each containing one complete flightline) totaling 54 GB compared to 1102 ASCII files (each containing a slice of a flightline) totaling 574 GB.