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

Coco Loco for LiDAR: Open Data and 3D GIS Carnival Cruise into Sao Paulo

The municipality of São Paulo just dropped a MASSIVE amount of open LiDAR upon the geospatial community – a whopping 33 billion points. São Paulo city now has their own Entwine Point Tiles (EPT) service available for everyone interested in downloading, visualizing, or otherwise accessing this amazing data set. The EPT point cloud service can be accessed here and an interactive 3D portal that allows you to explore this massive data with nothing but a Web browser can be accessed here. It simple serves these EPT tiles via Potree.


What you see above is the beautiful Matarazzo Building (Edifício Matarazzo in Portuguese) where the City Hall of Sao Paulo is located. On May 27 2019 we had a one day workshop in the similarly beautiful Martinelli building (Edifício Martinelli in Portuguese) to brainstorm how to process the data and I seem to remember being rather VOCAL about repeating the success of Open LiDAR in Guadalajara and also providing easy access to the data for Sao Paolo. KUDOS to those who made it happen.

This data was originally collected for the Urban Agriculture project “São Paulo: Growing Farmers’ Income, Shrinking Urban Sprawl” that was supported after winning the 2016 MAYORS CHALLENGE by Bloomburg Philanthropies. At this point I would like to thank former New York mayor Michael Rubens Bloomburg  (“Hey Mike, do they also constantly miss-spell your family name? So annoying!”) for ending his 2020 bid for the Democratic presidential nomination.

This project stands on the shoulders of many many giants such as Howard Butler and his team who worked for years to create Entwine, Markus Schuetz who came out of nowhere and suddenly dropped Potree upon our community, and the many contributors in person at OSGeo hackfests and codesprints or remotely via patches and pull requests on the software infrastructure behind this effort. Don’t forget to support these folks!

Above other visualizations of the the Matarazzo Building where the City Hall of Sao Paulo is located. Below an overview of the city after zooming out a bit.

Surprise Release of Airborne LiDAR in Germany’s “Closed Data State” Bavaria

You have guessed correctly. This is mostly fake news as our “Freistaat” (read “Free State“) of Bavaria continues to tightly guard all of its tax-payer funded geospatial basis data for no good reason. Our other “Free State” – that of Thuringia – has become one poster child of open data in Germany with North Rhine-Westphalia being the original one. But there is “some” open LiDAR in Bavaria now.

The authors of a recent paper on change detection in urban areas have published two interesting airborne LiDAR data sets from 2008 and 2009 for the town of Abenberg in [Hebel, Arens, Stilla, 2013]. What is interesting about these data sets is that they (a) were flown with a forward looking laser scanner with flight trajectories from 4 different directions (as illustrated in the image below) and (b) were surveyed again the same way in the following year for temporal change detection.


Our lasview is able to visualize how the four different flight lines scan this house from four different directions, once we have reconstructed a properly populated LAZ file with flight line information, return numbers, and unique GPS time stamps.

The data is provided for download here as zipped ASCII files that have one line per returns containing 11 comma-separated values. Below is a sample of the first 10 lines of the 2008 data:

1, 290.243, 28.663, -11.787, 0.060, -0.052, 0.997, 517.3170, -58.6934, 313.0817, 52
1, 290.208, 28.203, -11.825, 0.062, -0.056, 0.996, 517.3167, -58.6934, 313.0817, 49
1, 290.182, 27.739, -11.852, 0.063, -0.055, 0.997, 517.3164, -58.6935, 313.0817, 53
1, 290.165, 27.272, -11.866, 0.061, -0.058, 0.996, 517.3161, -58.6935, 313.0817, 53
1, 290.163, 26.800, -11.858, 0.061, -0.053, 0.997, 517.3157, -58.6935, 313.0817, 68
1, 290.152, 26.334, -11.864, 0.059, -0.054, 0.997, 517.3154, -58.6936, 313.0817, 57
1, 290.092, 25.882, -11.938, 0.050, -0.057, 0.997, 517.3151, -58.6936, 313.0817, 57
1, 290.103, 25.406, -11.911, 0.043, -0.058, 0.997, 517.3147, -58.6937, 313.0817, 63
1, 290.067, 24.947, -11.952, 0.043, -0.061, 0.997, 517.3144, -58.6937, 313.0817, 63
1, 290.034, 24.488, -11.989, 0.044, -0.063, 0.997, 517.3141, -58.6937, 313.0817, 56

The first number is either a classification into ground, vegetation, or other surface, or represents an identifier for a planar shape that the return is part of. The next three numbers in red are the x, y, and z coordinate of the LiDAR point in a local coordinate system. The next three numbers in green are the x, y, and z coordinates of an estimated surface normal. The next three numbers in blue are the x, y, and z coordinates of the sensor position. The last number is the intensity of the LiDAR return.

This textual representation makes it difficult to efficiently load the data into most LiDAR processing software. Also several attributes such as return number, number of returns, flight line ID, and GPS time stamps are missing.

We use this as an opportunity for a little exercise in converting ASCII to LAZ while preserving any “additional attributes” using the “extra bytes” functionality available since the LAS 1.4 specification. This is a timely experiment as the LAS Working Group of the ASPRS is currently contemplating how to standardize some useful “additional attributes”. Here you can download the resulting files:

In order to replicate these steps, please get your hands of the very latest version of LAStools. First we use txt2las to convert the TXT file to LAZ format as follows:

txt2las -i abenberg_data_2008.txt ^
        -set_point_type 1 ^
        -parse 0xyz123456i ^
        -set_scale 0.001 0.001 0.001 ^
        -add_attribute 3 "planar shape ID" "preliminary classification" ^
        -add_attribute 4 "normal x coord" "local normal direction estimate" 0.001 ^
        -add_attribute 4 "normal y coord" "local normal direction estimate" 0.001 ^
        -add_attribute 4 "normal z coord" "local normal direction estimate" 0.001 ^
        -add_attribute 6 "sensor x coord" "sensor position" 0.0001 ^
        -add_attribute 6 "sensor y coord" "sensor position" 0.0001 ^
        -add_attribute 6 "sensor z coord" "sensor position" 0.0001 ^
        -odix _temp1 -olaz

We use a back and forth of lasview and las2las with option ‘-subseq 12345 67890’ to interactively find the exact index of the return where each flightline is ending and the next one is starting. The command below allows you to visualize the the trajectories.

lasview -i abenberg_data_2008_temp1.laz ^
        -copy_attribute_into_x 4 ^
        -copy_attribute_into_y 5 ^
        -copy_attribute_into_z 6 ^
        -points 10000000 ^
        -point_size 5


By hovering with the mouse over a point where the trajectory end and pressing ‘i’ like info we can query the coordinates and attributes of a point. The console output also lists the index of the point in the file. We use this index as the start index for the manual search of the exact index where the flight lines really ends by piping the coordinates as text to stdout and looking for a “jump” in coordinates that indicates the start of a new flightline as shown below.

las2las -i abenberg_data_2008_temp1.laz ^
        -subseq 1213485 1213505 ^
        -oparse xyz -otxt -stdout

-279.846 -98.442 -11.973
-279.984 -98.900 -12.123
-280.150 -99.357 -12.311
-280.195 -99.825 -12.332
-280.229 -100.294 -12.340
-280.275 -100.763 -12.363
-280.302 -101.233 -12.361
-280.344 -101.700 -12.379
-280.371 -102.171 -12.376
-280.156 -102.661 -12.042
110.259 304.077 3.873
110.646 304.118 3.925
111.036 304.154 3.970
111.424 304.167 3.982
111.811 304.211 4.037
112.201 304.252 4.088
112.588 304.278 4.118
112.976 304.315 4.164
113.364 304.336 4.186
113.752 304.344 4.192

Then we use the found index to seperate the first flight line form the rest with with two more runs of las2las:

las2las -i abenberg_data_2008_temp1.laz ^
         -subseq 0 1213495 ^
         -odix _strip1 -olaz

las2las -i abenberg_data_2008_temp1.laz ^
        -subseq 1213495 100000000 ^
        -odix _rest234 -olaz

We then repeat this procedure for the rest until we have four individual strips. Next we use las2las to change the point type from 0 to 1 to have a GPS time attribute that we will then populate with “fake” but useful GPS time stamps:

las2las -i abenberg_data_2008_temp1_strip*.laz ^
        -set_point_type 1 ^
        -odix _pt1 -olaz

We noticed in the original text file that subsequent groups of points often have the exact same value for the three numbers in blue that are the x, y, and z coordinates of the sensor position. This suggests that those points are multiple returns from the same laser shot. We wrote a small tool using LASlib that exploits this observed pattern to recover the “return number” and the “number of returns” attribute for each point and to set a “fake” but unique GPS time for each such group of returns. You can download the source code for this tool here. We run this tool for each strip with a different start GPS time:

lasrecover -i abenberg_data_2008_temp1_strip1_pt1.laz ^
           -gpstime_start 1000000 ^
           -odix _rec -olaz

lasrecover -i abenberg_data_2008_temp1_strip2_pt1.laz ^
           -gpstime_start 2000000 ^
           -odix _rec -olaz

lasrecover -i abenberg_data_2008_temp1_strip3_pt1.laz ^
           -gpstime_start 3000000 ^
           -odix _rec -olaz

lasrecover -i abenberg_data_2008_temp1_strip4_pt1.laz ^
           -gpstime_start 4000000 ^
           -odix _rec -olaz

Finally we merged the four strips back into one file using lasmerge:

lasmerge -i abenberg_data_2008_temp1_strip1_pt1_rec.laz ^
         -i abenberg_data_2008_temp1_strip2_pt1_rec.laz ^
         -i abenberg_data_2008_temp1_strip3_pt1_rec.laz ^
         -i abenberg_data_2008_temp1_strip4_pt1_rec.laz ^
         -faf ^
         -o abenberg_data_2008_temp2.laz

With lasview we are now able to visualize not only the multiple returns per shot but also the different angles from which the laser scanner has observed the scene.

Finally we use las2las and a “filtered transform” with the custom classfication code provided in the first “additional attribute” to populate the official LAS classification codes as ratified by the ASPRS. First we turn code 1 (“Ground level”) into ground points.

las2las -i abenberg_data_2008_temp2.laz ^
        -keep_attribute_between 0 1 1 ^
        -filtered_transform ^
        -set_classification 2 ^
        -o abenberg_data_2008_temp3.laz


Then we turn code 5 (“Vegetation”) into high vegetation points, code 6 (“Other Surface”) into keypoints, and code 9 or higher (“Planar shape #”) into building points.

las2las -i abenberg_data_2008_temp3.laz ^
        -keep_attribute_between 0 5 5 ^
        -filtered_transform ^
        -set_classification 5 ^
        -o abenberg_data_2008_temp4.laz

las2las -i abenberg_data_2008_temp4.laz ^
        -keep_attribute_between 0 6 6 ^
        -filtered_transform ^
        -set_classification 8 ^
        -o abenberg_data_2008_temp5.laz

las2las -i abenberg_data_2008_temp5.laz ^
         -keep_attribute_above 0 8 ^
         -filtered_transform ^
         -set_classification 6 ^
         -o abenberg_data_2008.laz

Then we are done. Here you can download the resulting files:


Hebel, M., Arens, M., Stilla, U., 2013. Change detection in urban areas by object-based analysis and on-the-fly comparison of multi-view ALS data. ISPRS Journal of Photogrammetry and Remote Sensing 86, pp. 52-64. [DOI: 10.1016/j.isprsjprs.2013.09.005] [PDF]

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

Removing Noise from Single Photon LiDAR to Generate a Smooth DTM

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:

lasinfo -i 0_raw\*.laz ^
        -cd ^
        -histo intensity 1 ^
        -histo user_data 1 ^
        -histo point_source 1 ^
        -histo gps_time 10 ^
        -odir 1_quality -odix _info -otxt

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.

las2las -i 0_raw\*.laz ^
        -rescale 0.01 0.01 0.01 ^
        -auto_reoffset ^
        -set_global_encoding_gps_bit 1 ^
        -remove_vlr 1 ^
        -map_user_data beamlet_ID_map.txt ^
        -odir 2_fix_rescale_reoffset_remap -olaz ^
        -cores 3

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.

lassort -i 2_fix_rescale_reoffset_remap\*.laz ^
        -epsg 25830 ^
        -gps_time ^
        -user_data ^
        -return_number ^
        -odir 2_maximum_compression -olaz ^
        -cores 3

lassort -i 2_maximum_compression\*.laz ^
        -bucket_size 50 ^
        -odir 2_spatial_coherence -olaz ^
        -cores 3

The resulting optimized nine tiles are around 200 MB each and can be downloaded as one file here or as individual tiles here:

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.

lastile -i 2_spatial_coherence\*.laz ^
        -set_classification 0 ^
        -tile_size 500 -buffer 25 -flag_as_withheld ^
        -odir 3_buffered -o yecora.laz

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

las2las -i 3_buffered\*.laz ^
        -keep_intensity_below 1000 ^
        -filtered_transform ^
        -set_classification 7 ^
        -odir 4_intensity_denoised -olaz ^
        -cores 3

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.

 lasnoise -i 4_intensity_denoised\*.laz ^
         -ignore_class 7 ^
         -step_xy 1.0 -step_z 0.2 ^
         -isolated 5 ^
         -classify_as 9 ^
         -odir 4_isolation_denoised -olaz ^
         -cores 3

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.

lasthin -i 4_isolation_denoised\*.laz ^
        -ignore_class 7 9 ^
        -step 0.5 -percentile 20 5 ^
        -classify_as 8 ^
        -odir 5_thinned_p20_050cm -olaz ^
        -cores 3

lasthin -i 5_thinned_p20_050cm\*.laz ^
        -ignore_class 7 9 ^
        -step 1.0 -percentile 20 5 ^
        -classify_as 8 ^
        -odir 5_thinned_p20_100cm -olaz ^
        -cores 3

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.

lasnoise -i 5_thinned_p20_100cm\*.laz ^
         -ignore_class 0 7 9 ^
         -step_xy 2.0 -step_z 0.2 ^
         -isolated 1 ^
         -classify_as 6 ^
         -odir 5_thinned_p20_100cm_denoised -olaz ^
         -cores 3

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.

lasground_new -i 5_thinned_p20_100cm_denoised\*.laz ^
              -ignore_class 0 6 7 9 ^
              -town ^
              -odir 5_tiles_ground_050cm -olaz ^
              -cores 3

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.

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.

las2dem -i 5_tiles_ground_050cm\*.laz ^
        -keep_class 2 ^
        -step 1.0 ^
        -use_tile_bb ^
        -odir 6_tiles_dtm_100cm -olaz ^
        -cores 3

Finally we merged all RasterLAZ tiles into one and compute the final hillshaded DTM with blast2dem.

blast2dem -i 6_tiles_dtm_100cm\*.laz -merged ^
          -step 1.0 ^
          -hillshade ^
          -o yecora_dtm_100cm.png

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.

lastile -i 2_spatial_coherence\*.laz ^
        -tile_size 500 -buffer 25 -flag_as_withheld ^
        -odir 3_tiles_buffered_orig -o yecora.laz

las2dem -i 3_tiles_buffered_orig\*.laz ^
        -keep_class 2 22 ^
        -step 1.0 ^
        -use_tile_bb ^
        -odir 6_tiles_dtm_100cm_orig -olaz ^
        -cores 3

blast2dem -i 6_tiles_dtm_100cm_orig\*.laz -merged ^
          -step 1.0 ^
          -hillshade ^
          -o yecora_dtm_100cm_orig.png

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.

Another European Country Opens LiDAR: Welcome to the Party, Slovakia!

We got a little note from Vítězslav Moudrý from CULS pointing out that the Geodesy, Cartography and Cadastre Authority of the Slovak Republic has started releasing LiDAR as open data on their interactive Web portal. Congratulations, Slovakia!!! Welcome to the Open Data Party!!! We managed to download some data starting from this Web portal link and describe the process of obtaining LiDAR data from the Low Tatras mountain range in central Slovakia with pictures below.


(1) click the new “data export” link


(2) change the export selection to “Shape”


(3) change the file format to “LAZ”


(4) zoom to a colored area-of-interest


(5) zoom further and draw a nice polygon


(6) edit polygon into nice shape and realize heart is red because area is too big


(7) zoom further and draw polygon smaller than 2 square kilometer


(8) when polygon turns green, accept license, enter email address and export


(9) short wait and you get download link to such an archive


(10) license conditions: PDF auto-translated from Slovak to English



(11) LiDAR are spatially indexed flight lines clipped to area-of-interest


(12) all return density: blue = 20 and red = 50 returns per square meter

lasgrid -i LowTatras\*.laz -merged ^
        -step 2 -point_density_16bit ^
        -false -set_min_max 20 50 ^
        -o LowTatras\density_all_returns_20_50.png

(13) last return density: blue = 4 and red = 40 last returns per square meter

lasgrid -i LowTatras\*.laz -merged ^
        -keep_last ^
        -step 2 -point_density_16bit ^
        -false -set_min_max 4 40 ^
        -o LowTatras\density_last_returns_4_40.png

(14) ground return density: blue = 4 and red = 40 ground returns per square meter

lasgrid -i LowTatras\*.laz -merged ^
        -keep_classification 2 ^
        -step 2 -point_density_16bit ^
        -false -set_min_max 4 40 ^
        -o LowTatras\density_ground_returns_4_40.png

(15) flight line difference image: white <= +/- 10 cm and red/blue >= +/- 20 cm

lasoverlap -i LowTatras\*.laz -faf ^
           -drop_classification 7 18 ^
           -min_diff 0.1 -max_diff 0.2 ^
           -o LowTatras\overlap_10cm_20cm.png

Finally we compute a DSM and a corresponding DTM using the already existing ground classification with BLAST using the command sequence shown below.


lasthin -i LowTatras\*.laz -merged ^
        -drop_classification 7 18 ^
        -step 0.5 -highest ^
        -o LowTatras\highest_50cm.laz

blast2dem -i LowTatras\highest_50cm.laz ^
          -hillshade ^
          -o LowTatras -o dsm_1m_hillshaded.png

blast2dem -i LowTatras\*.laz -merged ^
          -keep_classification 2 ^
          -thin_with_grid 0.5 ^
          -hillshade ^
          -o LowTatras\dtm_1m_hillshaded.png

We thank the Geodesy, Cartography and Cadastre Authority of the Slovak Republic for providing their LiDAR as open data for both commercial and non-commercial purposes and name the source of the data used above (as the license requires) as the Office of Geodesy, Cartography and Cadastre of the Slovak Republic (GCCA SR) or – in Slovak – the Úrad geodézie, kartografie a katastra Slovenskej republiky (ÚGKK SR).

Which European country goes next? Czech Republic? Poland? Hungary? Switzerland?



Completeness and Correctness of Discrete LiDAR Returns per Laser Pulse fired

Again and again we have preached about the importance of quality checking when you first get your expensive LiDAR data from the vendor or your free LiDAR data from an open data portal. The minimal quality check we usually advocate consists of lasinfo, lasvalidate, lasoverlap, and lasgrid. The information computed by these LAStools can reassure you that the data contains the right information, is specification conform, has properly aligned flight lines, and has the density distribution you expect. For deliveries or downloads in LAZ format we in addition recommend running laszip with the option ‘-check’ to find the rare file that might have gotten bit-corrupted or truncated during the transfer or the download. Today we learn about a more advanced quality check that can be done by running lassort followed by lasreturn.

For every laser shot fired there are usually between one to five discrete LiDAR returns and some full-waveform systems may even deliver up to fifteen returns. Each of these one to fifteen returns is then given the exact same GPS time stamp that corresponds to the moment in time the laser pulse was fired. By having these unique GPS time stamps we can always recover the set of returns that come from the same laser shot. This makes it possible to check completeness (are all the returns in the file) and correctness (is the returns numbering correct) for the discrete returns of each laser pulse.


Showing all sets of returns in the file that do not have an unique GPS time stamp because the set has one or more duplicate returns (e.g. two first returns, two second returns, … ).

With LAStools we can do this by running lassort followed by lasreturn for any LiDAR that comes from a single beam system. For LiDAR that comes from some multi-beam system, such as the Velodyne 16, 32, 64, or 128, the Optech Pegasus, the RIEGL LMS 1560 (aka “crossfire”), or the Leica ALS70 or ALS80 we first need to seperate the files into one file per beam, which can be done with lassplit.  In the following we investigate data coming from an Optech Galaxy single-beam system. First we sort the returns by GPS time stamp using lassort (this step can be omitted if the data is already sorted in acquisition order (aka by increasing GPS time stamps)) and then we check the return numbering with lasreturn:

lassort -i L001-1-M01-S1-C1_r.laz -gps_time -odix _sorted -olaz

lasreturn -i L001-1-M01-S1-C1_r_sorted.laz -check_return_numbering
checked returns of 11809046 multi and 8585573 single return pulses. took 26.278 secs
missing: 0 duplicate: 560717 too large: 0 zero: 0
200543 returns with n = 1 and r = 1 are duplicate
80548 returns with n = 2 and r = 1 are duplicate
80548 returns with n = 2 and r = 2 are duplicate
41962 returns with n = 3 and r = 1 are duplicate
41962 returns with n = 3 and r = 2 are duplicate
41962 returns with n = 3 and r = 3 are duplicate
13753 returns with n = 4 and r = 1 are duplicate
13753 returns with n = 4 and r = 2 are duplicate
13753 returns with n = 4 and r = 3 are duplicate
13753 returns with n = 4 and r = 4 are duplicate
3636 returns with n = 5 and r = 1 are duplicate
3636 returns with n = 5 and r = 2 are duplicate
3636 returns with n = 5 and r = 3 are duplicate
3636 returns with n = 5 and r = 4 are duplicate
3636 returns with n = 5 and r = 5 are duplicate
WARNING: there are 59462 GPS time stamps that have returns with different number of returns

The output we see above indicates a problem in the return numbering. A recently added new options to lasreturn that allow to reclassify those returns that seem to be part of a problematic set of returns that either contains missing returns, duplicate returns, or returns with different values for the “numbers of returns of given pulse” attribute. This allows us to visualize the issue with lasview. All returns whose are part of a problematic set is shown in the image above.

lasreturn -i L001-1-M01-S1-C1_r_sorted.laz ^
          -check_return_numbering ^
          -classify_as 8 ^
          -classify_duplicate_as 9 ^
          -classify_violation_as 7 ^
          -odix _marked -olaz

This command will mark all sets of returns (i.e. returns that have the exact same GPS time stamp) that have missing returns as 8, that have duplicate returns as 9, and that have returns which different “number of returns per pulse” attribute as 7. The data we have here has no missing returns (no returns are classified as 8) but we have duplicate (9) and violating (7) returns. We look at them closely in single scan lines to conclude.

It immediately becomes obvious that the same GPS time stamp was assigned to the returns of pair of subsequent shots. If the subsequent shots have the same number of returns per shot they are classified as duplicate (9 or blue). If the subsequent shots have different number of returns per shot they are marked as violating (7 or violett) but the reason for the issue is the same. We can look at a few of these return sets in ASCII. Here two subsequent four return shots that have the same GPS time stamp.

237881.011730 4 1 691602.736 5878246.425 141.992 6 79
237881.011730 4 2 691602.822 5878246.415 141.173 6 89
237881.011730 4 3 691603.051 5878246.389 138.993 6 44
237881.011730 4 4 691603.350 5878246.356 136.150 6 169
237881.011730 4 1 691602.793 5878246.439 142.037 6 114
237881.011730 4 2 691602.883 5878246.429 141.185 6 96
237881.011730 4 3 691603.109 5878246.404 139.033 6 50
237881.011730 4 4 691603.414 5878246.370 136.129 6 137

Here a four return shot followed by a three return shot that have the same GPS time stamp.

237881.047753 4 1 691603.387 5878244.501 140.187 6 50
237881.047753 4 2 691603.602 5878244.476 138.141 6 114
237881.047753 4 3 691603.776 5878244.456 136.490 6 60
237881.047753 4 4 691603.957 5878244.436 134.767 6 116
237881.047753 3 1 691603.676 5878244.492 138.132 6 97
237881.047753 3 2 691603.845 5878244.473 136.534 6 90
237881.047753 3 3 691604.034 5878244.452 134.739 6 99

It appears the GPS time counter in the LMS export software did not store the GPS time with sufficient resolution to always distinguish subsequent shots. The issue was confirmed by Optech and was already fixed a few months ago.

We should point out that these double-used GPS time stamps have zero impact on the geometric quality of the point cloud or the distribution of returns. The drawback is that not all returns can easily be grouped into one unique set per laser shot and that the files are not entirely specification conform. Any software that relies on accurate and unique GPS time stamps (such as flight line alignment software) may potentially struggle as well. The bug of the twice-used GPS time stamps was a discovery that is probably of such low consequence that no user of Optech Galaxy data had noticed it in the 4 years that Galaxy had been sold … until we really really scrutinized some data from one of our clients. Optech reports that the issue has been fixed now. But there are other vendors out there with even more serious GPS time and return numbering issues … to be continued.