Plots to Stands: Producing LiDAR Vegetation Metrics for Imputation Calculations

Some professionals in remote sensing find LAStools a useful tool to extract statistical metrics from LiDAR that are used to make estimations about a larger area of land from a small set of sample plots. Common applications are prediction of the timber volume or the above-ground biomass for entire forests based on a number of representative plots where exact measurements were obtained with field work. The same technique can also be used to make estimations about animal habitat or coconut yield or to classify the type of vegetation that covers the land. In this tutorial we describe the typical workflow for computing common metrics for smaller plots and larger areas using LAStools.

Download these six LiDAR tiles (1, 2, 3, 4, 5, 6) from a Eucalyptus plantation in Brazil to follow along the step by step instructions of this tutorial. This data is courtesy of Suzano Pulp and Paper. Please also download the two shapefiles that delineate the plots where field measurements were taken and the stands for which predictions are to be made. You should download version 170327 (or higher) of LAStools due to some recent bug fixes.

Quality Checking

Before processing newly received LiDAR data we always perform a quality check first. This ranges from visual inspection with lasview, to printing textual content reports and attribute histograms with lasinfo, to flight-line alignment checks with lasoverlap, pulse density and pulse spacing checks with lasgrid and las2dem, and completeness-of-returns check with lassort followed by lasreturn.

lasinfo -i tiles_raw\CODL0003-C0006.laz ^
        -odir quality -odix _info -otxt

The lasinfo report tells us that there is no projection information. However, we remember that this Brazilian data was in the common SIRGAS 2000 projection and try for a few likely UTM zones whether the hillshaded DSM produced by las2dem falls onto the right spot in Google Earth.

las2dem -i tiles_raw\CODL0003-C0006.laz ^
        -keep_first -thin_with_grid 1 ^
        -hillshade -epsg 31983 ^
        -o epsg_check.png

Hillshaded DSM and Google Earth imagery align for EPSG code 31983

The lasinfo report also tells us that the xyz coordinates are stored with millimeter resolution which is a bit of an overkill. For higher and faster LASzip compression we will later lower this to a more appropriate centimeter resolution. It further tells us that the returns are stored using point type 0 and that is a bit unfortunate. This (older) point type does not have a GPS time stamp so that some quality checks (e.g. “completeness of returns” with lasreturn) and operations (e.g. “resorting of returns into acquisition order” with lassort) will not be possible. Fortunately the min-max range of the ‘point source ID’ suggests that this point attribute is correctly populated with flightline numbers so that we can do a check for overlap and alignment of the different flightlines that contribute to the LiDAR in each tile.

lasoverlap -i tiles_raw\*.laz ^
           -min_diff 0.2 -max_diff 0.4 ^
           -epsg 31983 ^
           -odir quality -opng ^
           -cores 3

We run lasoverlap to visualize the amount of overlap between flightlines and the vertical differences between them. The produced images (see below) color code the number of flightlines and the maximum vertical difference between any two flightlines as seen below. Most of the area is cyan (2 flightlines) except in the bottom left where the pilot was sloppy and left some gaps in the yellow seams (3 flightlines) so that some spots are only blue (1 flightline). We also see that two tiles in the upper left are partly covered by a diagonal flightline. We will drop that flightline later to create a more uniform density.across the tiles. The mostly blue areas in the difference are mostly aligned with features in the landscape and less with the flightline pattern. Closer inspection shows that these vertical difference occur mainly in the dense old growth forests with species of different heights that are much harder to penetrate by the laser than the uniform and short-lived Eucalyptus plantation that is more of a “dead forest” with little undergrowth or animal habitat.

Interesting observation: The vertical difference of the lowest return from different flightlines computed per 2 meter by 2 meter grid cell could maybe be used a new forestry metric to help distinguish mono cultures from natural forests.

lasgrid -i tiles_raw\*.laz ^
        -keep_last ^
        -step 2 -point_density ^
        -false -set_min_max 10 20 ^
        -epsg 31983 ^
        -odir quality -odix _d_2m_10_20 -opng ^
        -cores 3

lasgrid -i tiles_raw\*.laz ^
        -keep_last ^
        -step 5 -point_density ^
        -false -set_min_max 10 20 ^
        -epsg 31983 ^
        -odir quality -odix _d_5m_10_20 -opng ^
        -cores 3

We run lasgrid to visualize the pulse density per 2 by 2 meter cell and per 5 by 5 meter cell. The produced images (see below) color code the number of last return per square meter. The impact of the tall Eucalyptus trees on the density per cell computation is evident for the smaller 2 meter cell size in form of a noisy blue/red diagonal in the top right as well as a noisy blue/red area in the bottom left. Both of those turn to a more consistent yellow for the density per cell computation with larger 5 meter cells. Immediately evident is the higher density (red) for those parts or the two tiles in the upper left that are covered by the additional diagonal flightline. The blueish area left to the center of the image suggests a consistently lower pulse density whose cause remains to be investigated: Lower reflectivity? Missing last returns? Cloud cover?

The lasinfo report suggests that the tiles are already classified. We could either use the ground classification provided by the vendor or re-classify the data ourselves (using lastilelasnoise, and lasground). We check the quality of the ground classification by visually inspecting a hillshaded DTM created with las2dem from the ground returns. We buffer the tiles on-the-fly for a seamless hillshade without artifacts along tile boundaries by adding ‘-buffered 25’ and ‘-use_orig_bb’ to the command-line. To speed up reading the 25 meter buffers from neighboring tiles we first create a spatial indexing with lasindex.

lasindex -i tiles_raw\*.laz ^
         -cores 3

las2dem -i tiles_raw\*.laz ^
        -buffered 25 ^
        -keep_class 2 -thin_with_grid 0.5 ^
        -use_orig_bb ^
        -hillshade -epsg 31983 ^
        -odir quality -odix _dtm -opng ^
        -cores 3

hillshaded DTM tiles generated with las2dem and on-the-fly buffering

The resulting hillshaded DTM shows a few minor issues in the ground classification but also a big bump (above the mouse cursor). Closer inspection of this area (you can cut it from the larger tile using the command-line below) shows that there is a group of miss-classified points about 1200 meters below the terrain. Hence, we will start from scratch to prepare the data for the extraction of forestry metrics.

las2las -i tiles_raw\CODL0004-C0006.laz ^
        -inside_tile 207900 7358350 100 ^
        -o bump.laz

lasview -i bump.laz

bump in hillshaded DTM caused by miss-classified low noise

Data Preparation

Using lastile we first tile the data into smaller 500 meter by 500 meter tiles with 25 meter buffer while flagging all points in the buffer as ‘withheld’. In the same step we lower the resolution to centimeter and put nicer a coordinate offset in the LAS header. We also remove the existing classification and classify all points that are much lower than the target terrain as class 7 (aka noise). We also add CRS information and give each tile the base name ‘suzana.laz’.

lastile -i tiles_raw\*.laz ^
        -rescale 0.01 0.01 0.01 ^
        -auto_reoffset ^
        -set_classification 0 ^
        -classify_z_below_as 500.0 7 ^
        -tile_size 500 ^
        -buffer 25 -flag_as_withheld ^
        -epsg 31983 ^
        -odir tiles_buffered -o suzana.laz

With lasnoise we mark the many isolated points that are high above or below the terrain as class 7 (aka noise). Using two tiles we played around with the ‘step’ parameters until we find good parameter settings. See the README of lasnoise for the exact meaning and the choice of parameters for noise classification. We look at one of the resulting tiles with lasview.

lasnoise -i tiles_buffered\*.laz ^
         -step_xy 4 -step_z 2 ^
         -odir tiles_denoised -olaz ^
         -cores 3

lasview -i tiles_denoised\suzana_206000_7357000.laz ^
        -color_by_classification ^
        -win 1024 192

noise points shown in pink: all points (top), only noise points (bottom)

Next we use lasground to classify the last returns into ground (2) and non-ground (1). It is important to ignore the noise points with classification 7 to avoid the kind of bump we saw in the vendor-delivered classification. We again check the quality of the computed ground classification by producing a hillshaded DTM with las2dem. Here the las2dem command-line is sightly different as instead of buffering on-the-fly we use the buffers stored with each tile.

lasground -i tiles_denoised\*.laz ^
          -ignore_class 7 ^
          -nature -extra_fine ^
          -odir tiles_ground -olaz ^
          -cores 3

las2dem -i tiles_ground\*.laz ^
        -keep_class 2 -thin_with_grid 0.5 ^
        -hillshade ^
        -use_tile_bb ^
        -odir quality -odix _dtm_new -opng ^
        -cores 3

Finally, with lasheight we compute how high each return is above the triangulated surface of all ground returns and store this height value in place of the elevation value into the z coordinate using the ‘-replace_z’ switch. This height-normalizes the LiDAR in the sense that all ground returns are set to an elevation of 0 while all other returns get an elevation relative to the ground. The result are height-normalized LiDAR tiles that are ready for producing the desired forestry metrics.

lasheight -i tiles_ground\*.laz ^
          -replace_z ^
          -odir tiles_normalized -olaz ^
          -cores 3
Metric Production

The tool for computing the metrics for the entire area as well as for the individual field plots is lascanopy. Which metrics are well suited for your particular imputation calculation is your job to determine. Maybe first compute a large number of them and then eliminate the redundant ones. Do not use any point from the tile buffers for these calculations. We had flagged them as ‘withheld’ during the lastile operation, so they are easy to drop. We also want to drop the noise points that were classified as 7. And we were planning to drop that additional diagonal flightline we noticed during quality checking. We generated two lasinfo reports with the ‘-histo point_source 1’ option for the two tiles it was covering. From the two histograms it was easy to see that the diagonal flightline has the number 31.

First we run lascanopy on the 11 plots that you can download here. When running on plots it makes sense to first create a spatial indexing with lasindex for faster querying so that only those tiny parts of the LAZ file need to be loaded that actually cover the plots.

lasindex -i tiles_normalized\*.laz ^
         -cores 3

lascanopy -i tiles_normalized\*.laz -merged ^
          -drop_withheld ^
          -drop_class 7 ^
          -drop_point_source 31 ^
          -lop WKS_PLOTS.shp ^
          -cover_cutoff 3.0 ^
          -cov -dns ^
          -height_cutoff 2.0 ^
          -c 2.0 999.0 ^
          -max -avg -std -kur ^
          -p 25 50 75 95 ^
          -b 30 50 80 ^
          -d 2.0 5.0 10.0 50.0 ^
          -o plots.csv

The resulting ‘plots.csv’ file you can easily process in other software packages. It contains one line for each polygonal plot listed in the shapefile that lists its bounding box followed by all the requested metrics. But is why is there a zero maximum height (marked in orange) for plots 6 though 10? All height metrics are computed solely from returns that are higher than the ‘height_cutoff’ that was set to 2 meters. We added the ‘-c 2.0 999.0’ absolute count metric to keep track of the number of returns used in these calculations. Apparently in plots 6 though 10 there was not a single return above 2 meters as the count (also marked in orange) is zero for all these plots. Turns out this Eucalyptus stand had recently been harvested and the new seedlings are still shorter than 2 meters.

more plots.csv
index,min_x,min_y,max_x,max_y,max,avg,std,kur,p25,p50,p75,p95,b30,b50,b80,c00,d00,d01,d02,cov,dns
0,206260.500,7358289.909,206283.068,7358312.477,11.23,6.22,1.91,2.26,4.71,6.01,7.67,9.5,26.3,59.7,94.2,5359,18.9,41.3,1.5,76.3,60.0
1,206422.500,7357972.909,206445.068,7357995.477,13.54,7.5,2.54,1.97,5.32,7.34,9.65,11.62,26.9,54.6,92.2,7030,12.3,36.6,13.3,77.0,61.0
2,206579.501,7358125.909,206602.068,7358148.477,12.22,5.72,2.15,2.5,4.11,5.32,7.26,9.76,46.0,73.7,97.4,4901,24.8,28.7,2.0,66.8,51.2
3,206578.500,7358452.910,206601.068,7358475.477,12.21,5.68,2.23,2.64,4.01,5.14,7.18,10.04,48.3,74.1,95.5,4861,25.7,26.2,2.9,68.0,50.2
4,206734.501,7358604.910,206757.068,7358627.478,15.98,10.3,2.18,2.64,8.85,10.46,11.9,13.65,3.3,27.0,91.0,4946,0.6,32.5,44.5,91.0,77.5
5,207043.501,7358761.910,207066.068,7358784.478,15.76,10.78,2.32,3.43,9.27,11.03,12.49,14.11,3.2,20.7,83.3,4819,1.5,24.7,51.0,91.1,76.8
6,207677.192,7359630.526,207699.760,7359653.094,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.0,0.0,0.0,0,0.0,0.0,0.0,0.0,0.0
7,207519.291,7359372.366,207541.859,7359394.934,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.0,0.0,0.0,0,0.0,0.0,0.0,0.0,0.0
8,207786.742,7359255.850,207809.309,7359278.417,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.0,0.0,0.0,0,0.0,0.0,0.0,0.0,0.0
9,208159.017,7358997.344,208181.584,7359019.911,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.0,0.0,0.0,0,0.0,0.0,0.0,0.0,0.0
10,208370.909,7358602.565,208393.477,7358625.133,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.0,0.0,0.0,0,0.0,0.0,0.0,0.0,0.0

Then we run lascanopy on the entire area and produce one raster per tile for each metric. Here we remove the buffered points with the ‘-use_tile_bb’ switch that also ensures that the produced rasters have exactly the extend of the tiles without buffers. Again, it is imperative that you download the version 170327 (or higher) of LAStools for this to work correctly.

lascanopy -version
LAStools (by martin@rapidlasso.com) version 170327 (academic)

lascanopy -i tiles_normalized\*.laz ^
          -use_tile_bb ^
          -drop_class 7 ^
          -drop_point_source 31 ^
          -step 10 ^
          -cover_cutoff 3.0 ^
          -cov -dns ^
          -height_cutoff 2.0 ^
          -c 2.0 999.0 ^
          -max -avg -std -kur ^
          -p 25 50 75 95 ^
          -b 30 50 80 ^
          -d 2.0 5.0 10.0 50.0 ^
          -odir tile_metrics -oasc ^
          -cores 3

The resulting rasters in ASC format can easily be previewed using lasview for some “sanity checking” that our metrics make sense and to get a quick overview about what these metrics look like.

lasview -i tile_metrics\suzana_*max.asc
lasview -i tile_metrics\suzana_*p95.asc
lasview -i tile_metrics\suzana_*p50.asc
lasview -i tile_metrics\suzana_*p25.asc
lasview -i tile_metrics\suzana_*cov.asc
lasview -i tile_metrics\suzana_*d00.asc
lasview -i tile_metrics\suzana_*d01.asc
lasview -i tile_metrics\suzana_*b30.asc
lasview -i tile_metrics\suzana_*b80.asc

The maximum height rasters are useful to inspect more closely as they will immediately tell us if there was any high noise point that slipped through the cracks. And indeed it happened as we see a maximum of 388.55 meters for of the 10 by 10 meter cells. As we know the expected height of the trees we could have added a ‘-drop_z_above 70’ to the lascanopy command line. Careful, however, when computing forestry metrics in strongly sloped terrains as the terrain slope can significantly lift up returns to heights much higher than that of the tree. This is guaranteed to happen for LiDAR returns from branches that are extending horizontally far over the down-sloped part of the terrain as shown in this paper here.

We did not use the shapefile for the stands in this exercise. We could have clipped the normalized LiDAR points to these stands using lasclip as shown in the GUI below before generating the raster metrics. This would have saved space and computation time as many of the LiDAR points lie outside of the stands. However, it might be better to do that clipping step on the rasters in whichever GIS software or statistics package you are using for the imputation computation to properly account for partly covered raster cells along the stand boundary. This could be subject of another blog article … (-:

not all LiDAR was needed to compute metrics for

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 ^
        -stdout
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 ^
        -stdout
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 ^
        -stdout
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 ^
        -stdout
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 ^
        -color_by_return

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

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

LASmoons: Chloe Brown

Chloe Brown (recipient of three LASmoons)
Geosciences, School of Geography
University of Nottingham, UK

Background:
Malaysia’s North Selangor peat swamp forest is experiencing rapid and large scale conversion of peat swampland to oil palm agriculture, contrary to prevailing environmental guidelines. Given the global importance of tropical peat lands, and the uncertainties surrounding historical and future oil palm development, quantifying the spatial distribution of ecosystem service values, such as climate mitigation, is key to understanding the trade-offs associated with anthropogenic land use change.
The study explores the capabilities and methods of remote sensing and field-based data sets for extracting relevant metrics for the assessment of carbon stocks held in North Selangor peat swamp forest reserve, estimating both the current carbon stored in the above and below ground biomass, as well as the changes in carbon stock over time driven by anthropogenic land use change. Project findings will feed directly into peat land management practices and environmental accounting in Malaysia through the Tropical Catchments Research Initiative (TROCARI), and support the Integrated Management Plan of the Selangor State Forest Department (see here for a sample).

some clever caption

Goal:
LiDAR data is now seen as the practical option when assessing canopy height over large scales (Fassnacht et al., 2014), with Lucas et al., (2008) believing LiDAR data to produce more accurate tree height estimates than those derived from manual field based methods. At this stage of the project, the goal is to produce a high quality LiDAR-derived Canopy Height Model (CHM) following the “pit-free” algorithm of Khosravipour et.al., 2014 using the LAStools software.

Data:
+ LiDAR provided by the Natural Environment Research Council (NERC) Airborne Research and Survey Facility’s 2014 Malaysia Campaign.
+ covers 685 square kilometers (closed source)
+ collected with Leica ALS50-II LiDAR system
+ average pulse spacing < 1 meter, average pulse density 1.8 per square meter

LAStools processing:
1) Create 1000 meter tiles with 35 meter buffer to avoid edge artifacts [lastile]
2) Remove noise points (class 7) that are already classified [las2las]
3) Classify point clouds into ground (class 2) and non-ground (class 1) [lasground]
4) Generate normalized above-ground heights [lasheight]
5) Create DSM and DTM [las2dem]
6) Generate a pit-free Canopy Height Model (CHM) as described here [lasthin, las2dem, lasgrid]
7) Generate a spike-free Canopy Height Model (CHM) as described here for comparison [las2dem]

References:
Fassnacht, F.E., Hartig, F., Latifi, H., Berger, C., Hernández, J., Corvalán, and P., Koch, B. (2014). Importance of sample size, data type and prediction method for remote sensing-based estimations of above-ground forest biomass. Remote Sensing. Environment. 154, 102–114.
Khosravipour, A., Skidmore, A. K., Isenburg, M., Wang, T., and Hussin, Y. A. (2014). Generating pit-free canopy height models from airborne LiDAR. Photogrammetric Engineering & Remote Sensing, 80(9), 863-872.
Lucas, R. M., Lee, A. C., and Bunting, P. J., (2008). Retrieving forest biomass through integration of casi and lidar data. International Journal of Remote Sensing, 29 (5), 1553-1577.

LASmoons: Elia Palop-Navarro

Elia Palop-Navarro (recipient of three LASmoons)
Research Unit in Biodiversity (UO-PA-CSIC)
University of Oviedo, SPAIN.

Background:
Old-growth forests play an important role in biodiversity conservation. However, long history of human transformation of the landscape has led to the existence of few such forests nowadays. Its structure, characterized by multiple tree species and ages, old trees and abundant deadwood, is particularly sensible to management practices (Paillet et al. 2015) and requires long time to recover from disturbance (Burrascano et al. 2013). Within protected areas we would expect higher proportions of old-growth forests since these areas are in principle managed to ensure conservation of natural ecosystems and processes. Nevertheless, most protected areas in the EU sustained use and exploitation in the past, or even still do.

lasmoons_elia_palopnavarro_0

Part of the study area. Dotted area corresponds to forest surface under protection.

Goal:
Through the application of a model developed in the study area, using public LiDAR and forest inventory data (Palop-Navarro et al. 2016), we’d like to know how much of the forest in a network of mountain protected areas retains structural attributes compatible with old-growth forests. The LiDAR processing tasks which LAStools will be used for involve a total of 614,808 plots in which we have to derive height metrics, such as mean or median canopy height and its variability.

Vegetation profile colored by height in a LiDAR sample of the study area.

Vegetation profile colored by height in a LiDAR sample of the study area.

Data:
+ Public LiDAR data that can be downloaded here with mean pulse density 0.5 points per square meter. This data has up to 5 returns and is already classified into ground, low, mid or high vegetation, building, noise or overlapped.
+ The area covers forested areas within protected areas in Cantabrian Mountains, occupying 1,207 km2.

LAStools processing:
1) quality checking of the data as described in several videos and blog posts [lasinfo, lasvalidate, lasoverlap, lasgrid, las2dem]
2) use existing ground classification (if quality suffices) to normalize the elevations of to heights above ground using tile-based processing with on-the-fly buffers of 50 meters to avoid edge artifacts [lasheight]
3) compute height-based forestry metrics (e.g. ‘-avg’, ‘-std’, and ‘-p 50’) for each plot in the study area [lascanopy]

References:
Burrascano, S., Keeton, W.S., Sabatini, F.M., Blasi, C. 2013. Commonality and variability in the structural attributes of moist temperate old-growth forests: a global review. Forest Ecology and Management 291:458-479.
Paillet, Y., Pernot, C., Boulanger, V., Debaive, N., Fuhr, M., Gilg, O., Gosselin, F. 2015. Quantifying the recovery of old-growth attributes in forest reserves: A first reference for France. Forest Ecology and Management 346:51-64.
Palop-Navarro, E., Bañuelos, M.J., Quevedo, M. 2016. Combinando datos lidar e inventario forestal para identificar estados avanzados de desarrollo en bosques caducifolios. Ecosistemas 25(3):35-42.

Prototype for “native LAS 1.4 extension” of LASzip LiDAR Compressor Released

PRESS RELEASE (for immediate release)
February 13, 2017
rapidlasso GmbH, Gilching, Germany

Just in time for ILMF 2017 in Denver, the makers of the popular LiDAR processing software LAStools announce that the prototype for the “native LAS 1.4 extension” of their award-winning open source LASzip LiDAR compressor is ready for testing. An update to the compressed LAZ format had become necessary due to a core change in the ASPRS LAS 1.4 specification which had introduced several new point types.

A new feature of the updated LASzip compressor is the ability to selectively decompress of only those attributes of each point that really are needed by the application that is reading the LAZ file. Minimally this will be the x and y coordinate of each point and the return counts, which are sufficient to – for example – calculate the exact extend of the survey area. Most applications will also want to access z coordinate. However, the intensities, the GPS times, the RGB or NIR colors, and the new “Extra Bytes” are often not needed. As the updated LAZ format compresses these different attributes into separate layers, their decompression can then be skipped. Therefore sometimes only 40% of a compressed LAZ file needs to be decompressed to access the coordinates of points with many attributes.

percentage of bytes in a compressing LAZ file corresponding to different point attributes

The percentages of a compressed LAZ file used to encode different point attributes for two example LAS 1.4 files.

The new LASzip prototype is currently being crowd-tested. Interested parties who already have holdings of LAS 1.4 files with point types 6 to 10 may send an email to ‘lasproto@rapidlasso.com’ to participate in these tests.

The release of the new LASzip compressor comes more than a year late as development had been intentionally delayed to give ESRI an opportunity to contribute their needs and ideas to create a joint open format with the community and avoid LiDAR format fragmentation. Sadly, this effort ultimately failed.

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 http://rapidlasso.com for more information.

LASmoons: Jesús García Sánchez

Jesús García Sánchez (recipient of three LASmoons)
Landscapes of Early Roman Colonization (LERC) project
Faculty of Archaeology, Leiden University, The Netherlands

Background:
Our project Landscapes of Early Roman Colonization (LERC) has been studying the hinterland of the Latin colony of Aesernia (Molise region, Italy) using several non-destructive techniques, chiefly artefactual survey, geophysics, and interpretation of aerial photographs. Nevertheless large areas of the territory are covered by the dense forests of the Matese mountains, a ridge belonging the Apennine chain, or covered by bushes due to the abandonment of the countryside. The project won’t be complete without integrating the marginal, remote and forested areas into our study of the Roman hinterland. Besides, it’s also relevant to discuss the feasibility of LiDAR data sets in the study of Mediterranean landscapes and its role within contemporary Landscape Archaeology.

some clever caption

LiDAR coverage in Molise region, Italy.

Goal:
+ to study in detail forested areas in the colonial hinterland of Aesernia.
+ to found the correct parameters of the classification algorithm to be able to locate possible archaeological structures or to document appropriately those we already known.
+ to document and create new visualization of hill-top fortified sites that belong to the indigenous population and are currently poorly studied due to inaccessibility and forest coverage (Monte San Paolo, Civitalla, Castelriporso, etc.)
+ to demonstrate the archaeological potential of LiDAR data in Italy and help other scholars to work with that kind of data, explaining basic information about data quality, where and how to acquire imagery and examples of application in archaeology. A paper entitled “Working with ALS – LiDAR data in Central South Italy. Tips and experiences”, will be presented in the International Mediterranean Survey Workshop by the end of February in Athens.

Civitella hillfort (Longano, IS) and its local context: ridges and forest belonging to the Materse mountains and the Appenines.

Data:
Recently the LERC project has acquired a large LiDAR dataset created by the Italian Geoportale Nazionale and the Minisstero dell’Ambiente e della Tutella del Territorio e del Mare. The data was produced originally to monitor land-slides and erosive risk.
The average point resolution is 1 meter.
+ The data sets were cropped originally in 1 sq km. tiles by the Geoportale Nazionale for distribution purposes.

LAStools processing:
1) data is provided in *.txt files thus the first step is to create appropriate LAS files to work with [txt2las]
2) combine areas of circa 16 sq km (still fewer than 20 million points to be processed in one piece with LAStools) in the surroundings of the colony of Aesernia and in the Matese mountains [lasmerge]
3) assign the correct projection to the data [lasmerge or las2las]
4) extract the care-earth with the best-fitting parameters [lasground or lasground_new]
5) create bare-earth terrain rasters as a first step to visualize and analyze the area [lasdem]