Background: Hydrological models require various input data for flood vulnerability mapping. An important input data for flood vulnerability mapping is the DTM over which flow is being routed. DTMs are generated using cartography, ground surveying, digital aerial photogrammetry, interferometric SAR (InSAR), LiDAR amongst other means. The accuracy of high resolution DTMs minimize errors that may emanate from input data when conducting hydrological modelling, especially in small built-up catchment areas. This research involves the application of digital aerial photogrammetry to generate point clouds which can subsequently be utilized for flood vulnerability mapping.
photogrammetry point cloud
photogrammetic DSM
Goal:
To consolidate on previous gains in using LAStools to generate DTMs required for flood vulnerability mapping. The suitability of these DTMs will be subsequently validated for flood vulnerability analysis. These results will be compared with other DTMs in order to determine the uncertainty associated with the use of such DTMs for flood vulnerability mapping.
photogrammetry point cloud
photogrammetric DSM
Data:
+ high-resolution photogrammetry point cloud and DSM for Lagos Island, Ikorodu and Ajah Nigeria – – – imagery obtained with an Ebee Sensefly drone flight – – – photogrammetry point cloud generated with Photoscan by AgiSoft + rainfall data + classified LiDAR point cloud with a resolution of 1 pulse per square meter obtained for the study area from the Lagos State Government
Photogrammetry point cloud
Photogrammetric DSM
LAStools processing: 1) tile large photogrammetry point cloud into tiles with buffer [lastile] 2) mark set of points whose z coordinate is a certain percentile of that of their neighbors [lasthin] 3) remove isolated low points from the set of marked points [lasnoise] 4) classify marked points into ground and non-ground [lasground] 5) pull in points close above and below the ground [lasheight] 6) create Digital Terrain Model (DTM) from ground points [las2dem] 7) merge and hillshade individual raster DTMs [blast2dem]
In the last article about the Livox MID-40 LiDAR scan of Samara we quality checked the data, aligned the flight lines and cleaned the remaining spurious scan lines. In this article we will process this data into the standard products. A focus will be on generating a smooth “median ground” surface from the “fluffy” scanner data. You can get the flight lines here and follow along with the processing after downloading LAStools. Below is one result of this work.
The first processing step will be to tile the strips into tiles that contain fewer points for faster and also parallel processing. One quick “flat terrain” trick first. Often there are spurious points that are far above or below the terrain. For a relatively flat area these can be easily be identified by computing a histogram elevation values with lasinfo and then eliminated with simple drop-filters on the z coordinate.
The relevant excerpts of the output of the lasinfo report are shown below:
[…] z coordinate histogram with bin size 1.000000 bin -104 has 1 bin 5 has 1 bin 11 has 273762 bin 12 has 1387999 bin 13 has 5598767 bin 14 has 36100225 bin 15 has 53521371 […] bin 59 has 60308 bin 60 has 26313 bin 61 has 284 bin 65 has 10 bin 66 has 31 bin 67 has 12 bin 68 has 1 bin 83 has 3 bin 84 has 4 bin 93 has 31 bin 94 has 93 bin 95 has 17 […]
The few points below 11 meters and above 61 meters in elevation can be considered spurious. In the initial tiling step with lastile we add simple elevation filters to drop those points on-the-fly from the buffered tiles. The importance of buffers when processing LiDAR in tiles is discussed in this article. With lastile we create tiles of size 125 meters with a buffer of 20 meters, while removing the points identified as spurious with the appropriate filters. Because the input strips have their “file source ID” in the LAS header correctly set, we use ‘-apply_file_source_ID’ to set the “point source ID” of every point to this value. This preserves the information of which point comes from which flight line.
This produces 49 buffered tiles that will now be processed similarly to the workflow outlined for another lower-priced system that generates similarly “fluffy” point clouds on hard surfaces, the Velodyne HLD-32E, described here and here. What do we mean with “fluffy”? We cut out a 1 meter slice across the road with the new ‘-keep_profile’ filter and las2las and inspect it with lasview.
In the view below we pressed hot key twice ‘]’ to exaggerate the z scale. The “fuzziness” is that thickness of the point cloud in the middle of this flat tar road. It is around 20 to 25 centimeters and is equally evident in both flight lines. What is the correct ground surface through this 20 to 25 centimeter “thick” road? We will compute a “mean ground” that roughly falls into the middle of this “fluffy” surface,
slice of 1 meter width across the tar road in front of Omael store
The next three lasthin runs mark a sub set of low candidate points for our lasground filtering. In every 25 cm by 25 cm, every 33 cm by 33 cm and every 50 cm by 50 cm area we reclassify the point closest to the 10th percentile as class 8. In the first call to lasthin we put all other points into class 1.
Below you can see the resulting points of the 10th percentile classified as class 8 in red.
Operating only on the points classified as 8 (i.e. ignoring those classified as 1) we then run a ground classification with lasground using the following command line, which creates a “low ground” classification. .
Since this is an open road this classifies most of the red points as ground points.
Using lasheight we then create a “thick ground” by pulling all those points into the ground surface that are between 5 centimeter below and 17 centimeter above the “low ground”. For visualization purposes we temporarily use class 6 to capture this thickened ground.
We go back to lasthin and reclassify in every 50 cm by 50 cm area the point closest to the 50th percentile as class 8. This is what we call the “median ground”.
The final “median ground” points are shown in red below. These are the points we will use to eventually compute the DTM.
We complete the fully automated classification available in LAStools by running lasclassify with the following options. See the README file for what these options mean. Note that we move the “thick ground” from the temporary class 6 to the proper class 2. The “median ground” continues to be in class 8.
Before the resulting tiles are published or shared with others we should remove the temporary buffers, which is done with lastile – the same tool that created the buffers.
Below a screenshot of the resulting Potree 3D Web portal rendered with Potree Desktop. Inspecting the classification will reveal a number of errors that could be tweaked manually with lasview. How the point colors were generated is not described here but I used Google satellite imagery and mapped it with lascolor to the points. The elevation colors are mapped from 14 meters to 25 meters. The intensity image may help us understand why the black tar road on the left hand side that runs from the “Las Palmeras Condos” to the beach in “Cangrejal” has no samples. It seems the intensity is lower on this side which indicates that the drone may have flown higher here – too high to for the road to reflect enough photons. The yellow view of return type indicates that despite it’s multi-return capability, the Livox MID-40 LiDAR is mostly collecting single returns.
colored by classification
colored by Google Earth RGB
colored by elevation
corored by intensity
colored by return type
The penetration capability of the Livox MID-40 LiDAR was less good than we had hoped for. Below thick vegetation we have too few points on the ground to give us a good digital terrain model. In the visualization below you can see that below the dense vegetation there are large black areas which are completely void of points.
all returns
only building, thick and median ground points
Now we produce the standard product DTM and DSM at a resolution of 50 cm. Because the total area is not that big we generate temporary tiles in “raster LAZ” with las2dem and merge them into a single GeoTiff with blast2dem.
After weeks of planning the helicopter finally came. Since we have been cooperating with Stereocarto‘s San Jose office for quite some time, we were finally able to lure their helicopter down to Samara. Taking off from Liberia it flew about an hour back and forth across Samara and the adjacent beaches. Processing this data will take weeks (if not months) to complete and analyzing it even longer. One of the first applications will be the usual extraction of the bare earth terrain topography followed by the detection of large changes in comparison to our 2012 data. Such changes might be the result of permitted construction or earthworks, but are often illegal remodeling or removal of mountains to “create” building lots or filling in wetland areas with dirt and debris to “gain” land. Finding and quantifying such environmental damages and reporting about the process will be the first primary objective of this project.
flight planning on the beach
actual flight flown
rapidlasso’s quaint office
full house in Samara with LiDAR processing on multiple machines
Before the flight arrived I raked the word “LAStools” into the sand …… it came out nicely in the images taken with a Hasselblad camera.
Background: As a low-lying coastal nation, the Republic of the Marshall Islands (RMI) is at the forefront of exposure to climate change impacts. RMI has a strong dependence on natural resources and biodiversity not only for food and income but also for culture and livelihood. However, these resources are threatened by rising sea levels and associated coastal hazards (king tides, storm surges, wave run-up, saltwater intrusion, erosion). This project aims at addressing the lack of technical capacity and available data to implement effective risk reduction and adaptation measures, with a particular focus on inundation mapping and local evacuation planning in population centers.
Typical low-lying coastal area of the Republic of the Marshall
Goal:
This project intends to use LAStools to generate a DEM of the inhabited sections of 3 remote atolls (Aur, Ebon, Likiep) and 1 island (Mejit). The resulting DEM will be used to produce an inundation exposure model (and map) under variable sea level rise projections for each site. The ultimate goal is to integrate the results into each site’s disaster risk reduction strategy (long-term outcome) and present it through community consultations in schools, community centers, and council houses.
Data:
+ Aerial imagery of 11.5 square kilometers of land (6.3% of total national landmass) using DJI Matrice 200 V2 & DJI Zenmuse X5S with a minimum overlap of 75/75 and maximum altitude of 120m.
LAStools processing: 1) tile large point cloud into tiles with buffer [lastile] 2) remove noise points [lasthin, lasnoise] 3) classify points into ground and non-ground [lasground] 4) create Digital Terrain Models and Digital Surface Models [lasthin, las2dem]
UPDATE: As of January 30th also orthophotos and classified LAZ tiles are available for download.
It took just a few years of nagging, a vulcanic eruption, and then a few more weeks of nagging but now it has happened. The Philippines have become the first country in Asia to offering LiDAR as open data for free and unencumbered download. The portal created by the UP Training Center for Applied Geodesy and Photogrammetry (UP TCAGP) and their DREAM and PHIL LiDAR program already offers LiDAR-derived 1 meter DTM and DSM data flown between 2013 and 2017 as part of a national mission to aquire flood mapping data for a certain area around the Taal Vulcano. In the coming days orthophotos and the classified LiDAR point cloud will be added (at the moment the data is still undergoing another quality assurance review process).
As a quick test we went to the new online portal and downloaded the 34 DTM raster tiles that cover the Taal Vulcano Lake as seen in the screenshot below.
Downloading the area-of-interest is easy with LiPAD’s nice download portal.
The downloaded 1 meter DTM tiles are in TIF format and each cover an area of 1000 by 1000 meter. However, they are overlapping because they have a 50 meter buffer, so that each raster contains elevation samples organized in 1100 columns by 1100 rows plus “no data” values. We use two LAStools commands to remove the buffers. First we use our new demzip to turn the TIF to RasterLAZ format. Use demzip from version 200131 of LAStools (or newer) as older releases did not handle “no data” values correctly.
demzip -i Taal\DTM\*.tif ^
-olaz
The conversion from TIF to RasterLAZ also reduces the total file size for the 34 files from 157 MB to 27 MB. Next we remove the buffers using a new functionality in lasgrid (make sure you have the latest LAStools version 200112 or newer).
laspublish -i Taal\DTM_unbuffered ^
-elevation ^
-odir Taal\DTM_portal ^
-o TaalVulcanoIsland.html ^
-title "DTM of Taal Vulcano Island" ^
-description "DTM of Taal Vulcano Island" ^
-olaz -overwrite
Below we see the result visualized with the Desktop version of Potree. You can access the interactive portal we have created here with any Web browser.
Visualizing the 1 meter DTM of Taal Vulcano Island as RasterLAZ point cloud with Potree to instantly create interactive portal allowing simple measurements that give an intuition about the height and the size of the vulcanic formation that makes up Taal Vulcano Island.
Four 2000 by 2000 meter tiles of the Brandenburg 1m DTM.
One particularity about most official German and Austrian rasters (anywhere else?) is that they sample the elevations in the corners rather than in the center of each raster cell. Here a one square kilometer raster tile of 1 meter resolution will have 1001 columns by 1001 rows instead of the more familiar 1000 by 1000 layout. While this corner-based representation does have some benefits, we convert these rasters in to the more common area-based representation using new functionality recently added to lasgrid.
After downloading one sample DTM tile such as dgm_33250-5886.zip we find three files in the zip folder. Two files with meta data and license information and the actual data file, which is a 2 km by 2km corner-based raster tile called “dgm_33250-5886.xyz” with 2001 columns by 2001 rows. Here is how the 4004001 lines looks:
The first step is to convert these XYZ rasters to LAZ format. We do this with txt2las as shown below. In case the vertical datum is the “Deutsches Haupthoehennetz 2016” we should also add ‘-vertical_dhhn2016’ but not sure at the moment:
For 84 files this reduces the size by a factor of 31 or compresses it down to 3.2 percent of the original, namely from 8.45 GB for raw XYZ to 277 MB for LAZ. So far we have really just converted a list of x, y and z coordinates from verbose ASCII to more compact LAZ. We can easily go back to ASCII with las2txt whenever needed:
Next we use lasgrid to convert from a corner-based raster to an area-based raster using the new option ‘-subsquare 0.2’ which replaces each input point by four points that are displaced by all possibilities of adding +/- 0.2 in x and y. We then average the exactly four points that fall into each relevant raster cell with option ‘-average’ and clip the output to the meaningful 2000 columns by 2000 rows with ‘-use_tile_size 2000’. You need to get the most recent version of LAStools to have these options.
Instead of RasterLAZ you can also choose the TIF, BIL, IMG, or ASC format here. The final result are standard 1 meter elevation products with 2000 columns by 2000 rows with the averaged elevation sample being associated with the center of the raster cell. The lasinforeport for a sample tile is shown at the end of this article.
You may proceed to optimize the RasterLAZ for area-of-interest queries by reordering the raster into a space-filling curve with lassort or lasoptimize and compute a spatial index. You may also classify the RasterLAZ elevation samples, for example, into building, high, medium, and low vegetation, ground, and other common classifications with lasclip or lascolor. You may also add RGB or intensity values to the RasterLAZ elevation samples using the orthophotos that are also available as open data with lascolor. These are some of the benefits of RasterLAZ beyond efficient storage and access.
We like to acknowledge the LGB (Landesvermessung und Geobasisinformation Brandenburg) for providing state-wide coverage of their geospatial data holdings as easily downloadable open data with the user-friendly Deutschland Namensnennung 2.0 license. But we also would like to ask to please add the raw LiDAR point clouds to the open data portal. The storage savings in going from ASCII XYZ to LAZ for the DTM and DSM rasters should free enough space to host the LiDAR … (-;
lasinfo (200112) report for 'dgm_33\DGM_33250-5886.laz'
reporting all LAS header entries:
file signature: 'LASF'
file source ID: 0
global_encoding: 0
project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
version major.minor: 1.2
system identifier: 'raster compressed as LAZ points'
generating software: 'LAStools (c) by rapidlasso GmbH'
file creation day/year: 13/20
header size: 227
offset to point data: 455
number var. length records: 2
point data format: 0
point data record length: 20
number of point records: 4000000
number of points by return: 4000000 0 0 0 0
scale factor x y z: 0.5 0.5 0.001
offset x y z: 200000 5800000 0
min x y z: 250000.5 5886000.5 13.419
max x y z: 251999.5 5887999.5 33.848
variable length header record 1 of 2:
reserved 0
user ID 'Raster LAZ'
record ID 7113
length after header 80
description 'by LAStools of rapidlasso GmbH'
ncols 2000
nrows 2000
llx 250000
lly 5886000
stepx 1
stepy 1
sigmaxy <not set>
variable length header record 2 of 2:
reserved 0
user ID 'LASF_Projection'
record ID 34735
length after header 40
description 'by LAStools of rapidlasso GmbH'
GeoKeyDirectoryTag version 1.1.0 number of keys 4
key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
key 3072 tiff_tag_location 0 count 1 value_offset 25833 - ProjectedCSTypeGeoKey: ETRS89 / UTM 33N
key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter
LASzip compression (version 3.4r3 c2 50000): POINT10 2
reporting minimum and maximum for all LAS point record entries ...
X 100001 103999
Y 172001 175999
Z 13419 33848
intensity 0 0
return_number 1 1
number_of_returns 1 1
edge_of_flight_line 0 0
scan_direction_flag 0 0
classification 0 0
scan_angle_rank 0 0
user_data 0 0
point_source_ID 0 0
number of first returns: 4000000
number of intermediate returns: 0
number of last returns: 4000000
number of single returns: 4000000
overview over number of returns of given pulse: 4000000 0 0 0 0 0 0
histogram of classification of points:
4000000 never classified (0)
Background: Structure from motion (SfM) photogrammetry, has emerged as an effective tool to accurately extract three-dimensional (3D) structures from a series of overlapping two-dimensional (2D) Unmanned aerial vehicles (UAVs) images. The bid to switch from the current labour-intensive, and time consuming forestry inventory practices has seen a lot of interest geared towards understanding the use of SfM photogrammetry to derive forest metrics (Iglhaut et al., 2019). There are a range of commercial, free and open source SfM photogrammetric software packages that can be used to process UAV images into 3D point clouds. Selection of the most appropriate package has become an important issue for most projects (Turner, Lucieer, & Wallace, 2013). A comparison of software performance in terms of accuracy, processing times and related costs would help foresters in deciding the best tool for the job.
Typical point cloud derived with SfM software from UAV imagery.
Goal:
The study will generate 3D point clouds of images of a young forest trial and LAStools will be used to derive canopy height models (CHM) for computing tree heights. Tree heights from LiDAR data will serve as a baseline for accuracy assessment of heights derived from the point clouds.
Data:
+ 422 UAV images processed into 3D point clouds using ten (10) different commercial and open source SfM software packages
LAStools processing: 1) tile large point cloud into tiles with buffer [lastile] 2) remove noise points [lasthin, lasnoise] 3) classify points into ground and non-ground [lasground] 4) create Digital Terrain Modelsand Digital Surface Models [lasthin, las2dem] 5) produce Canopy Height Models for computing tree heights [lasheight, las2dem]
References:
Iglhaut, J., Cabo, C., Puliti, S., Piermattei, L., O’Connor, J., & Rosette, J. (2019). Structure from motion photogrammetry in forestry: A review. Current Forestry Reports, 5(3), 155-168. doi:https://doi.org/10.1007/s40725-019-00094-3
Turner, D., Lucieer, A., & Wallace, L. (2013). Direct georeferencing of ultrahigh-resolution UAV imagery. EEE Transactions on Geoscience and Remote Sensing, 52(5), 2738-2745. doi:10.1109/TGRS.2013.2265295
A while back we had a first look at the Single Photon LiDAR from Leica’s SPL100 sensor (that eventually turned out just to be an SPL99 because one beamlet or one receiver in the 10 by 10 array was broken and did not produce any returns). Today we are taking a closer look at a strategy to remove the excessive noise in the raw Single Photon LiDAR data from a “proper” SPL100 sensor (where all of the 100 beamlets are firing) that was flown in 2017 in Navarra, Spain.
Profile through original points on top of generated DTM.
The data was provided as open data by the cartography section of Navarra’s Government and is available via a simple download FTP portal. We describe the LAStools processing steps that were used to eliminate the excessive noise and to generate a smooth DTM. In the following we are using the originally released version of the data, that we obtained shortly after the portal went online that seems to be a bit more “raw” than the current files available now. One starndard quality check with lasinfo was done with:
Upon inspecting the lasinfo report we suggest a few changes in how to store this Single Photon LiDAR data for more efficient hosting via an online portal. We perform these changes here before starting the actual processing. First we use the las2las call shown below to fix an error in the global encoding bits, remove an irrelevant VLR, re-scale the coordinates from millimeter to centimeters, re-offset the coordinates to nice numbers, and – what is by far the most crucial change for better compression – remap the beamlet ID stored in the ‘user data’ field as described in an earlier article.
Then we use two lassort calls, one to maximize compression and one to improve spatial coherence. One lassort call rearranges the points in increasing order first based on the GPS time stamps, then breaks ties based on the user data field (that stores the beamlet ID), and finally stores the returns of every beamlet ordered by return number. We also add spatial reference information in this step. The other lassort call rearranges the points into a spatially coherent layout. It uses a Z-order sort with the granularity of 50 meter by 50 meter buckets of points. Within each bucket the point order from the prior sort is kept.
Now we start the usual processing workflow by tiling the data with lastile into smaller 500 meter by 500 meter tiles with a 25 meter buffer. We also set the pre-existing point classification in the data to zero as we will compute our own later.
We notice that a large amount of the noise has intensity values below 1000. We are still a bit puzzled where those intensity values come from and what exactly they mean in a Single Photon LiDAR system. But it works. We run las2las with a “filtered transform” to set classification of all points whose intensity value is 1000 or less to the classification code 7 (aka “noise”).
We then ignore this “easy-to-identify” noise and go after the remaining one with lasnoise by ignoring classification code 7 and setting the newly identified noise to classification code 9 – not because it’s “water” (the usual meaning of class 9) but because these points are drawn with a distinct blue color when checking the result with lasview.
Of the surviving non-noise points we then use lasthin to reclassify the point closest to the 20th elevation percentile per 50 cm by 50 cm area with classification code 8 (for all areas that have more than 5 non-noise points per 50 cm by 50 cm area. We repeat the same for every 1 meter by 1 meter area.
We then perform a more agressive second noise removal step one with lasnoise using only those points with classification code 8, namely those non-noise points that were the 20th elevation percentile in either a 50 cm by 50 cm cell or a 1 meter by 1 meter cell. This can be done by ignoring classification code 0, 7, and 9. We mark those noise points as 6 so they appear orange in the point cloud with lasview.
The 20th elevation percentile points that survive the last noise removal are then classified into ground (2) and non-ground (1) points with lasground_new by ignoring all other points, namely those with classification codes 0, 6, 7, and 9.
These images below illustrate the steps we took. They also show that not all data was used and might give you ideas where to tweak our workflow for even better results.
First returns are red, intermediate returns are green, last returns are blue and single returns are yellow.
Colored by intensity. Most noise has low intensity values.
Unused points. Another workflow could try to use more of this data.
Final ground points used for generating the smooth DTM.
Finally we raster the ground points into 1 meter Digital Terrain Model (DTM) rasters with las2dem and store the result (without buffers) to the RasterLAZ format.
The hillshaded DTM that is result of the entire sequence of processing steps described above is shown below.
DTM from ground classification created with LAStools
For comparison we generate the same DTM using the originally provided classification. According to the README file the original ground points are classified with code 22 in areas of flight line overlap and as the usual code 2 elsewhere. Hence we must use both classification codes to construct the DTM. We do this analogue to the earlier processing steps with the three LAStools commands lastile, las2dem, and blast2dem below.
Below the hillshaded DTM generated from the ground classification that was provided with the LiDAR when it was originally released as open data.
DTM from ground classification of originally released data.
In the meantime Andorra’s SPL data have been updated with a newer version in the open data portal. The new version of the data contains a much better ground classification that might have been improved manually as the new files now have the the string ‘cam’ instead of ‘ca’ in the file name, which probably means ‘classified automatically and manually’ instead of the original ‘classified automatically’. We decided not to switch to the new data release as it seemed less “raw” than the original release. For example there are suddenly points with GPS times and returns counts and numbers of zero in the file that seem synthetic. But we also computed the hillshaded DTM for the new release which is shown below.
DTM from ground classification of newly released data.
We thank the cartography section of Navarra’s Government for providing their LiDAR as open data. This not only allows re-purposing expensive data paid for by public taxes but also generates additional value, encourages citizen science, and provides educational opportunity and insights such as this blog article.
Currently available via the online portal are the LiDAR-derived raster Digital Terrain Model (DTM) at 1 meter resolution (DGM 1m) for everything flown since 2015 and and at 2 meter resolution (DGM 2m) or 20 meter resolution (DGM 20m) for the entire state. The horizontal coordinates use UTM zone 33 with ETRS89 (aka EPSG code 25833) and the vertical coordinate uses the “Deutsche Haupthöhennetz 2016” or “DHHN2016” (aka EPSG code 7837). Also available are orthophotos at 20 cm (!!!) resolution (DOP 20cm).
Overview of current LiDAR holdings. Areas flown 2015 or later have LAS files and 1 meter rasters. Others have LiDAR as ASCII files and lower resolution rasters.
Offline – by ordering through either this online form or that online form – you can also get the 5 meter DTM and the 10 meter DTM, the raw LiDAR point clouds, LiDAR intensity rasters, hill-shaded DTM rasters, as well as the 1 meter and the 2 meter Digital Surface Model (DSM) for a small administrative fee that ranges between 25 EUR and 500 EUR depending on the effort involved.
Our immediate thought is to get a copy on the entire raw LiDAR points clouds (available as LAS 1.2 files for all data acquired since 2015 and as ASCII text for earlier acquisitions) and find some portal willing to hosts this data online. We are already in contact with the land survey of Saxony to discuss this option and/or alternate plans.
Products available in Saxony’s new open data portal.
Downloading the 1m DTM raster.
Downloading the 20 cm ortophoto.
Let’s have a look at the data. First we download four 2 km by 2 km tiles of the 1 meter DTM raster for an area surrounding the so called “Greifensteine” using the interactive map of the download portal, which are provided as simple XYZ text. Here a look at the contents of one ot these tiles:
Note that the elevation are not sampled in the center of every 1 meter by 1 meter cell but exactly on the full meter coordinate pair, which seems especially common in German-speaking countries. Using txt2las we convert these XYZ rasters to LAZ format and add geo-referencing information for more efficient subsequent processing.
Below you see that going from XYZ to LAZ reduces the amount of data from 366 MB to 10.4 MB, meaning that the data on disk becomes over 35 times smaller. The ability of LASzip to compress elevation rasters was first noted during the search for missing airliner MH370 and resulted in our new LAZ-based compressor for height grid called DEMzip. The resulting LAZ files now also include geo-referencing information.
Below the result that nicely shows how the penetrating laser of the LiDAR allows us to strip away the forest to see interesting geological features in the bare-earth terrain.
In a second exercise we use the available RGB orthophoto images to color one of the DTM tiles and explore it using lasview. For this we download the image for the top left of the four tiles that covers the area containing the “Greifensteine” from the interactive download portal for orthophotos. As the resolution of the TIF image is 20 cm and that of the DTM is only 1 meter, we first down-sample the TIF using gdalwarp of GDAL.
If you are not yet using GDAL today is a good day to start. It nicely complements the point cloud processing functionality of LAStools for raster inputs. Next we use lascolor to give each elevation pixel of the DTM stored in LAZ format its corresponding color from the orthophoto.
Now we can view the colored DTM in LAZ format interactively with lasview or any other LiDAR viewing software and turn on the RGB colors from the orthophoto as needed to understand the scene.
Latvia, officially the Republic of Latvia, is a country in the Baltic region of Northern Europe has around 2 million inhabitants, a territory of 65 thousand square kilometers and – since recently – also a fabulous open LiDAR policy. Here is a list of 65939 tiles in LAS format available for free download that cover the entire country with airborne LiDAR with a density from 4 to 6 pulses per square meters. The data is classified into ground, building, vegetation, water, low noise, and a few other classifications. It is licensed Creative Commons CC0 1.0 – meaning that you can copy, modify, and distribute the data, even for commercial purposes, all without asking permission. And there is a simple and functional interactive download portal where you can easily download individual tiles.
Interactive open LiDAR download portal of Latvia.
We downloaded the 5 by 5 block of square kilometertiles matching “4311-32-XX.las” for checking the quality and creating a 1m DTM and a 1m DSM raster. You can follow along after downloading the latest version of LAStools.
Quality Checking
We first run lasvalidate and lasinfo on the downloaded LAS files and then immediately compress them with laszip because multi-core processing of uncompressed LAS files will quickly overwhelm our file system, make processing I/O bound, and result in overall longer processing times with CPUs waiting idly for data to be loaded from the drives.
Despite already excluding a missing Coordinate Reference System (CRS) from being a reason to fail (the lasinfo reports show that the downloaded LAS files do not have any geo-referencing information) lasvalidate still reports a few failing files, but scrutinizing the resulting XML file ‘report.xml’ shows only minor issues.
Usually during laszip compression we do not alter the contents of a file, but here we also add the EPSG code 3059 for CRS “LKS92 / Latvia TM” as we turn bulky LAS files into slim LAZ files so we don’t have to specify it in all future processing steps.
Compression reduces the total size of the 25 tiles from over 4.1 GB to below 0.6 GB.
Next we use lasgrid to visualize the last return density which corresponds to the pulse density of the LiDAR survey. We map each 2 by 2 meter pixel where the last return density is 2 or less to blue and each 2 by 2 meter pixel it is 8 or more to red.
This we follow by the mandatory lasoverlap check for flight line overlap and alignment where we map the number of overlapping swaths as well as the worst vertical difference between overlapping swaths to a color that allows for quick visual quality checking.
The results of the quality checks with lasgrid and lasoverlap are shown below.
last return density: blue = 2 per square meter, red = 8 per square meter
flight line overlap: blue – 1 flight line, aqua – 2 flight lines, yellow – 3 flight lines, orange – 4 flight lines, red – 5 or more flight lines
flight line alignment: white – error less than 10 cm, red or blue – error more than 20 cm
Raster Derivative Generation
Now we use first las2dem to create a Digital Terrain Model (DTM) and a Digital Surface Model (DSM) in RasterLAZ format and then use blast2dem to create merged and hill-shaded versions of both. Because we will use on-the-fly buffering to avoid edge effects along tile boundaries we first spatially index the data using lasindex for more efficient access to the points from neighboring tiles.
Because the overlaid imagery does not look as nice in our new Google Earth installation, below are the DTM and DSM at versions down-sampled to 25% of their original size.
hill-shaded DTM down-sampled to 25%
hill-shaded DSM down-sampled to 25%
Many thanks to SunGIS from Latvia who tweeted us about the Open LiDAR after we chatted about it during the Foss4G 2019 gala dinner. Kudos to the Latvian Geospatial Information Agency (LGIA) for implementing a modern national geospatial policy that created opportunity for maximal return of investment by opening the expensive tax-payer funded LiDAR data for re-purposing and innovation without barriers. Kudos!