LASmoons: Manuel Jurado

Manuel Jurado (recipient of three LASmoons)
Departamento de Ingeniería Topográfica y Cartografía
Universidad Politécnica de Madrid, SPAIN

The availability of LiDAR data is creating a lot of innovative possibilities in different fields of science, education, and other field of interests. One of the areas that has been deeply impacted by LiDAR is cartography and in particular one highly specialized sub-field of cartography in the domain of recreational and professional orienteering running: the production of high-quality maps for orienteering races (Ditz et al., 2014). These are thematic maps with a lot of fine detail which demands many hours of field work for the map maker. In order to reduce the fieldwork, LiDAR data obtained from Airborne Research Australia (ARA) is going to be used in order to obtain DEM and to extract features that must be included in these maps. The data will be filtered and processed with the help of LAStools.

Final map with symbolism typical for use in orienteering running

The goal of this project is to extract either point (boulders, mounds), linear (contours, erosion gullies, cliffs) and area features (vegetation density) that should be drawn in a orienteering map derived from high-resolution LiDAR. Different LiDAR derived raster images are being created: 0.5m DTM, vegetation density (J. Ryyppo, 2013), slope, Sky-View factor (Ž. Kokalj et al., 2011), and shaded relief. The area used is in Renmark, South Australia and the produced map is going to be used for the Australian Orienteering Championships 2018.

Sky-View factor of DTM for same area as shown above.

4 square kilometers of airborne LiDAR data produced by Airborne Research Australia at 18 pulses per square meter using the full waveform scanning LiDAR Q680i-S laser scanner from RIEGL
+ 60 hours of check and validation work in the field

LAStools processing:
1) tile into 500 by 500 meter tiles with 20 meter buffer [lastile]
2) classify isolated points as noise [lasnoise]
3) classify point clouds into ground and non-ground [lasground]
4) create a Digital Terrain Model (DTM) [las2dem]
5) normalize height of points above the ground [lasheight]
6) compute vegetation density metrics [lascanopy]
7) create hillshades of the raster DTMs [blast2dem or GDAL]

Ditz, Robert, Franz Glaner, and Georg Gartner. (2014). “Laser Scanning and Orienteering Maps.” Scientific Journal of Orienteering 19.1.
JRyyppo, Jarkko. (2013). “Karttapullautin vegetation mapping guide”.
Kokalj, Žiga, Zaksek, Klemen, and Oštir, Krištof. (2011). Application of sky-view factor for the visualization of historic landscape features in lidar-derived relief models. Antiquity. 85. 263-273.

LASmoons: Chris J. Chandler

Chris J. Chandler (recipient of three LASmoons)
School of Geography
University of Nottingham, UNITED KINGDOM

Wetlands provide a range of important ecosystem services: they store carbon, regulate greenhouse gas emissions, provide flood protection as well as water storage and purification. Preserving these services is critical to achieve sustainable environmental management. Currently, mangrove forests are protected in Mexico, however, fresh water wetland forests, which also have high capacity for storing carbon both in the trees and in the soil, are not protected under present legislation. As a result, coastal wetlands in Mexico are threatened by conversion to grazing areas, drainage for urban development and pollution. Given these threats, there is an urgent need to understand the current state and distribution of wetlands to inform policy and protect the ecosystem services provided by these wetlands.
In this project we will combine field data collection, satellite data (i.e. optical remote sensing, radar and LiDAR remote sensing) and modelling to provide an integrated technology for assessing the value of a range of ecosystem services, tested to proof of concept stage based on carbon storage. The outcome of the project will be a tool for mapping the value of a range of ecosystem services. These maps will be made directly available to local stakeholders including policy makers and land users to inform policy regarding forest protection/legislation and aid development of financial incentives for local communities to protect these services.

Wetland classification in the Chiapas region of Mexico

At this stage of the project we have characterized wetlands for three priority areas in Mexico (Pantanos de Centla, La Encrucijada and La Mancha). Next stage is the up scaling of the field data at the three study sites using LiDAR data for producing high quality Canopy Height Model (CHM), which has been of great importance for biomass estimation (Ferraz et al., 2016). A high quality CHM will be achieved using LAStools software.

LiDAR provided by the Mexican National Institute of Statistics and Geography (INEGI)
+ average height: 5500 m, mirror angle: +/- 30 degrees, speed: 190 knots
+ collected with Cessna 441, Conquest II system at 1 pts/m².

LAStools processing:
create 1000 meter tiles with 35 meter buffer to avoid edge artifacts [lastile]
2) classify point clouds into ground and non-ground [lasground]
3) normalize height of points above the ground [lasheight]
4) create a Digital Terrain and Surface Model (DTM and DSM) [las2dem]
5) generate a spike-free Canopy Height Model (CHM) as described here and here [las2dem]
6) compute various metrics for each plot and the normalized tiles [lascanopy]

Ferraz, A., Saatchi, S., Mallet, C., Jacquemoud S., Gonçalves G., Silva C.A., Soares P., Tomé, M. and Pereira, L. (2016). Airborne Lidar Estimation of Aboveground Forest Biomass in the Absence of Field Inventory. Remote Sensing, 8(8), 653.

Processing Drone LiDAR from YellowScan’s Surveyor, a Velodyne Puck based System

Points clouds from UAVs have become a common sight. Cheap consumer drones equipped with cameras produce points from images with increasing quality as photogrammetry software is improving. But vegetation is always a show stopper for point clouds generated from imagery data. Only an active sensing technique such as laser scanning can penetrate through the vegetation and generate points on the ground under the canopy of a forested area. Advances in UAV technology and the miniaturization of LiDAR systems have allowed lasers-scanning solutions for drones to enter the market.

Last summer we attended the LiDAR for Drone 2017 Conference by YellowScan and processed some data sets flown with their Surveyor system that is built around the Velodyne VLP-16 Puck LiDAR scanner and the Applanix APX15 single board GNSS-Inertial solution. One common challenge observed in LiDAR data generated by the Velodyne Puck is that surfaces are not as “crisp” as those generated by other laser scanners. Flat and open terrain surfaces are described by a layer of points with a “thickness” of a few centimeter as you can see in the images below. This visualization uses a 10 meter by 5 meter cut-out of from this data set with the coordinate range [774280,774290) in x and [6279463,6279468) in y. Standard ground classification routines will “latch onto” the lowermost envelope of these thick point layers and therefore produce a sub-optimal Digital Terrain Model (DTM).

In part this “thickness” can be reduced by using fewer flightlines as the “thickness” of each flightline by itself is lower but it is compounded when merging all flightlines together. However, deciding which (subset of) flightlines to use for which part of the scene to generate the best possible ground model is not an obvious tasks either and even per flightline there will be a remaining “thickness” to deal with as can be seen in the following set of images.

In the following we show how to deal with “thickness” in a layer of points describing a ground surface. We first produce a “lowest ground” which we then widen into a “thick ground” from which we then derive “median ground” points that create a plausible terrain representation when interpolated by a Delaunay triangulation and rasterized onto a DTM. Step by step we process this example data set captured in a “live demo” during the LiDAR for Drone 2017 Conference – the beautiful Château de Flaugergues in Montpellier, France where the event took place. You can download this data via this link if you would like to repeat these processing steps:

Once you decompress the RAR file (e.g. with the UnRar.exe freeware) you will find six raw flight strips in LAS format and the trajectory of the UAV in ASCII text format as it was provided by YellowScan.

E:\LAStools\bin>dir Flaugergues
06/27/2017 08:03 PM 146,503,985 Flaugergues_test_demo_ppk_L1.las
06/27/2017 08:02 PM  91,503,103 Flaugergues_test_demo_ppk_L2.las
06/27/2017 08:03 PM 131,917,917 Flaugergues_test_demo_ppk_L3.las
06/27/2017 08:03 PM 219,736,585 Flaugergues_test_demo_ppk_L4.las
06/27/2017 08:02 PM 107,705,667 Flaugergues_test_demo_ppk_L5.las
06/27/2017 08:02 PM  74,373,053 Flaugergues_test_demo_ppk_L6.las
06/27/2017 08:03 PM   7,263,670 Flaugergues_test_demo_ppk_traj.txt

As usually we start with quality checking by visual inspection with lasview and by creating a textual report with lasinfo.

E:\LAStools\bin>lasview Flaugergues_test_demo_ppk_L1.las

The raw LAS file “Flaugergues_test_demo_ppk_L1.las” colored by elevation.

E:\LAStools\bin>lasinfo Flaugergues_test_demo_ppk_L1.las
lasinfo (171011) report for Flaugergues_test_demo_ppk_L1.las
reporting all LAS header entries:
 file signature: 'LASF'
 file source ID: 1
 global_encoding: 1
 project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
 version major.minor: 1.2
 system identifier: 'YellowScan Surveyor'
 generating software: 'YellowReader by YellowScan'
 file creation day/year: 178/2017
 header size: 227
 offset to point data: 297
 number var. length records: 1
 point data format: 3
 point data record length: 34
 number of point records: 4308932
 number of points by return: 4142444 166488 0 0 0
 scale factor x y z: 0.001 0.001 0.001
 offset x y z: 774282 6279505 92
 min x y z: 774152.637 6279377.623 82.673
 max x y z: 774408.344 6279541.646 116.656
variable length header record 1 of 1:
 reserved 0
 user ID 'LASF_Projection'
 record ID 34735
 length after header 16
 description ''
 GeoKeyDirectoryTag version 1.1.0 number of keys 1
 key 3072 tiff_tag_location 0 count 1 value_offset 2154 - ProjectedCSTypeGeoKey: RGF93 / Lambert-93
reporting minimum and maximum for all LAS point record entries ...
 X -129363 126344
 Y -127377 36646
 Z -9327 24656
 intensity 0 65278
 return_number 1 2
 number_of_returns 1 2
 edge_of_flight_line 0 0
 scan_direction_flag 0 0
 classification 0 0
 scan_angle_rank -120 120
 user_data 75 105
 point_source_ID 1 1
 gps_time 219873.160527 219908.550379
 Color R 0 0
 G 0 0
 B 0 0
number of first returns: 4142444
number of intermediate returns: 0
number of last returns: 4142444
number of single returns: 3975956
overview over number of returns of given pulse: 3975956 332976 0 0 0 0 0
histogram of classification of points:
 4308932 never classified (0)

Nicely visible are the circular scanning patterns of the Velodyne VLP-16 Puck. We also notice that the trajectory of the UAV can be seen in the lasview visualization because the Puck was scanning the drone’s own landing gear. The lasinfo report tells us that point coordinates are stored with too much resolution (mm) and that points do not need to be stored using point type 3 (with RGB colors) because all RGB values are zero. We fix this with an initial run of las2las and also compress the raw strips to the LAZ format on 4 CPUs in parallel.

las2las -i Flaugergues\*.las ^
        -rescale 0.01 0.01 0.01 ^
        -auto_reoffset ^
        -set_point_type 1 ^
        -odir Flaugergues\strips_raw -olaz ^
        -cores 4

Next we do the usual check for flightline alignment with lasoverlap (README) which we consider to be by far the most important quality check. We compare the lowest elevation from different flightline per 25 cm by 25cm cell in all overlap areas. We consider a vertical difference of up to 5 cm as acceptable (color coded as white) and mark differences of over 30 cm (color coded as saturated red or blue).

lasoverlap -i Flaugergues\strips_raw\*.laz -faf ^
           -step 0.25 ^
           -min_diff 0.05 -max_diff 0.3 ^
           -odir Flaugergues\quality -o overlap.png

The vertical difference in open areas between the flightlines is slightly above 5 cm which we consider acceptable in this example. Depending on the application we recommend to investigate further where these differences come from and what consequences they may have for post processing. We also create a color-coded visualization of the last return density per 25 cm by 25 cm cell using lasgrid (README) with blue meaning less than 100 returns per square meter and red meaning more than 4000 returns per square meter.

lasgrid -i Flaugergues\strips_raw\*.laz -merged ^
        -keep_last ^
        -step 0.25 ^
        -point_density ^
        -false -set_min_max 100 4000 ^
        -odir Flaugergues\quality -o density_100_4000.png

Color coded density of last returns per square meter for each 25 cm by 25 cm cell. Blue means 100 or less last returns per square meter. Red means 4000 or more last returns per square meter

As usual we start the LiDAR processing by reorganizing the flightlines into square tiles. Because of the variability in the density that is evident in the visualization above we use lastile (README) to create an adaptive tiling that starts with 200 m by 200 m tiles and then iterate to refine those tiles with over 10 million points down to smaller 25 m by 25 m tiles.

lastile -i Flaugergues\strips_raw\*.laz ^
        -apply_file_source_ID ^
        -tile_size 200 -buffer 8 -flag_as_withheld ^
        -refine_tiling 10000000 ^
        -odir Flaugergues\tiles_raw -o flauge.laz

lastile -i Flaugergues\tiles_raw\flauge*_200.laz ^
        -refine_tiles 10000000 ^
        -olaz ^
        -cores 4

lastile -i Flaugergues\tiles_raw\flauge*_100.laz ^
        -refine_tiles 10000000 ^
        -olaz ^
        -cores 4

lastile -i Flaugergues\tiles_raw\flauge*_50.laz ^
        -refine_tiles 10000000 ^
        -olaz ^
        -cores 4

Subsequent processing is faster when the points have a spatially coherent order. Therefore we rearrange the points into standard space-filling z-order using a call to lassort (README). We run this in parallel on as many cores as it makes sense (i.e. not using more cores than there are physical CPUs).

lassort -i Flaugergues\tiles_raw\flauge*.laz ^
        -odir Flaugergues\tiles_sorted -olaz ^
        -cores 4

Next we classify those points as noise that are isolated on a 3D grid of 1 meter cell size using lasnoise. See the README file of lasnoise for a description on the exact manner in which the isolated points are classified. We do this to eliminate low noise points that would otherwise cause trouble in the subsequent processing.

lasnoise -i Flaugergues\tiles_sorted\flauge*.laz ^
         -step 1 -isolated 5 ^
         -odir Flaugergues\tiles_denoised -olaz ^
         -cores 4

Next we mark the subset of lowest points on a 2D grid of 10 cm cell size with classification code 8 using lasthin (README) while ignoring the noise points with classification code 7 that were marked as noise in the previous step.

lasthin -i Flaugergues\tiles_denoised\flauge*.laz ^
        -ignore_class 7 ^
        -step 0.1 -lowest ^
        -classify_as 8 ^
        -odir Flaugergues\tiles_lowest -olaz ^
        -cores 4

Considering only the resulting points marked with classification 8 we then create a temporary ground classification that we refer to as the “lowest ground”. For this we run lasground (README) with a set of suitable parameters that were found by experimentation on two of the most complex tiles from the center of the survey.

lasground -i Flaugergues\tiles_lowest\flauge*.laz ^
          -ignore_class 0 7 ^
          -step 5 -hyper_fine -bulge 1.5 -spike 0.5 ^
          -odir Flaugergues\tiles_lowest_ground -olaz ^
          -cores 4

We then “thicken” this “lowest ground” by classifying all points that are between 2 cm below and 15 cm above the lowest ground to a temporary classification code 6 using the lasheight (README) tool. Depending on the spread of points in your data set you may want to tighten this range accordingly, for example when processing the flightlines acquired by the Velodyne Puck individually. We picked our range based on the visual experiments with “drop lines” and “rise lines” in the lasview viewer that are shown in images above.

lasheight -i Flaugergues\tiles_lowest_ground\flauge*.laz ^
          -do_not_store_in_user_data ^
          -classify_between -0.02 0.15 6 ^
          -odir Flaugergues\tiles_thick_ground -olaz ^
          -cores 4

The final ground classification is obtained by creating the “median ground” from the “thick ground”. This uses a brand-new option in the lasthin (README) tool of LAStools. The new ‘-percentile 50 10’ option selects the point that is closest to the specified percentile of 50 of all point elevations within a grid cell of a specified size given there are at least 10 points in that cell. The selected point either survives the thinning operation or gets marked with a specified classification code or flag.

lasthin -i Flaugergues\tiles_thick_ground\flauge*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.1 -percentile 50 10 ^
        -classify_as 8 ^
        -odir Flaugergues\tiles_median_ground_10_10cm -olaz ^
        -cores 4

lasthin -i Flaugergues\tiles_median_ground_10_10cm\%NAME%*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.2 -percentile 50 10 ^
        -classify_as 8 ^
        -odir Flaugergues\tiles_median_ground_10_20cm -olaz ^
        -cores 4

lasthin -i Flaugergues\tiles_median_ground_10_20cm\%NAME%*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.4 -percentile 50 10 ^
        -classify_as 8 ^
        -odir Flaugergues\tiles_median_ground_10_40cm -olaz ^
        -cores 4

lasthin -i Flaugergues\tiles_median_ground_10_40cm\flauge*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.8 -percentile 50 10 ^
        -classify_as 8 ^
        -odir Flaugergues\tiles_median_ground_10_80cm -olaz ^
         -cores 4

We now compare a triangulation of the median ground points with a triangulation of the highest and the lowest points per 10 cm by 10 cm cell to demonstrate that – at least in open areas – we really have computed a median ground surface.

Finally we raster the tiles with the las2dem (README) tool onto binary elevation grids in BIL format. Here we make the resolution dependent on the tile size, giving the 25 meter and 50 meter tiles the highest resolution of 10 cm and rasterize the 100 meter and 200 meter tiles at 20 cm and 40 cm respectively.

las2dem -i Flaugergues\tiles_median_ground_10_80cm\*_25.laz ^
        -i Flaugergues\tiles_median_ground_10_80cm\*_50.laz ^
        -keep_class 8 ^
        -step 0.1 -use_tile_bb ^
        -odir Flaugergues\tiles_dtm -obil ^
        -cores 4

las2dem -i Flaugergues\tiles_median_ground_10_80cm\*_100.laz ^
        -keep_class 8 ^
        -step 0.2 -use_tile_bb ^
        -odir Flaugergues\tiles_dtm -obil ^
        -cores 4

las2dem -i Flaugergues\tiles_median_ground_10_80cm\*_200.laz ^
        -keep_class 8 ^
        -step 0.4 -use_tile_bb ^
        -odir Flaugergues\tiles_dtm -obil ^
        -cores 4

Because all LAStools can read BIL files via on the fly conversion from rasters to points we can visually inspect the resulting elevation rasters with the lasview (README) tool. By adding the ‘-faf’ or ‘files_are_flightlines’ argument we treat the BIL files as if they were different flightlines which allows us to assign different color to points from different files to better inspect the transitions between tiles. The ‘-points 10000000’ argument instructs lasview to load up to 10 million points into memory instead of the default 5 million.

lasview -i Flaugergues\tiles_dtm\*.bil -faf ^
        -points 10000000

Final raster tiles in BIL format of three different sizes form seamless DTM.

For visual comparison we also produce a DSM and create hillshades. Note that the workflow for DSM creation shown below produces a “highest DSM” that will always be a few centimeter above the “median DTM”. This will be noticeable only in open areas of the terrain where the DSM and the DTM should coincide and their elevation should be identical.

lasthin -i Flaugergues\tiles_denoised\flauge*.laz ^
        -keep_z_above 110 ^
        -filtered_transform ^
        -set_classification 18 ^
        -ignore_class 7 18 ^
        -step 0.1 -highest ^
        -classify_as 5 ^
        -odir Flaugergues\tiles_highest -olaz ^
        -cores 4

las2dem -i Flaugergues\tiles_highest\*_25.laz ^
        -i Flaugergues\tiles_highest\*_50.laz ^
        -keep_class 5 ^
        -step 0.1 -use_tile_bb ^
        -odir Flaugergues\tiles_dsm -obil ^
        -cores 4

las2dem -i Flaugergues\tiles_highest\*_100.laz ^
        -keep_class 5 ^
        -step 0.2 -use_tile_bb ^
        -odir Flaugergues\tiles_dsm -obil ^
        -cores 4

las2dem -i Flaugergues\tiles_highest\*_200.laz ^
        -keep_class 5 ^
        -step 0.4 -use_tile_bb ^
        -odir Flaugergues\tiles_dsm -obil ^
        -cores 4

We thank YellowScan for challenging us to process their drone LiDAR with LAStools in order to present results at their LiDAR for Drone 2017 Conference and for sharing several example data sets with us, including the one used here.

Scotland’s LiDAR goes Open Data (too)

Following the lead of England and Wales, the Scottish LiDAR is now also open data. The implementation of such an open geospatial policy in the United Kingdom was spear-headed by the Environment Agency of England who started to make all of their LiDAR holdings available as open data. In September 2015 they opened DTM and DSM raster derivatives down to 25 cm resolution and in March 2016 also the raw point clouds went online our compressed and open LAZ format (more info here) – all with the very permissible Open Government Licence v3. This treasure cove of geospatial data was collected by the Environment Agency Geomatics own survey aircraft mainly for flood mapping purposes. The data that had been access restricted for the past 17 years of operation and was made open only after it was shown that restricting access in order to recover costs to finance future operations – a common argument for withholding tax-payer funded data – was nothing but an utter myth. This open data policy has resulted in an incredible re-use of the LiDAR and the Environment Agency has literally been propelled into the role of a “champion for open data” inspiring Wales (possibly the German states of North-Rhine Westfalia and Thuringia) and now also Scotland to open up their geospatial archives as well …

Huge LAS files available for download from the Scottish Open Data portal.

We went to the nice online portal of Scotland to download three files from the Phase II LiDAR for Scotland that are provided as uncompressed LAS files, namely LAS_NN45NE.las, LAS_NN55NE.las, and LAS_NN55NW.las, whose sizes are listed as 1.2 GB, 2.8 GB, and 4.7 GB in the screenshot above. Needless to say that it took quite some time and several restarts (using wget with option ‘-c’) to successfully download these very large LAS files.

laszip -i LAS_NN45NE.las -odix _cm -olaz -rescale 0.01 0.01 0.01 
laszip -i LAS_NN45NE.las -odix _mm -olaz
laszip -i LAS_NN55NE.las -odix _cm -olaz -rescale 0.01 0.01 0.01 
laszip -i LAS_NN55NE.las -odix _mm -olaz
laszip -i LAS_NN55NW.las -odix _cm -olaz -rescale 0.01 0.01 0.01 
laszip -i LAS_NN55NW.las -odix _mm -olaz

After downloading we decided to see how well these files compress with LASzip by running the six commands shown above creating LAZ files when re-scaling of coordinate resolution to centimeter (cm) and LAZ files with the original millimeter (mm) coordinate resolution (i.e. the original scale factors are 0.001 which is somewhat excessive for aerial LiDAR where the error in position per coordinate is typically between 5 cm and 20 cm). Below you see the resulting file sizes for the three different files.

 1,164,141,247 LAS_NN45NE.las
   124,351,690 LAS_NN45NE_cm.laz (1 : 9.4)
   146,651,719 LAS_NN45NE_mm.laz (1 : 7.9)
 2,833,123,863 LAS_NN55NE.las
   396,521,115 LAS_NN55NE_cm.laz (1 : 7.1)
   474,767,495 LAS_NN55NE_mm.laz (1 : 6.0)
 4,664,782,671 LAS_NN55NW.las
   531,454,473 LAS_NN55NW_cm.laz (1 : 8.8)
   629,141,151 LAS_NN55NW_mm.laz (1 : 7.4)

The savings in download time and storage space of storing the LiDAR in LAZ versus LAS are sixfold to tenfold. If I was a tax payer in Scotland and if my government was hosting open data on in the Amazon cloud (i.e. paying for AWS cloud services with my taxes) I would encourage them to store their data in a more compressed format. Some more details on the data.

According to the provided meta data, the Scottish Public Sector LiDAR Phase II dataset was commissioned by the Scottish Government in response to the Flood Risk Management Act (2009). The project was managed by Sniffer and the contract was awarded to Fugro BKS. Airborne LiDAR data was collected for 66 sites (the dataset does not have full national coverage) totaling 3,516 km^2 between 29th November 2012 and 18th April 2014. The point density was a minimum of 1 point/sqm, and approximately 2 points/sqm on average. A DTM and DSM were produced from the point clouds, with 1m spatial resolution. The Coordinate reference system is OSGB 1936 / British National Grid (EPSG code 27700). The data is licensed under an Open Government Licence. However, under the use constraints section it now only states that the following attribution statement must be used to acknowledge the source of the information: “Copyright Scottish Government and SEPA (2014)” but also that Fugro retain the commercial copyright, which is somewhat disconcerting and may require more clarification. According to this tweet a lesser license (NCGL) applies to the raw LiDAR point clouds. Below a lasinfo report for the large LAS_NN55NW.las as well as several visualizations with lasview.

lasinfo (170915) report for LAS_NN55NW.las
reporting all LAS header entries:
 file signature: 'LASF'
 file source ID: 0
 global_encoding: 1
 project ID GUID data 1-4: 00000000-0000-0000-0000-000000000000
 version major.minor: 1.2
 system identifier: 'Riegl LMS-Q'
 generating software: 'Fugro LAS Processor'
 file creation day/year: 343/2016
 header size: 227
 offset to point data: 227
 number var. length records: 0
 point data format: 1
 point data record length: 28
 number of point records: 166599373
 number of points by return: 149685204 14102522 2531075 280572 0
 scale factor x y z: 0.001 0.001 0.001
 offset x y z: 250050 755050 270
 min x y z: 250000.000 755000.000 203.731
 max x y z: 254999.999 759999.999 491.901
reporting minimum and maximum for all LAS point record entries ...
 X -50000 4949999
 Y -50000 4949999
 Z -66269 221901
 intensity 39 2046
 return_number 1 4
 number_of_returns 1 4
 edge_of_flight_line 0 1
 scan_direction_flag 1 1
 classification 1 11
 scan_angle_rank -30 30
 user_data 0 3
 point_source_ID 66 91
 gps_time 38230669.389034 38402435.753789
number of first returns: 149685204
number of intermediate returns: 2813604
number of last returns: 149687616
number of single returns: 135599244
overview over number of returns of given pulse: 135599244 23122229 6754118 1123782 0 0 0
histogram of classification of points:
 287819 unclassified (1)
 109019874 ground (2)
 14476880 low vegetation (3)
 3487218 medium vegetation (4)
 39141518 high vegetation (5)
 165340 building (6)
 13508 rail (10)
 7216 road surface (11)

Kudos to the Scottish government for opening their data. We hereby acknowledge the source of the LiDAR that we have used in the experiments above as “Copyright Scottish Government and SEPA (2014)”.

LAStools Win Big at INTERGEO Taking Home Two Innovation Awards

PRESS RELEASE (for immediate release)
October 2, 2017
rapidlasso GmbH, Gilching, Germany

At INTERGEO 2017 in Berlin, rapidlasso GmbH – the makers of the popular LiDAR processing software LAStools – were awarded top honors in both of the categories they had been nominated for: most innovative software and most innovative startup. The third award for most innovative hardware went to Leica Geosystems for the BLK360 terrestrial scanner. The annual Wichman Innovation Awards have been part of INTERGEO for six years now. Already at the inaugural event in 2012 the open source LiDAR compressor LASzip of rapidlasso GmbH had been nominated, coming in as runner-up in second place.

Dr. Martin Isenburg, the founder and CEO of rapidlasso GmbH, receives the two innovation awards at the ceremony during INTERGEO 2017 in Berlin.

After receiving the two awards Dr. Martin Isenburg, the founder and CEO of rapidlasso GmbH, was quick to thank the “fun, active, and dedicated user community” of the LAStools software for their “incredible support in the online voting”. He pointed out that it was its users who make LAStools more than just an efficient software for processing point clouds. Since 2011, the community surrounding LAStools has constantly grown to several thousand users who help and motivate each other in designing workflows and in solving format issues and processing challenges. They are an integral part of what makes these tools so valuable, so Dr. Isenburg.

About rapidlasso GmbH:
Technology powerhouse rapidlasso GmbH specializes in efficient LiDAR processing tools that are widely known for their high productivity. They combine robust algorithms with efficient I/O and clever memory management to achieve high throughput for data sets containing billions of points. The company’s flagship product – the LAStools software suite – has deep market penetration and is heavily used in industry, government agencies, research labs, and educational institutions. Visit for more information.

LASmoons: Huaibo Mu

Huaibo Mu (recipient of three LASmoons)
Environmental Mapping, Department of Geography
University College London (UCL), UK

This study is a part of the EU-funded Metrology for Earth Observation and Climate project (MetEOC-2). It aims to combine terrestrial and airborne LiDAR data to estimate biomass and allometry for woodland trees in the UK. Airborne LiDAR can capture large amounts of data over large areas, while terrestrial LiDAR can provide much more details of high quality in specific areas. The biomass and allometry for individual specific tree species in 1 ha of Wytham Woods located about 5km north west of the University of Oxford, UK are estimated by combining both airborne and terrestrial LiDAR. Then the bias will be evaluated when estimation are applied on different levels: terrestrial LiDAR level, tree level, and area level. The goal are better insights and a controllable error range in the bias of biomass and allometry estimates for woodland trees based on airborne LiDAR.

The study aims to find the suitable parameters of allometric equations for different specific species and establish the relationship between the tree height and stem diameter and crown diameter to be able to estimate the above ground biomass using airborne LiDAR. The biomass estimates under different levels are then compared to evaluate the bias and for the total 6ha of Wytham Woods for calibration and validation. Finally the results are to be applied to other woodlands in the UK. The LiDAR processing tasks for which LAStools are used mainly center around the creation of suitable a Canopy Height Model (CHM) from the airborne LiDAR.

+ Airborne LiDAR data produced by Professor David Coomes (University of Cambridge) with Airborne Research and Survey Facility (ARSF) Project code of RG13_08 in June 2014. The average point density is about 5.886 per m^2.
+ Terrestrial LiDAR data collected by UCL’s team leader by Dr. Mat Disney and Dr. Kim Calders in order to develop very detailed 3D models of the trees.
+ Fieldwork from the project “Initial Results from Establishment of a Long-term Broadleaf Monitoring Plot at Wytham Woods, Oxford, UK” by Butt et al. (2009).

LAStools processing:
1) check LiDAR quality as described in these videos and articles [lasinfo, lasvalidate, lasoverlap, lasgrid, las2dem]
2) classify into ground and non-ground points using tile-based processing  [lastile, lasground]
3) generate a Digital Terrain Model (DTM) [las2dem]
4) compute height of points and delete points higher than maximum tree height obtained from terrestrial LiDAR [lasheight]
5) convert points into disks with 10 cm diameter to conservatively account for laser beam width [lasthin]
6) generate spike-free Digital Surface Model (DSM) based on algorithm by Khosravipour et al. (2016) [las2dem]
7) create Canopy Height Model (CHM) by subtracting DTM from spike-free DSM [lasheight].

Butt, N., Campbell, G., Malhi, Y., Morecroft, M., Fenn, K., & Thomas, M. (2009). Initial results from establishment of a long-term broadleaf monitoring plot at Wytham Woods, Oxford, UK. University Oxford, Oxford, UK, Rep.
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.
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).

Removing Excessive Low Noise from Dense-Matching Point Clouds

Point clouds produced with dense-matching by photogrammetry software such as SURE, Pix4D, or Photoscan can include a fair amount of the kind of “low noise” as seen below. Low noise causes trouble when attempting to construct a Digital Terrain Model (DTM) from the points as common algorithm for classifying points into ground and non-ground points – such as lasground – tend to “latch onto” those low points, thereby producing a poor representation of the terrain. This blog post describes one possible LAStools workflow for eliminating excessive low noise. It was developed after a question in the LAStools user forum by LASmoons holder Muriel Lavy who was able to share her noisy data with us. See this, this, this, thisthis, and this blog post for further reading on this topic.

Here you can download the dense matching point cloud that we are using in the following work flow:

We leave the usual inspection of the content with lasinfolasview, and lasvalidate that we always recommend on newly obtained data as an exercise to the reader. Note that a check for proper alignment of flightlines with lasoverlap that we consider mandatory for LiDAR data is not applicable for dense-matching points.

With lastile we turn the original file with 87,261,083 points into many smaller 500 by 500 meter tiles for efficient multi-core processing. Each tile is given a 25 meter buffer to avoid edge artifacts. The buffer points are marked as withheld for easier on-the-fly removal. We add a (terser) description of the WGS84 UTM zone 32N to each tile via the corresponding EPSG code 32632:
lastile -i muriel\20161127_Pancalieri_UTM.laz ^
        -tile_size 500 -buffer 25 -flag_as_withheld ^
        -epsg 32632 ^
        -odir muriel\tiles_raw -o panca.laz
Because dense-matching points often have a poor point order in the files they get delivered in we use lassort to rearrange them into a space-filling curve order as this will speed up most following processing steps:
lassort -i muriel\tiles_raw\panca*.laz ^
        -odir muriel\tiles_sorted -olaz ^
        -cores 7
We then run lasthin to reclassify the highest point of every 2.5 by 2.5 meter grid cell with classification code 8. As the spacing of the dense-matched points is around 40 cm in both x and y, around 40 points will fall into each such grid cell from which the highest is then classified as 8:
lasthin -i muriel\tiles_sorted\panca*.laz ^
        -step 2.5 ^
        -highest -classify_as 8 ^
        -odir muriel\tiles_thinned -olaz ^
        -cores 7
Considering only those points classified as 8 in the last step we then run lasnoise to find points that are highly isolated in wide and flat neighborhoods that are then reclassified as 7. See the README file of lasnoise for a detailed explanation of the different parameters:
lasnoise -i muriel\tiles_thinned\panca*.laz ^
         -ignore_class 0 ^
         -step_xy 5 -step_z 0.1 -isolated 4 ^
         -classify_as 7 ^
         -odir muriel\tiles_isolated -olaz ^
         -cores 7
Now we run a temporary ground classification of only (!!!) on those points that are still classified as 8 using the default parameters of lasground. Hence we only use the points that were the highest points on the 2.5 by 2.5 meter grid and that were not classified as noise in the previous step. See the README file of lasground for a detailed explanation of the different parameters:
lasground -i muriel\tiles_isolated\panca*.laz ^
          -city -ultra_fine -ignore_class 0 7 ^
          -odir muriel\tiles_temp_ground -olaz ^
          -cores 7
The result of this temporary ground filtering is then merely used to mark all points that are 0.5 meter below the triangulated TIN of these temporary ground points with classification code 12 using lasheight. See the README file of lasheight for a detailed explanation of the different parameters:
lasheight -i muriel\tiles_temp_ground\panca*.laz ^
          -do_not_store_in_user_data ^
          -classify_below -0.5 12 ^
          -odir muriel\tiles_temp_denoised -olaz ^
          -cores 7
In the resulting tiles the low noise (but also many points above the ground) are now marked and in a final step we produce properly classified denoised tiles by re-mapping the temporary classification codes to conventions that are more consistent with the ASPRS LAS specification using las2las:
las2las -i muriel\tiles_temp_denoised\panca*.laz ^
        -change_classification_from_to 1 0 ^
        -change_classification_from_to 2 0 ^
        -change_classification_from_to 7 0 ^
        -change_classification_from_to 12 7 ^
        -odir muriel\tiles_denoised -olaz ^
        -cores 7
Let us visually check what each of the above steps has produced by zooming in on a 300 meter by 100 meter strip of points with the bounding box (388500,4963125) to (388800,4963225) in tile ‘panca_388500_4963000.laz’:
The final classification of all points that are not already classified as noise (7) into ground (2) or non-ground (1) was done with a final run of lasground. See the README file of lasground for a detailed explanation of the different parameters:
lasground -i muriel\tiles_denoised\panca*.laz ^
          -ignore_class 7 ^
          -city -ultra_fine ^
          -odir muriel\tiles_ground -olaz ^
          -cores 7
Then we create a seamless hill-shaded DTM tiles by triangulating all the points classified as ground into a temporary TIN (including those in the 25 meter buffer) and then rasterizing only the inner 500 meter by 500 meter of each tile with option ‘-use_tile_bb’ of las2dem. For more details on the importance of buffers in tile-based processing see this blog post here.
las2dem -i muriel\tiles_ground\panca*.laz ^
        -keep_class 2 ^
        -step 1 -hillshade ^
        -use_tile_bb ^
        -odir muriel\tiles_dtm -opng ^
        -cores 7

And here the original DSM side-by-side with resulting DTM after low noise removal. One dense forested area near the center of the data was not entirely removed due to the lack of ground points in this area. Integrating external ground points or manual editing with lasview are two possible way to rectify these few remaining errors …