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

Gabriele Garnero (recipient of three LASmoons)
Interuniversity Department of Regional and Urban Studies and Planning
Politecnico e Università degli Studi, Torino

Last spring, the LARTU research group produced a laser scanner survey of the Abbey of Sant’Andrea in Vercelli, on the occasion of the VIII centenary of the dedication (1219). The database produced with a topographic tool that integrates the potential of a total station with laser scanner and photogrammetric sensors (Trimble SX 10), has been used to produce representations that can be consulted in interactive mode, navigating within the point clouds and producing a consultation platform that can also be accessed by non-specialist users such as art historians or archaeologists.


The LAStools software will be used to improve both the point cloud produced by eliminating the remaining noises, and check other ways of publishing the data, so as to make it usable from outside, to the community of researchers.

laser scanner and photogrammetric acquisitions of the interior of the building (150 millions of points)
+ laser scanner and photogrammetric acquisitions of the outside of the building (210 millions of points)
+ drone-based shooting of outdoor areas processed with Pix4D (23 millions of points)

LAStools processing:
1) tile large point cloud into tiles with buffer [lastile]
2) mark set of points whose z coordinate is a certain percentile of that of their neighbors [lasthin]
3) remove isolated low points from the set of marked points [lasnoise]
4) classify marked points into ground and non-ground [lasground]
5) creates a LiDAR portal for 3D visualization of LAS files [laspublish]

City of Guadalajara creates first Open LiDAR Portal of Latin America

Small to medium sized LiDAR data sets can easily be published online for exploration and download with laspublish of LAStools, which is an easy-to-use wrapper around the powerful Potree open source software for which rapidlasso GmbH has been a major sponsor. During a workshop on LiDAR processing at CICESE in Ensenada, Mexico we learned that Guadalajara – the city with five “a” in its name – has recently published its LiDAR holdings online for download using an interactive 3D portal based on Potree.

There is a lot more data available in Mexico but only Guadalajara seems to have an interactive download portal at the moment with open LiDAR. Have a look at the map below to get an idea of the LiDAR holdings that are held in the archives of the Instituto Nacional de Estadística y Geografía (INEGI). You can request this data either by filling out this form or by sending an email to You will need to explain the use of the information, but apparently INEGI has a fast response time. I was given the KML files you see below and told that each letter in scale 1: 50,000 is divided into 6 regions (a-f) and each region subdivided into 4 parts. Contact me if you want the KML files or if you can provide further clarification on this indexing scheme and/or the data license.

LiDAR available at the Instituto Nacional de Estadística y Geografía (INEGI)

But back to Guadalajara’s open LiDAR. The tile names become visible when you zoom in closer on the map with the tiling overlay as seen below. An individual tile can easily be downloaded by first clicking so that it becomes highlighted and then pressing the “D” button in the lower left corner. We download the two tiles called ‘F08C04.laz’ and ‘F08C05.laz’ and use lasinfo to determine that their average density is 9.0 and 8.9 last returns per per square meter. This means on average 9 laser pulses were fired at each square meter in those two tiles.

lasinfo -i F08C04.laz -cd
lasinfo -i F08C05.laz -cd

Selecting a tile on the map and pressing the “D” button will download the highlighted tile.

The minimal quality check that we recommend doing for any newly obtained LiDAR data is to verify proper alignment of the flightlines using lasoverlap. For tiles with properly populated ‘point source ID’ fields this can be done using the command line shown below.

lasoverlap -i F08C04.laz F08C05.laz ^
           -min_diff 0.1 -max_diff 0.3 ^
           -odir quality -opng ^
           -cores 2

We notice some slight miss-alignments in the difference image (see other tutorials such as this one for how to interpret the resulting color images). We suggest you follow the steps done there to take a closer look at some of the larger strip-like areas that exhibit some systematic disscolorization (compared to other areas) into overly blueish or reddish tones of with lasview. Overlaying one of the resulting *_diff.png files in the GUI of LAStools makes it easy to pick a suspicious area.

We use the “pick” functionality to view only the building of interest.

Unusual are also the large red and blue areas where some of the taller buildings are. Usually those are just one pixel wide which has to do with the laser of one flightline not being able to see the lower area seen by the laser of the other flightline because the line-of-sight is blocked by the structure. We have a closer look at one of these unusual building colorization by picking the building shown above and viewing it with the different visualization options that are shown in the images below.

No. Those are not the “James Bond movie” kind of lasers that burn holes into the building to get ground returns through several floors. The building facade is covered with glass so that the lasers do not scatter photons when they hit the side of the building. Instead they reflect by the usual rule “incidence angle equals reflection angle” of perfectly specular surfaces and eventually hit the ground next to the building. Some of the photons travel back the same way to the receiver on the plane where they get registered as returns. The LiDAR system has no way to know that the photons did not travel the usual straight path. It only measures the time until the photons return and generates a return at the range corresponding to this time along the direction vector that this laser shot was fired at. If the specular reflection of the photons hits a truck or a tree situated next to to building, then we should find that truck or that tree – mirrored by the glossy surface of the building – on the inside of the building. If you look careful at the “slice” through the building below you may find an example … (-:

Some objects located outside the building are mirrored into the building due to its glossy facade.

Kudos to the City of Guadalajara for becoming – to my knowledge – the first city in Latin America to both open its entire LiDAR holdings and also making it available for download in form of a nice and functional interactive 3D portal.

LASmoons: Marzena Wicht

Marzena Wicht (recipient of three LASmoons)
Department of Photogrammetry, Remote Sensing and GIS
Warsaw University of Technology, Poland.

More than half of human population (Heilig 2012) suffers from many negative effects of living in cities: increased air pollution, limited access to the green areas, Urban Heat Island (UHI) and many more. To mitigate some of these effects, many ideas came up over the years: reducing the surface albedo, the idea of the Garden City, green belts, and so on. Increasing horizontal wind speed might actually improve both, the air pollution dispersion and the thermal comfort in urban areas (Gál & Unger 2009). Areas of low roughness promote air flow – discharging the city from warm, polluted air and supplying it with cool and fresh air – if they share specific parameters, are connected and penetrate the inner city with a country breeze. That is why mapping low roughness urban areas is important in better understanding urban climate.

The goal of this study is to derive buildings (outlines and height) and high vegetation using LAStools and to use that data in mapping urban ventilation corridors for our case study area in Warsaw. There are many ways to map these; however using ALS data has certain advantages (Suder& Szymanowski 2014) in this case: DSMs can be easily derived, tree canopy (incl. height) can be joined to the analysis and buildings can be easily extracted. The outputs are then used as a basis for morphological analysis, like calculating frontal area index. LAStools has the considerable advantage of processing large quantities of data (~500 GB) efficiently.

Frontal area index calculation based on 3D building database

+ LiDAR provided by Central Documentation Center of Geodesy and Cartography
+ average pulse density 12 p/m^2
+ covers 517 km^2 (whole Warsaw)

LAStools processing:
1) quality checking of the data as described in several videos and blog posts [lasinfo, lasvalidate, lasoverlap, lasgrid, lasduplicate, lasreturnlas2dem]
2) reorganize data into sufficiently small tiles with buffers to avoid edge artifacts [lastile]
3) classify point clouds into vegetation and buildings [lasground, lasclassify]
4) normalize LiDAR heights [lasheight]
5) create triangulated, rasterized derivatives: DSM / DTM / nDSM / CHM [las2dem, blast2dem]
6) compute height-based metrics (e.g. ‘-avg’, ‘-std’, and ‘-p 50’) [lascanopy]
7) generate subsets during the workflow [lasclip]
8) generate building footprints [lasboundary]

Heilig, G. K. (2012). World urbanization prospects: the 2011 revision. United Nations, Department of Economic and Social Affairs (DESA), Population Division, Population Estimates and Projections Section, New York.
Gal, T., & Unger, J. (2009). Detection of ventilation paths using high-resolution roughness parameter mapping in a large urban area. Building and Environment, 44(1), 198-206.
Suder, A., & Szymanowski, M. (2014). Determination of ventilation channels in urban area: A case study of Wroclaw (Poland). Pure and Applied Geophysics, 171(6), 965-975.

Leaked: “Classified LiDAR” of Pentagon in LAS 1.4 Format

LiDAR leaks have happened! Black helicopters are in the sky!  A few days ago a tiny tweet leaked the online location of “classified LiDAR” for Washington, DC. This LiDAR really is “classified” and includes an aerial scan of the Pentagon. For rogue scientists world-wide we offer a secret download link. It links to a file code-named ‘pentagon.laz‘ that contains the 8,044,789 “classified” returns of the Pentagon shown below. This “classified file” can be deciphered by any software with native LAZ support. It was encrypted with the “LAS 1.4 compatibility mode” of LASzip. The original LAS 1.4 content was encoded into a inconspicuous-looking LAZ file. New point attributes (such as the scanner channel) were hidden as “extra bytes” for fully lossless encryption. Use ‘laszip‘ to fully decode the original “classified” LAS 1.4 file … (-;

Seriously, a tiled LiDAR data set for the District of Columbia flown in 2015 is available for anyone to use on Amazon S3 with a very permissive open data license, namely the Creative Commons Attribution 3.0 License. The LiDAR coverage can be explored via this interactive map. The tiles are provided in LAS 1.4 format and use the new point type 6. We downloaded a few tiles near the White House, the Capitol, and the Pentagon to test the “native LAS 1.4 extension” of our LASzip compressor which will be released soon (a prototype for testing is already available). As these uncompressed LAS files are YUUUGE we use the command line utility ‘wget‘ for downloading. With option ‘-c’ the download continues where it left off in case the transfer gets interrupted.

LiDAR pulse density from 20 or less (blue) to 100 or more (red) pulses per square meter.

We use lasboundary to create labeled bounding boxes for display in Google Earth and lasgrid to a create false color visualization of pulse density with the command lines shown below. Pulse densities of 20 or below are mapped to blue. Pulse densities of 100 or above are mapped to red. We picked the min value 20 and the max value 100 for this false color mapping by running lasinfo with the ‘-cd’ option to compute an average pulse density and then refining the numbers experimentally. We also use lasoverlap to visualize how flightlines overlap and how well they align. Vertical differences of up to 20 cm are mapped to white and differences of 40 cm or more are mapped to saturated blue or red.

lasboundary -i *.las ^
            -use_bb ^
            -labels ^
            -odir quality -odix _bb -okml

lasgrid -i *.las ^
        -keep_last ^
        -point_density -step 2 ^
        -false -set_min_max 20 100 ^
        -odir quality -odix _d_20_100 -opng ^
        -cores 2

lasoverlap -i *.las ^
           -min_diff 0.2 -max_diff 0.4 ^
           -odir quality -opng ^
           -cores 2

The visualization of the pulse density and of the flightline overlap both show that there is no LiDAR for the White House or Capitol Hill. We will never know how tall the tomato and kale plants had grown in Michelle Obama’s organic garden on that day. Note that the White House and Capitol Hill were not simply “cut out”. Instead the flight plan of the survey plane was carefully designed to avoid these areas. Surprisingly, the Pentagon did not receive the same treatment and is (almost) fully included in the open LiDAR as mentioned in the dramatic first paragraph. Interesting is how the varying (tidal?) water level of the Potomac River shows up in the visualization of flightline miss-alignments.

There are a number of issues in these LiDAR files. The most serious ones are reported at the very end of this article. We will now scrutinize the partly-filled tile 2016.las close to the White House with only 11,060,334 returns. A lasvalidate check immediately reports three deviations from the LAS 1.4 specification:

lasvalidate -i 2016.las -o 2016_check.xml
  1. For proper LAS 1.4 files containing point type 6 through 10 all ‘legacy’ point counts in the LAS header should be set to 0. The following six fields in the LAS header should be zero for tile 2016.las (and all other tiles):
    + legacy number of point records
    + legacy number of points by return[0]
    + legacy number of points by return[1]
    + legacy number of points by return[2]
    + legacy number of points by return[3]
    + legacy number of points by return[4]
  2. There should not be any LiDAR return in a valid LAS file whose ‘number of returns of given pulse’ attribute is zero but there are 8 such points in tile 2016.las (and many more in various other tiles).
  3. There should not be any LiDAR return whose ‘return number’ attribute is larger than their ‘number of returns of given pulse’ attribute but there are 8 such points in tile 2016.las (and many more in various other tiles).

The first issue is trivial. There is an efficient in-place fix that does not require to rewrite the entire file using lasinfo with the following command line:

lasinfo -i 2016.las ^
        -nh -nv -nc ^
        -set_number_of_point_records 0 ^
        -set_number_of_points_by_return 0 0 0 0 0 ^

A quick check with las2txt shows us that the second and third issue are caused by the same eight points. Instead of writing an 8 for the ‘number of returns’ attribute the LAS file exporter must have written a 0 (marked in red for all eight returns) and instead of writing an 8 for the ‘return number’ attribute the LAS file exporter must have written a 1 (also marked in red). We can tell it from the true first return via its z coordinate (marked in blue) as the last return should be the lowest of all.

las2txt -i 2016.las ^
        -keep_number_of_returns 0 ^
        -parse xyzrnt ^
397372.70 136671.62 33.02 4 0 112813299.954811
397372.03 136671.64 28.50 5 0 112813299.954811
397371.28 136671.67 23.48 6 0 112813299.954811
397370.30 136671.68 16.86 7 0 112813299.954811
397369.65 136671.70 12.50 1 0 112813299.954811
397374.37 136671.58 44.17 3 0 112813299.954811
397375.46 136671.56 51.49 1 0 112813299.954811
397374.86 136671.57 47.45 2 0 112813299.954811

With las2las we can change the ‘number of returns’ from 0 to 8 using a ‘-filtered_transform’ as illustrated in the command line below. We suspect that higher number of returns such as 9 or 10 might have been mapped to 1 and 2. Fixing those as well as repairing the wrong return numbers will require a more complex tool. We would recommend to check all tiles with more scrutiny using the lasreturn tool. But wait … more return numbering issues are to come.

las2las -i 2016.las ^
        -keep_number_of_returns 0 ^
        -filtered_transform ^
        -set_extended_number_of_returns 8 ^
        -odix _fixed -olas

A closer look at the scan pattern reveals that the LiDAR survey was flown with a dual-beam system where two laser beams scan the terrain simultaneously. This is evident in the textual representation below as there are multiple “sets of returns” for the same GPS time stamp such as 112813952.110394. We group the returns from the two beams into an orange and a green group. Their coordinates show that the two laser beams point into different directions when they are simultaneously “shot” and therefore hit the terrain far apart from another.

las2txt -i 2016.las ^
        -keep_gps_time 112813952.110392 112813952.110396 ^
        -parse xyzlurntp ^
397271.40 136832.35 54.31 0 0 1 1 112813952.110394 117
397277.36 136793.35 38.68 0 1 1 4 112813952.110394 117
397277.35 136793.56 32.89 0 1 2 4 112813952.110394 117
397277.34 136793.88 24.13 0 1 3 4 112813952.110394 117
397277.32 136794.25 13.66 0 1 4 4 112813952.110394 117

The information about which point is from which beam is currently stored into the generic ‘user data’ attribute instead of into the dedicated ‘scanner channel’ attribute. This can be fixed with las2las as follows.

las2las -i 2016.las ^
        -copy_user_data_into_scanner_channel ^
        -set_user_data 0 ^
        -odix _fixed -olas

Unfortunately the LiDAR files have much more serious issues in the return numbering. It’s literally a “Total Disaster!” and “Sad!” as the US president will tweet shortly. After grouping all returns with the same GPS time stamp into an orange and a green group there is one more set of returns left unaccounted for.

las2txt -i 2016.las ^
        -keep_gps_time 112813951.416451 112813951.416455 ^
        -parse xyzlurntpi ^
397286.02 136790.60 45.90 0 0 1 4 112813951.416453 117 24
397286.06 136791.05 39.54 0 0 2 4 112813951.416453 117 35
397286.10 136791.51 33.34 0 0 3 4 112813951.416453 117 24
397286.18 136792.41 21.11 0 0 4 4 112813951.416453 117 0
397286.12 136791.75 30.07 0 0 1 1 112813951.416453 117 47
397291.74 136750.70 45.86 0 1 1 1 112813951.416453 117 105
las2txt -i 2016.las ^
        -keep_gps_time 112813951.408708 112813951.408712 ^
        -parse xyzlurntpi ^
397286.01 136790.06 45.84 0 0 1 4 112813951.408710 117 7
397286.05 136790.51 39.56 0 0 2 4 112813951.408710 117 15
397286.08 136790.96 33.33 0 0 3 4 112813951.408710 117 19
397286.18 136792.16 17.05 0 0 4 4 112813951.408710 117 0
397286.11 136791.20 30.03 0 0 1 2 112813951.408710 117 58
397286.14 136791.67 23.81 0 0 2 2 112813951.408710 117 42
397291.73 136750.16 45.88 0 1 1 1 112813951.408710 117 142

This can be visualized with lasview and the result is unmistakably clear: The return numbering is messed up. There should be one shot with five returns (not a group of four and a single return) in the first example. And there should be one shot with six returns (not a group of four and a group of two returns) in the second example. Such a broken return numbering results in extra first (or last) returns. These are serious issues that affect any algorithm that relies on the return numbering such as first-return DSM generation or canopy cover computation. Those extra returns will also make the pulse density appear higher and the pulse spacing appear tighter than they really are. The numbers from 20 (blue) to 100 (red) pulses per square meters in our earlier visualization are definitely inflated.

lasview -i 2016.las ^
        -keep_gps_time 112813951.416451 112813951.416455 ^

lasview -i 2016.las ^
        -keep_gps_time 112813951.408708 112813951.408712 ^

After all these troubles here something nice. Side-by-side a first-return TIN and a spike-free TIN (using a freeze of 0.8 m) of the center court cafe in the Pentagon. Especially given all these “fake first returns” in the Washington DC LiDAR we really need the spike-free algorithm to finally “Make a DSM great again!” … (-;

We would like to acknowledge the District of Columbia Office of the Chief Technology Officer (OCTO) for providing this data with a very permissive open data license, namely the Creative Commons Attribution 3.0 License.


NRW Open LiDAR: Merging Points into Proper LAS Files

In the first part of this series we downloaded, compressed, and viewed some of the newly released open LiDAR data for the state of North Rhine-Westphalia. In the second part we look at how to merge the multiple point clouds provided back into single LAS or LAZ files that are as proper as possible. Follow along with a recent version of LAStools and a pair of DGM and DOM files for your area of interest. For downloading the LiDAR we suggest using the wget command line tool with option ‘-c’ that after interruption in transmission will restart where it left off.

In the first part of this series we downloaded the pair of DGM and DOM files for the City of Bonn. The DGM file and the DOM file are zipped archives that contain the points in 1km by 1km tiles stored as x, y, z coordinates with centimeter resolution. We had already converted these textual *.xyz files into binary *.laz files. We did this with the open source LASzip compressor that is distributed with LAStools as described in that blog post. We continue now with the assumption that you have converted all of the *.xyz files to *.laz files as described here.

Mapping from tile names of DGM and DOM archives to classification and return type of points.

The mapping from tile names in DGM and DOM archives to the classification and return type of points: lp = last return. fp = first return, ab,aw,ag = synthetic points

There are multiple tiles for each square kilometer as the LiDAR has been split into different files based on classification and return type. Furthermore there are also synthetic points that were created by the land survey department to replace LiDAR under bridges and along buildings for generating higher quality rasters. We want to combine all points of a square kilometer into a single LAZ tile as it is usually expected. Simply merging the multiple files per tile is not an option as this would result in loosing point classifications and return type information as well as in duplicating all single returns that are stored in more than one file. The folks at OpenNRW offer this helpful graphic to know what the acronyms above mean:

Illustration of how acronyms used in tile names correspond to point classification and type.

Illustration of how acronyms used in tile names correspond to point classification and type.

In the following we’ll be looking at the set of files corresponding to the UTM tile 32366 / 5622. We wanted an interesting area with large buildings, a bridge, and water. We were looking for a suitable area using the KML overlays generated in part one. The tile along the Rhine river selected in the picture below covers the old city hall, the opera house, and the “Kennedy Bridge” has a complete set of DGM and DOM files:

      3,501 dgm1l-ab_32366_5622_1_nw.laz
     16,061 dgm1l-ag_32366_5622_1_nw.laz
      3,269 dgm1l-aw_32366_5622_1_nw.laz
    497,008 dgm1l-brk_32366_5622_1_nw.laz
  7,667,715 dgm1l-lpb_32366_5622_1_nw.laz
 12,096,856 dgm1l-lpnb_32366_5622_1_nw.laz
     15,856 dgm1l-lpub_32366_5622_1_nw.laz

      3,269 dom1l-aw_32366_5622_1_nw.laz
 21,381,106 dom1l-fp_32366_5622_1_nw.laz

We find the name of the tiles that cover the "Kennedy Bridge" using the KML overlays generated in part one.

We find the name of the tile that covers the “Kennedy Bridge” using the KML overlays generated in part one.

We now assign classification codes and flags to the returns from the different files using las2las, merge them together with lasmerge, and recover single, first, and last return information with lasduplicate. We set classifications to bridge deck (17), ground (2), to unclassified (1), and to noise (7) for all returns in the files with the acronym ‘brk’ (= bridge points), the acronym ‘lpb’ (= last return ground), the acronym ‘lpnb’ (= last return non-ground), and the acronym ‘lpub’ (= last return under ground). with las2las and store the resulting files to a temporary folder.

las2las -i dgm1l-brk_32366_5622_1_nw.laz ^
        -set_classification 17 ^
        -odir temp -olaz

las2las -i dgm1l-lpb_32366_5622_1_nw.laz ^
        -set_classification 2 ^
        -odir temp -olaz

las2las -i dgm1l-lpnb_32366_5622_1_nw.laz ^
        -set_classification 1 ^
        -odir temp -olaz

las2las -i dgm1l-lpub_32366_5622_1_nw.laz ^
        -set_classification 7 ^
        -odir temp -olaz

Next we use the synthetic flag of the LAS format specification to flag any additional points that were added (no measured) by the survey department to generate better raster products. We set classifications to ground (2) and the synthetic flag for all points of the files with the acronym ‘ab’ (= additional ground) and the acronym ‘ag’ (= additional building footprint). We set classifications to water (9) and the synthetic flag for all points of the files with the acronym ‘aw’ (= additional water bodies). Files with acronym ‘aw’ appear both in the DGM and DOM archive. Obviously we need to keep only one copy.

las2las -i dgm1l-ab_32366_5622_1_nw.laz ^
        -set_classification 2 ^
        -set_synthetic_flag 1 ^
        -odir temp -olaz

las2las -i dgm1l-ag_32366_5622_1_nw.laz ^
        -set_classification 2 ^
        -set_synthetic_flag 1 ^
        -odir temp -olaz

las2las -i dgm1l-aw_32366_5622_1_nw.laz ^
        -set_classification 9 ^
        -set_synthetic_flag 1 ^
        -odir temp -olaz

Using lasmerge we merge all returns from files with acronyms ‘brk’ (= bridge points), ‘lpb’ (= last return ground),  ‘lpnb’ (= last return non-ground), and ‘lpub’ (= last return under ground) into a single file that will then contain all of the (classified) last returns for this tile.

lasmerge -i temp\dgm1l-brk_32366_5622_1_nw.laz ^
         -i temp\dgm1l-lpb_32366_5622_1_nw.laz ^
         -i temp\dgm1l-lpnb_32366_5622_1_nw.laz ^
         -i temp\dgm1l-lpub_32366_5622_1_nw.laz ^
         -o temp\dgm1l-lp_32366_5622_1_nw.laz

Next we run lasduplicate three times to recover which points are single returns and which points are the first and the last return of a pair of points generated by the same laser shot. First we run lasduplicate with option ‘-unique_xyz’ to remove any xyz duplicates from the last return file. We also mark all surviving returns as the second of two returns. Similarly, we remove any xyz duplicates from the first return file and mark all survivors as the first of two returns. Finally we run lasduplicate with option ‘-single_returns’ with the unique last and the unique first return files as ‘-merged’ input. If a return with the exact same xyz coordinates appears in both files only the first copy is kept and marked as a single return. In order to keep the flags and classifications from the last return file, the order in which the last and first return files are listed in the command line is important.

lasduplicate -i temp\dgm1l-lp_32366_5622_1_nw.laz ^
             -set_return_number 2 -set_number_of_returns 2 ^
             -unique_xyz ^
             -o temp\last_32366_5622_1_nw.laz

lasduplicate -i dom1l-fp_32366_5622_1_nw.laz ^
             -set_return_number 1 -set_number_of_returns 2 ^
             -unique_xyz ^
             -o temp\first_32366_5622_1_nw.laz

lasduplicate -i temp\last_32366_5622_1_nw.laz ^
             -i temp\first_32366_5622_1_nw.laz ^
             -merged ^
             -single_returns ^
             -o temp\all_32366_5622_1_nw.laz

We then add the synthetic points with another call to lasmerge to obtain a LAZ file containing all points of the tile correctly classified, flagged, and return-numbered.

lasmerge -i temp\dgm1l-ab_32366_5622_1_nw.laz ^
         -i temp\dgm1l-ag_32366_5622_1_nw.laz ^
         -i temp\dgm1l-aw_32366_5622_1_nw.laz ^
         -i temp\all_32366_5622_1_nw.laz ^
         -o temp\merged_32366_5622_1_nw.laz

Optional: For more efficient use of this file in subsequent processing – and especially to accelerate area-of-interest queries with lasindex – it is often of great advantage to reorder the points in a spatially coherent manner. A simple call to lassort will rearrange the points along a space-filling curve such as a Hilbert curve or a Z-order curve.

lassort -i temp\merged_32366_5622_1_nw.laz ^
        -o bonn_32366_5622_1_nw.laz

Note that we also renamed the file because a good name can be useful if you find that file again in two years from now. Let’s have a look at the result with lasview.

lasview -i bonn_32366_5622_1_nw.laz

In lasview you can press <c> repeatedly to switch through all available coloring modes until you see the yellow (single) / red (first) / last (blue) coloring of the returns. The recovered return types are especially evident under vegetation, in the presence of wires, and along building edges. Press <x> to select an area of interest and press <x> again to inspect it more closely. Press <i> while hovering above a point to show its coordinates, classification, and return type.

We did each processing in separate steps to illustrate the overall workflow. The above sequence of LAStools command line calls can be shortened by combining multiple processing steps into one operation. This is left as an exercise for the advanced LAStools user … (-;

Acknowledgement: The LiDAR data of OpenNRW comes with a very permissible license. It is called “Datenlizenz Deutschland – Namensnennung – Version 2.0” or “dl-de/by-2-0” and allows data and derivative sharing as well as commercial use. It only requires us to name the source. We need to cite the “Land NRW (2017)” with the year of the download in brackets and specify the Universal Resource Identification (URI) for both the DOM and the DGM. Done. So easy. Thank you, OpenNRW … (-:

NRW Open LiDAR: Download, Compression, Viewing

UPDATE: (March 6th): Second part merging Bonn into proper LAS files

This is the first part of a series on how to process the newly released open LiDAR data for the entire state of North Rhine-Westphalia that was announced a few days ago. Again, kudos to OpenNRW for being the most progressive open data state in Germany. You can follow this tutorial after downloading the latest version of LAStools as well as a pair of DGM and DOM files for your area of interest from these two download pages.

We have downloaded the pair of DGM and DOM files for the Federal City of Bonn. Bonn is the former capital of Germany and was host to the FOSS4G 2016 conference. As both files are larger than 10 GB, we use the wget command line tool with option ‘-c’ that will restart where it left off in case the transmission gets interrupted.

The DGM file and the DOM file are zipped archives that contain the points in 1km by 1km tiles stored as x, y, z coordinates in ETRS89 / UTM 32 projection as simple ASCII text with centimeter resolution (i.e. two decimal digits).

>> more
360000.00 5613026.69 164.35
360000.00 5613057.67 164.20
360000.00 5613097.19 164.22
360000.00 5613117.89 164.08
360000.00 5613145.35 164.03

There is more than one tile for each square kilometer as the LiDAR points have been split into different files based on their classification and their return type. Furthermore there are also synthetic points that were used by the land survey department to replace certain LiDAR points in order to generate higher quality DTM and DSM raster products.

The zipped DGM archive is 10.5 GB in size and contains 956 *.xyz files totaling 43.5 GB after decompression. The zipped DOM archive is 11.5 GB in size and contains 244 *.xyz files totaling 47.8 GB. Repeatedly loading these 90 GB of text data and parsing these human-readable x, y, and z coordinates is inefficient with common LiDAR software. In the first step we convert the textual *.xyz files into binary *.laz files that can be stored, read and copied more efficiently. We do this with the open source LASzip compressor that is distributed with LAStools using these two command line calls:

laszip -i dgm1l_05314000_Bonn_EPSG5555_XYZ\*.xyz ^
       -epsg 25832 -vertical_dhhn92 ^
       -olaz ^
       -cores 2
laszip -i dom1l_05314000_Bonn_EPSG5555_XYZ\*.xyz ^
       -epsg 25832 -vertical_dhhn92 ^
       -olaz ^
       -cores 2

The point coordinates are is in EPSG 5555, which is a compound datum of horizontal EPSG 25832 aka ETRS89 / UTM zone 32N and vertical EPSG 5783 aka the “Deutsches Haupthoehennetz 1992” or DHHN92. We add this information to each *.laz file during the LASzip compression process with the command line options ‘-epsg 25832’ and ‘-vertical_dhhn92’.

LASzip reduces the file size by a factor of 10. The 956 *.laz DGM files compress down to 4.3 GB from 43.5 GB for the original *.xyz files and the 244 *.laz DOM files compress down to 4.8 GB from 47.8 GB. From here on out we continue to work with the 9 GB of slim *.laz files. But before we delete the 90 GB of bulky *.xyz files we make sure that there are no file corruptions (e.g. disk full, truncated files, interrupted processes, bit flips, …) in the *.laz files.

laszip -i dgm1l_05314000_Bonn_EPSG5555_XYZ\*.laz -check
laszip -i dom1l_05314000_Bonn_EPSG5555_XYZ\*.laz -check

One advantage of having the LiDAR in an industry standard such as the LAS format (or its lossless compressed twin, the LAZ format) is that the header of the file stores the number of points per file, the bounding box, as well as the projection information that we have added. This allows us to very quickly load an overview for example, into lasview.

lasview -i dgm1l_05314000_Bonn_EPSG5555_XYZ\*.laz -GUI

The bounding boxes of the DGM files quickly display a preview of the data in the GUI when the files are in LAS or LAZ format.

The bounding boxes of the DGM files quickly give us an overview in the GUI when the files are in LAS or LAZ format.

Now we want to find a particular site in Bonn such as the World Conference Center Bonn where FOSS4G 2016 was held. Which tile is it in? We need some geospatial context to find it, for example, by creating an overview in form of KML files that we can load into Google Earth. We use the files from the DOM folder with “fp” in the name as points on buildings are mostly “first returns”. See what our previous blog post writes about the different file names or look at the second part of this series. We can create the KML files with lasboundary either via the GUI or in the command line.

lasboundary -i dom1l_05314000_Bonn_EPSG5555_XYZ\dom1l-fp*.laz ^

Only the "fp" tiles from the DOM folder loaded the GUI into lasboundary.

Only the “fp” tiles from the DOM folder loaded the GUI into lasboundary.

lasboundary -i dom1l_05314000_Bonn_EPSG5555_XYZ\dom1l-fp*.laz ^
            -use_bb -labels -okml

We zoom in and find the World Conference Center Bonn and load the identified tile into lasview. Well, we did not expect this to happen, but what we see below will make this series of tutorials even more worthwhile. There is a lot of “high noise” in the particular tile we picked. We should have noticed the unusually high z range of 406.42 meters in the Google Earth pop-up. Is this high electromagnetic radiation interfering with the sensors? There are a number of high-tech government buildings with all kind of antennas nearby (such as the United Nations Bonn Campus the mouse cursor points at).

Significant amounts of high noise are in the first returns of the DOM tile we picked.

Significant amounts of high noise are in the first returns of the DOM tile we picked.

But the intended area of interest was found. You can see the iconic “triangulated” roof of the building that is across from the World Conference Center Bonn.

The World Conference Center Bonn is across from the building with the "triangulated" roof.

The World Conference Center Bonn is across from the building with the “triangulated” roof.

Please don’t think it is the responsibility of OpenNRW to remove the noise and provide cleaner data. The land survey has already processed this data into whatever products they needed and that is where their job ended. Any additional services – other than sharing the raw data – are not in their job description. We’ll take care of that … (-:

Acknowledgement: The LiDAR data of OpenNRW comes with a very permissible license. It is called “Datenlizenz Deutschland – Namensnennung – Version 2.0” or “dl-de/by-2-0” and allows data and derivative sharing as well as commercial use. It only requires us to name the source. We need to cite the “Land NRW (2017)” with the year of the download in brackets and specify the Universal Resource Identification (URI) for both the DOM and the DGM. Done. So easy. Thank you, OpenNRW … (-:

Generating Spike-Free Digital Surface Models from LiDAR

A Digital Surface Model (DSM) represents the elevation of the landscape including all vegetation and man-made objects. An easy way to generate a DSM raster from LiDAR is to use the highest elevation value from all points falling into each grid cell. However, this “binning” approach only works when then the resolution of the LiDAR is higher than the resolution of the raster. Only then sufficiently many LiDAR points fall into each raster cell to prevent “empty pixels” and “data pits” from forming. For example, given LiDAR with an average pulse spacing of 0.5 meters one can easily generate a 2.5 meter DSM raster with simple “binning”. But to generate a 0.5 meter DSM raster we need to use an “interpolation” method.

Returns of four fightlines on two trees.

Laser pulses and discrete returns of four fightlines.

For the past twenty or so years, GIS textbooks and LiDAR tutorials have recommened to use only the first returns to construct the interpolating surface for DSM generation. The intuition is that the first return is the highest return for an airborne survey where the laser beams come (more or less) from above. Hence, an interpolating surface of all first returns is constructed – usually based on a 2D Delaunay triangulation – and the resulting Triangular Irregular Network (TIN) is rasterized onto a grid at a user-specified resolution to create the DSM raster. The same way a Canopy Height Model (CHM) is generated except that elevations are height-normalized either before or after the rasterization step. However, using a first-return interpolation for DSM/CHM generation has two critical drawbacks:

(1) Using only first returns means not all LiDAR information is used and some detail is missing. This is particularly the case for off-nadir scan angles in traditional airborne surveys. It becomes more pronounced with new scanning systems such as UAV or hand-held LiDAR where laser beams no longer come “from above”. Furthermore, in the event of clouds or high noise the first returns are often removed and the remaining returns are not renumbered. Hence, any laser shot whose first return reflects from a cloud or a bird does not contribute its highest landscape hit to the DSM or CHM.

(2) Using all first returns practically guarantees the formation of needle-shaped triangles in vegetated areas and along building roofs that appear as spikes in the TIN. This is because at off-nadir scan angles first returns are often generated far below other first returns as shown in the illustration above. The resulting spikes turn into “data pits” in the corresponding raster that not only look ugly but impact the utility of the DSM or CHM in subsequent analysis, for example, in forestry applications when attempting to extract individual trees.

In the following we present results and command-line examples for the new “spike-free” algorithm by (Khosravipour et. al, 2015, 2016) that is implemented (as a slow prototype) in the current LAStools release. This completely novel method for DSM generation triangulates all relevant LiDAR returns using Contrained Delaunay algorithm. This constructs a “spike-free” TIN that is in turn rasterized into “pit-free” DSM or CHM. This work is both a generalization and an improvement of our previous result of pit-free CHM generation.

We now compare our “spike-free” DSM to a “first-return” DSM on the two small urban data sets “france.laz” and “zurich.laz” distributed with LAStools. Using lasinfo with options ‘-last_only’ and ‘-cd’ we determine that the average pulse spacing is around 0.33 meter for “france.laz” and 0.15 meter for “zurich.laz”. We decide to create a hillshaded 0.25 meter DSM for “france.laz” and a 0.15 meter DSM for “zurich.laz” with the command-lines shown below.

las2dem -i ..\data\france.laz ^
        -keep_first ^
        -step 0.25 ^
        -hillshade ^
        -o france_fr.png
las2dem -i ..\data\france.laz ^
        -spike_free 0.9 ^
        -step 0.25 ^
        -hillshade ^
        -o france_sf.png

las2dem -i ..\data\zurich.laz ^
           -keep_first ^
           -step 0.15 ^
           -hillshade ^
           -o zurich_fr.png
las2dem -i ..\data\zurich.laz ^
        -spike_free 0.5 ^
        -step 0.15 ^
        -hillshade ^
        -o zurich_sf.png

The differences between a first-return DSM and a spike-free DSM are most drastic along building roofs and in vegetated areas. To inspect in more detail the differences between a first-return and our spike-free TIN we use lasview that allows to iteratively visualize the construction process of a spike-free TIN.

lasview -i ..\data\france.laz -spike_free 0.9

Pressing <f> and <t> constructs the first-return TIN. Pressing <SHIFT> + <t> destroys the first-return TIN. Pressing <SHIFT> + <y> constructs the spike-free TIN. Pressing <y> once destroys the spike-free TIN. Pressing <y> many times iteratively constructs the spike-free TIN.

One crucial piece of information is still missing. What value should you use as the freeze constraint of the spike-free algorithm that we set to 0.9 for “france.laz” and to 0.5 for “zurich.laz” as the argument to the command-line option ‘-spike_free’. The optimal value is related to the expected edge-length and we found the 99th percentile of a histogram of edge lengths of the last-return TIN to be useful. Or simpler … try a value that is about three times the average pulse spacing.

Khosravipour, A., Skidmore, A.K., Isenburg, M. and Wang, T.J. (2015) Development of an algorithm to generate pit-free Digital Surface Models from LiDAR, Proceedings of SilviLaser 2015, pp. 247-249, September 2015.
Khosravipour, A., Skidmore, A.K., Isenburg, M (2016) Generating spike-free Digital Surface Models using raw LiDAR point clouds: a new approach for forestry applications, (journal manuscript under review).

LASmoons: Andreas Konring and Susanne Bjerg Petersen

Andreas Konring and Susanne Bjerg Petersen (recipients of three LASmoons)
Department of Environmental Engineering
Technical University of Denmark, Lyngby, DENMARK

Copenhagen has in the recent years experienced severe floodings due to cloudbursts which has increased the focus of climate adaption and the implementation of green infrastructure. The use of sustainable urban drainage system (SUDS) solutions to divert stormwater from the existing drainage system will be a central measure to increase the climate resilience while greenifying the city and Copenhagen municipality is investing 700 million euros in SUDS projects alone. Additionally, the city has decided to plant 100.000 new trees in the next 10 years as another measure to enhance natural amenities but also because of air cleansing and cooling effects. However, it has not been investigated what effect the current canopy cover has on the rainwater retention due to increased evaporation and soil infiltration and if planting more trees could help improve the pluvial flooding issues.

Example of a pit-free CHM in an urban environment.

Example of a pit-free CHM in an urban environment.

This study aims to estimate the current number of trees and extract tree metrics such as volume, canopy cover and densities with the use of the national LIDAR dataset and NIR ortophotos from summer and spring. These canopy metrics will be used to inform a simple tree model which will be implemented in a 2-D overland flow model to assess the effect of trees on flood mitigation. The created CHM could also be used in further analysis of the urban heat island effect.

100 square kilometers of the Danish national LiDAR dataset collected in November 2014 covering the municipality of Copenhagen.
+  density of 4 – 5 last-returns per square meter
+  classified into surface (1), ground (2), vegetation (3,4,5), buildings (6), noise (7) and water (9).

LAStools processing:
create square tiles with buffer to avoid edge artifacts [lastile]
2) generate DTMs and DSMs with only buildings and terrain [las2dem]
3).normalize height, remove outliers and keep classes 2, 5 and 6 [lasheight]
4) create rasters with forest metrics [lascanopy]
5) calculate the pit-free Canopy Height Model (CHM) proposed by Khosravipour et al. (2014) [lasthin, las2dem, lasgrid]

Copenhagen Municipality, 2011. Copenhagen Climate Adaption Plan.
Geodatastyrelsen, 2014. Danmarks højdemodel, DHM/Punktsky – Dataversion 2.0 januar 2015. Product specification.
Khosravipour, A., Skidmore, A.K., Isenburg, M., Wang, T.J., Hussin, Y.A., 2014. Generating pit-free Canopy Height Models from Airborne LiDAR. PE&RS = Photogrammetric Engineering and Remote Sensing 80, 863-872.