LASmoons: Leonidas Alagialoglou

Leonidas Alagialoglou (recipient of three LASmoons)
Multimedia Understanding Group, Aristotle University of Thessaloniki
Thessaloniki, GREECE

Canopy height is a fundamental geometric tree parameter in supporting sustainable forest management. Apart from the standard height measurement method using LiDAR instruments, other airborne measurement techniques, such as very high-resolution passive airborne imaging, have also shown to provide accurate estimations. However, both methods suffer from high cost and cannot be regularly repeated.

Preliminary results of predicted CHE based on multi-temporal satellite images against ground-truth LiDAR measurements. The 3rd column depicts pixel-wise absolute error of prediction. Last column depicts pixel-wise uncertainty estimation of the prediction (in means of 3 standard deviations).

In our study, we attempt to substitute airborne measurements with widely available satellite imagery. In addition to spatial and spectral correlations of a single-shot image, we seek to exploit temporal correlations of sequential lower resolution imagery. For this we use a convolutional variant of a recurrent neural network based model for estimating canopy height, based on a temporal sequence of Sentinel-2 images. Our model’s performance using sequential space borne imagery is shown to outperform the compared state-of-the-art methods based on costly airborne single-shot images as well as satellite images.

Digital Terrain Model of a part of the study area

The experimental study area of approximately 940 squared km is includes two national parks, Bavarian Forest National Park and Šumava National Park, which are located at the border between Germany and Czech Republic. LiDAR measurements of the area from 2017 and 2019 will be used as ground truth height measurements that have been provided by the national park’s authorities. Temporal sequences of Sentinel-2 imagery will be acquired from the Copernicus hub for canopy height estimation.

LAStools processing:
Accurate conversion of LAS files into DEM and DSM in order to acquire ground truth canopy height model.
1) Remove noise [lasthin, lasnoise]
2) Classify points into ground and non-ground [lasground, lasground_new]
3) Create DTMs and DSMs [lasthin, las2dem]

LASmoons: Zak Kus

Zak Kus (recipient of three LASmoons)
Topology Enthusiast
San Francisco, USA

While LiDAR data enables a lot of research and innovation in a lot of fields, it can also be used to create unique and visceral art. Using the high resolution data available, a 3D printer, and a long tool chain, we can create a physical, 3D topological map of the San Francisco bay area that shows off both the city’s hilly geology, and its unique skyline.


Test print of San Francisco’s Golden Gate Park.


Test print of San Francisco’s Golden Gate Park.

The ultimate goal of this project is to create an accurate, unique physical map of San Francisco, and the surrounding areas, which will be given to a loved one as a birthday gift. Using the data from the 2010 ARRA-CA GoldenGate survey, we can filter and process the raw lidar data into a DEM format using LAStools, which can be converted using a python script into a “water tight” 3D printable STL file.

While the data works fairly well out of the box, it does require a lot of manual editing, to remove noise spikes, and to delineate the coast line from the water in low lng areas. Interestingly, while many sophisticated tools exist to edit STLs that could in theory be used to clean up and prepare the files at the STL stage, few are capable of even opening files with so much detailed data. Using LAStools to manually classify, and remove unwanted data is the only way to achieve the desired level of detail in the final piece.

LiDAR data provided through USGS OpenTopography, using the ARRA-CA GoldenGate 2010 survey
+ Average point density of 3.33 pts/m^2 (though denser around SF)
+ Covers 2638 km^2 in total (only a ~100 km^2 subset is used)

LAStools processing:
Remove noise [lasnoise]
2) Manually clean up shorelines and problematic structures [lasview, laslayers]
3) Combine multiple tiles (to fit 3d printer) [lasmerge]
4) Create DEMs (asc format) for external tool to process [las2dem]

LASmoons: Martin Romain

Martin Romain (recipient of three LASmoons)
Marshall Islands Conservation Society
Majuro, Republic of the MARSHALL ISLANDS

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

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.

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]

Potential LAStools pipelines:
Removing Excessive Low Noise from Dense-Matching Point Clouds
Digital Pothole Removal: Clean Road Surface from Noisy Pix4D Point Cloud
Creating DTMs from dense-matched points of UAV imagery from SenseFly’s eBee

Philippines use Taal Vulcano Eruption as Opportunity to become Very First Asian Country with Open LiDAR

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 ^

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).

lasgrid -i Taal\DTM\*.laz ^
        -step 1 ^
        -use_tile_size 1000 ^
        -odir Taal\DTM_unbuffered ^

Without buffers the total file size in RasterLAZ format shrinks to 22 MB. Now we have the data in a format that can either be treated as a raster or as a point cloud. Hence we can use laspublish and quickly create a visualization of the Taal Vulcano Island with Potree which we then copied onto our university Web space for you to play with.  This was he are able to instantly create an 3D visualization portal that lets anyone do various simple and also more complex measurements.

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.

We would like to acknowledge the UP Training Center for Applied Geodesy and Photogrammetry (UP TCAGP) and their DREAM and PHIL LiDAR program for providing easy and unencumbered open access to this data with a license that encourages data reuse and repurposing. Kudos for being first in Asia to make open LiDAR happen!!!

Converting Rasters from inefficient ASCII XYZ to more compact LAZ or TIF Formats

The German state of Brandenburg has recently started to provide many of their basic geospatial data as open data, such as digital ortophotos in TIF and JPG formats, vertical and horizontal control points in gzipped XML format, LOD1 and LOD2 building models in zipped GML format, topographic maps from 1:10000 to 1:100000 in zipped TIF and PDF formats, cadastral data in zipped XML and TIF formats, as well as LiDAR-derived 1m DTM rasters and image-derived 1m DSM rasters both in zipped XYZ ASCII format. All this data is provided with the user-friendly license called “Datenlizenz Deutschland Namensnennung 2.0“. In this article we show how to convert the 1m DTM rasters and the 1m DSM rasters  from verbose XYZ ASCII to more compact LAZ or TIF rasters.


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 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 “” with 2001 columns by 2001 rows. Here is how the 4004001 lines looks:

250000.0 5886000.0 15.284
250001.0 5886000.0 15.277
250002.0 5886000.0 15.273
250003.0 5886000.0 15.275
250004.0 5886000.0 15.289
250005.0 5886000.0 15.314
251994.0 5888000.0 13.565
251995.0 5888000.0 13.567
251996.0 5888000.0 13.565
251997.0 5888000.0 13.565
251998.0 5888000.0 13.564
251999.0 5888000.0 13.564
252000.0 5888000.0 13.565

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:

txt2las -i dgm\*.xyz ^
        -set_scale 1.0 1.0 0.001 ^
        -epsg 25833 ^
        -odir temp -olaz ^
        -cores 4

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:

txt2las -i temp\*.laz ^
        -odir ascii -otxt ^
        -cores 4

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.

lasgrid -i temp\*.laz ^
        -subsquare 0.25 ^
        -step 1 -average ^
        -use_tile_size 2000 ^
        -odir dgm -olaz ^
        -cores 4

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)

LASmoons: Volga Lipwoni

Volga Lipwoni (recipient of three LASmoons)
Department of Geography, School of Earth and Environment
University of Canterbury, NEW ZEALAND

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.

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.

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]

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:
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

Another German State Goes Open LiDAR: Saxony

Finally some really good news out of Saxony. 😊 After North Rhine-Westphalia and Thuringia released the first significant amounts of open geospatial data in Germany in a one-two punch in January 2017, we now have a third German state opening their entire tax-payer-funded geospatial data holdings to the tax-paying public via a simple and very easy-to-use online download portal. Welcome to the open data party, Saxony!!!

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.

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:

more Greifensteine\
352000 5613999 636.26
352001 5613999 636.27
352002 5613999 636.28
352003 5613999 636.27
352004 5613999 636.24

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.

txt2las -i greifensteine\333* ^
        -set_scale 1 1 0.01 ^
        -epsg 25833 ^

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.

384,000,000 bytes

2,684,820 333525610_dgm1.laz
2,590,516 333525612_dgm1.laz
2,853,851 333545610_dgm1.laz
2,795,430 333545612_dgm1.laz
10,924,617 bytes

Using blast2dem we then create a hill-shaded version of the 1 meter DTM in order to overlay a visual representation of the DTM onto Google Earth.

blast2dem -i greifensteine\333*_dgm1.laz ^
          -merged ^
          -step 1 ^
          -hillshade ^
          -o greifensteine.png

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.

gdalwarp -tr 1 1 ^
         -r cubic ^
         greifensteine\dop20c_33352_5612.tif ^

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.

lascolor -i greifensteine\333525612_dgm1.laz ^
         -image greifensteine\dop1m_33352_5612.tif ^
         -odix _rgb -olaz

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.

lasview -i greifensteine\333525612_dgm1_rgb.laz

We thank the “Staatsbetrieb Geobasisinformation und Vermessung Sachsen (GeoSN)” for giving us easy access to the 1 meter DTM and the 20 cm orthophoto that we have used in this article through their new open geodata portal as open data under the user-friendly license “Datenlizenz Deutschland – Namensnennung – Version 2.0.

National Open LiDAR Strategy of Latvia humiliates Germany, Austria, and other European “Closed Data” States

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 kilometer tiles 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.

lasinfo -i 00_tiles_raw\*.las ^
        -compute_density ^
        -histo z 5 ^
        -histo intensity 256 ^
        -histo user_data 1 ^
        -histo scan_angle 1 ^
        -histo point_source 1 ^
        -histo gps_time 10 ^
        -odir 01_quality -odix _info -otxt ^
        -cores 3
lasvalidate -i 00_tiles_raw\*.las ^
            -no_CRS_fail ^
            -o 01_quality\report.xml

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.

laszip -i 00_tiles_raw\*.las ^
       -epsg 3059 ^
       -cores 2

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.

lasgrid -i 00_tiles_raw\*.laz ^
        -keep_last ^
        -step 2 ^
        -density_16bit ^
        -false -set_min_max 2 8 ^
        -odir 01_quality -odix _d_2_8 -opng ^
        -cores 3

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.

lasoverlap -i 00_tiles_raw\*.laz ^
           -step 2 ^
           -min_diff 0.1 -max_diff 0.2 ^
           -odir 01_quality -opng ^
           -cores 3

The results of the quality checks with lasgrid and lasoverlap are shown below.

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.

lasindex -i 00_tiles_raw\*.laz ^
         -cores 3

las2dem -i 00_tiles_raw\*.laz ^
        -keep_class 2 9 ^
        -buffered 25 ^
        -step 1 ^
        -use_orig_bb ^
        -odir Latvia\02_dtm_1m -olaz ^
        -cores 3

blast2dem -i 02_dtm_1m\*.laz ^
          -merged ^
          -hillshade ^
          -step 1 ^
          -o dtm_1m.png

las2dem -i 00_tiles_raw\*.laz ^
        -drop_class 1 7 ^
        -buffered 10 ^
        -spike_free 1.5 ^
        -step 1 ^
        -use_orig_bb ^
        -odir 03_dsm_1m -olaz ^
        -cores 3

blast2dem -i 03_dsm_1m\*.laz ^
          -merged ^
          -hillshade ^
          -step 1 ^
          -o dsm_1m.png

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.

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!

Clean DTM from Agisoft Photogrammetric Points of Urban Scene

We occasionally get permission to distribute a nice data sets and blog about how to best process it with LAStools because this gets around having to pay our “outrageous” consulting fees. (-: This time we received a nice photogrammetric point cloud of the Tafawa Balewa Square in Lagos Island, Lagos, Nigeria. This area is part of the central business district of Lagos and characterized by high-rise buildings. The Tafawa Balewa Square was constructed in 1972 over the site of a defunct track for horse racing and is bounded by Awolowo road, Cable street, Force road, Catholic Mission street and the 26-story Independence House. We want to create a nice Digital Terrain Model from the dense-matching point cloud that was generated with Photoscan by AgiSoft and – as always with photogrammetry – we have to take special care of low noise points. The final result is shown below. All processing commands used are here.

After downloading the data it is useful to familiarize yourself with the file, which can be done with lasview, lasinfo, and lasgrid using the command lines shown below. According to the lasinfo report there are around 47 million points points with RGB colors in the file and their average density is around 100 points per square meter.

lasview -i 0_raw\TafawaBalewa.laz

lasinfo -i 0_raw\TafawaBalewa.laz ^
        -cd -histo intensity 256 ^
        -histo z 1 ^
        -odir 1_quality -odix _info -otxt

lasgrid -i 0_raw\TafawaBalewa.laz ^
        -step 1 ^
        -density ^
        -false -set_min_max 50 150 ^
        -odir 1_quality -odix _d_50_150 -opng

The average point density value of 100 from the lasinfo report suggests that 50 as the minimum and 150 as the maximum are good false color ramp values for a map showing how the point density per square meter is distributed.

Color-coded point density: blue equals 50 or less and red means 150 or more points per square meter.

We use lastile to create a buffered tiling for the 47 million points. We use a tile size of 200 meters and request a large buffer of 50 meters around every tile because there are large buildings in the survey areas. We also flag buffer points as withheld so they can be easily be dropped later.

lastile -i 0_raw\TafawaBalewa.laz ^
        -tile_size 200 -buffer 50 -flag_as_withheld ^
        -odir 2_tiles_raw -o tafawa.laz

If you inspect the resulting tiles – such as ‘tafawa_544000_712600.laz’ as shown here – with lasview you will see the kind of low noise that is shown below and that may cause a ground classification algorithm. While our lasground software is able to deal with a certain amount of low noise – if there are too many it will likely latch onto them. Therefore we will first generate a subset of points that has as few as possible of such low noise points.

Typical low noise in dense-matching photogrammetry points in urban scene.

Next we use a sequence of three LAStools modules, namely lasthinlasground, and lasheight to classify this photogrammtric point cloud into ground and non-ground points. All processing commands used are here. First we use lasthin to give the point the classification code 8 that is closest to the 50th percentile in elevation within every 50 centimeter by 50 centimeter cell (but only if the cells containing at least 20 points).

lasthin -i 2_tiles_raw\tafawa*.laz ^
        -step 0.5 ^
        -percentile 50 20 ^
        -classify_as 8 ^
        -odir 3_tiles_median_50cm -olaz ^
        -cores 3

Next we use lasground to ground-classify only the points that have classification code 8 (i.e. by ignoring those with classification codes 0) and set their classification code to ground (2) or non-ground (1). Because of the large buildings in this urban scene we use ‘-metro’ which uses a large step size of 50 meters for the pre-processing. This also sets the internally used bulge parameter to 5.0 which you can see if you run the tool in verbose ‘-v’ mode. In three different trial runs we determined that forcing the bulge parameter to be 0.5 instead gave better results. The bulge and the spike parameters can be useful to vary in order to improve ground classification results (also see the README file).

lasground -i 3_tiles_median_50cm\tafawa*.laz ^
          -ignore_class 0 ^
          -metro -bulge 0.5 ^
          -odir 4_tiles_ground_50cm -olaz ^
          -cores 3

The resulting ground points are a subset with a resolution of 50 centimeter that is good enough to create a DTM with meter resolution, which we do with las2dem command line shown below. We really like storing DTM elevation rasters to the LAZ point format because it is a more compressed way of storing elevation rasters compared to ASC, BIL, TIF, or IMG. It also makes the raster output a natural input to subsequent LAStools processing steps.

las2dem -i 4_tiles_ground_50cm\tafawa*.laz ^
        -keep_class 2 ^
        -step 1 -kill 100 ^
        -use_tile_bb ^
        -odir 5_tiles_dtm_1m -olaz ^
        -cores 3

Finally we use blast2dem to create a seamless hill-shaded version of our 1 meter DTM from on-the-fly merged elevation rasters. This is the DTM pictured at the beginning of this article.

blast2dem -i 5_tiles_dtm_1m\tafawa*.laz -merged ^
          -step 1 ^
          -hillshade ^
          -o dtm_1m.png

The corresponding DSM pictured at the beginning of this article was generated with the two command lines below by first keeping only the 95th percentile highest elevation of every 50 cm by 50 cm cell with lasthin (which remove spurious high noise points) and then by triangulating the surviving points with blast2dem into a seamless TIN that is also hill-shaded and rasterized with 1 meter resolution. Running the 64 bit version of lasthin (note the ‘-cpu64‘ switch) allows us to work on the entire data set (rather than its tiles version) at once, where the standard 32 bit version may run out of memory.

lasthin -i 0_raw\TafawaBalewa.laz ^
        -cpu64 ^
        -step 0.5 ^
        -percentile 95 20 ^
        -o 0_raw\TafawaBalewa_p95_50.laz

blast2dem -i 0_raw\TafawaBalewa_p95_50.laz ^
          -step 1 ^
          -hillshade ^
          -o dsm_1m.png

In order to generate the final DTM at higher resolution we use lasheight to pull all points into the ground class that lie within a 5 cm distance vertically below or a 15 cm distance vertically above the triangulated surface of ground points computed in the previous step. You could experiment with other values here to be less or more conservative about pulling detail into the ground class.

lasheight -i 4_tiles_ground_50cm\tafawa*.laz ^
          -classify_between -0.05 0.15 2 ^
          -odir 6_tiles_ground -olaz ^
          -cores 3

We repeat the same processing step as before las2dem to create the raster DTM tiles, but this time with a resolution of 25 cm.

las2dem -i 6_tiles_ground\tafawa*.laz ^
        -keep_class 2 ^
        -step 0.25 -kill 100 ^
        -use_tile_bb ^
        -odir 7_tiles_dtm_25cm -olaz ^
        -cores 3

And we again use blast2dem to create a seamless hill-shaded version of the DTM from on-the-fly merged elevation rasters, but this time with a resolution of 25 cm. This is the DTM shown below. All processing commands used are here.

blast2dem -i 7_tiles_dtm_25cm\tafawa*.laz -merged ^
          -step 0.25 ^
          -hillshade ^
          -o dtm_25cm.png

Hill-shade of final DTM with resolution of 25 cm.

Using Open LiDAR to Remove Low Noise from Photogrammetric UAV Point Clouds

We collected drone imagery of the restored “Kance” tavern during the lunch stop of the UAV Tartu summer school field trip (actually the organizer Marko Kohv did that, as I was busy introducing students to SUP boarding). With Agisoft PhotoScan we then processed the images into point clouds below the deck of the historical “Jõmmu” barge on the way home (actually Marko did that, because I was busy enjoying the view of the wetlands in the afternoon sun). The resulting data set with 7,855,699 points is shown below and can be downloaded here.

7,855,699 points produced with Agisoft Photoscan

Generating points using photogrammetric techniques in scenes containing water bodies tends to be problematic as dense blotches of noise points above and below the water surface are common as you can see in the picture below. Especially the low points are troublesome as they adversely affect ground classification which results in poor Digital Elevation Models (DTMs).

Clusters of low noise points nearly 2 meters below the actual surface in water areas.

In a previous article we have described a LAStools workflow that can remove excessive low noise. In this article here we use external information about the topography of the area to clean our photogrammetry points. How convenient that the Estonian Land Board has just released their entire LiDAR archives as open data.

Following these instructions here you can download the available open LiDAR for this area, which has the map sheet index 475681. Alternatively you can download the four currently available data sets here flown in spring 2010, in summer 2013, in spring 2014, and in summer 2017. In the following we will use the one flown in spring 2014.

We can view both data sets simultaneously in lasview. By adding ‘-faf’ to the command-line we can switch back and forth between the two data sets by pressing ‘0’ and ‘1’.

lasview -i Kantsi.laz ^
        -i 475681_2014_tava.laz ^
        -points 10000000 ^

We find cut the 1 km by 1 km LiDAR tile down to a 250 m by 250 m tile that nicely surrounds our photogrammetric point set using the following las2las command-line:

las2las -i 475681_2014_tava.laz ^
        -inside_tile 681485 6475375 250 ^
        -o LiDAR_Kantsi.laz

lasview -i Kantsi.laz ^
        -i LiDAR_Kantsi.laz ^
        -points 10000000 ^

Scrutinizing the two data sets we quickly find that there is a miss-alignment between the dense imagery-derived and the comparatively sparse LiDAR point clouds. With lasview we investigate the difference between the two point clouds by hovering over a point from one point cloud and pressing <i> and then hovering over a somewhat corresponding point from the other point cloud and pressing <SHIFT>+<i>. We measure displacements of around 2 meters vertically and of around 3 to 3.5 meter in total.

Before we can use the LiDAR points to remove the low noise from the photogrammetric points we must align them properly. For simple translation errors this can be done with a new feature that was recently added to lasview. Make sure to download the latest version (190404 or newer) of LAStools to follow the steps shown in the image sequence below.

las2las -i Kantsi.laz ^
        -translate_xyz 0.89 -1.90 2.51 ^
        -o Kantsi_shifted.laz

lasview -i Kantsi_shifted.laz ^
        -i LiDAR_Kantsi.laz ^
        -points 10000000 ^

The result looks good in the sense that both sides of the photogrammetric roof are reasonably well aligned with the LiDAR. But there is still a shift along the roof so we repeat the same thing once more as shown in the next image sequence:

We use a suitable displacement vector and apply it to the photogrammetry points, shifting them again:

las2las -i Kantsi_shifted.laz ^
        -translate_xyz -1.98 -0.95 0.01 ^
        -o Kantsi_shifted_again.laz

lasview -i Kantsi_shifted_again.laz ^
        -i LiDAR_Kantsi.laz ^
        -points 10000000 ^

The result is still not perfect as there is also some rotational error and you may find another software such as Cloud Compare more suited to align the two point clouds, but for this exercise the alignment shall suffice. Below you see the match between the photogrammetry points and the LiDAR TIN before and after shifting the photogrammetry points with the two (interactively determined) displacement vectors.

The final steps of this exercise use las2dem and the already ground-classified LiDAR compute a 1 meter DTM, which we then use as input to lasheight. We classify the photogrammetry points using their height above this set of ground points with 1 meter spacing: points that are 40 centimeter or more below the LiDAR DTM are classified as noise (7), points that are between 40 below to 1 meter above the LiDAR DTM are classified to a temporary class (here we choose 8) that has those points that could potentially be ground points. This will help, for example, with subsequent ground classification as large parts of the photogrammetry points – namely those on top of buildings and in higher vegetation – can be ignored from the start by a ground classification algorithm such as lasground.

las2dem -i LiDAR_Kantsi.laz ^
        -keep_class 2 ^
        -kill 1000 ^
        -o LiDAR_Kantsi_dtm_1m.bil
lasheight -i Kantsi_shifted_again.laz ^
          -ground_points LiDAR_Kantsi_dtm_1m.bil ^
          -classify_below -0.4 7 ^
          -classify_between -0.4 1.0 8 ^
          -o Kantsi_cleaned.laz

Below the results we have achieved after “roughly” aligning the two point clouds with some new lasview tricks and then using the LiDAR elevations to classify the photogrammetry points into “low noise”, “potential ground”, and “all else”.

We thank the Estonian Land Board for providing open data with a permissive license. Special thanks also go to the organizers of the UAV Summer School in Tartu, Estonia and the European Regional Development Fund for funding this event. Especially fun was the fabulous excursion to the Emajõe-Suursoo Nature Reserve and through to Lake Peipus aboard, overboard and aboveboard the historical barge “Jõmmu”. If you look carefully you can also find the barge in the photogrammetry point cloud. The photogrammetry data used here was acquired during our lunch stop.

Fun aboard and overboard the historical barge “Jõmmu”.