Integrating External Ground Points in Forests to Improve DTM from Dense-Matching Photogrammetry

The biggest problem of generating a Digital Terrain Model (DTM) from the photogrammetric point clouds that are produced from aerial imagery with dense-matching software such as SURE, Pix4D, or Photoscan is dense vegetation: when plants completely cover the terrain not a single point is generated on the ground. This is different for LiDAR point clouds as the laser can even penetrate dense multi-level tropical forests. The complete lack of ground points in larger vegetated areas such as closed forests or dense plantations means that the many processing workflows for vegetation analysis that have been developed for LiDAR cannot be used for photogrammetric point clouds  … unless … well unless we are getting those missing ground points some other way. In the following we see how to integrate external ground points to generate a reasonable DTM under a dense forest with LAStools. See this, this, this, this, and this article for further reading.

Here you can download the dense matching point cloud, the manually collected ground points, and the forest stand delineating polygon that we are using in the following example work flow:

We leave the usual inspection of the content with lasinfo and lasview that we always recommend on newly obtained data as an exercise to the reader. Using las2dem and lasgrid we created the Google Earth overlays shown above to visualize the extent of the dense matched point cloud and the distribution of the manually collected ground points:

las2dem -i DenseMatching.laz ^
        -thin_with_grid 1.0 ^
        -extra_pass ^
        -step 2.0 ^
        -hillshade ^
        -odix _hill_2m -opng

lasgrid -i ManualGround.laz ^
        -set_RGB 255 0 0 ^
        -step 10 -rgb ^
        -odix _grid_10m -opng

Attempts to ground-classify the dense matching point cloud directly are futile as there are no ground points under the canopy in the heavily forested area. Therefore 558 ground points were manually surveyed in the forest of interest that are around 50 to 120 meters apart from another. We show how to integrate these points into the dense matching point cloud such that we can successfully extract bare-earth information from the data.

In the first step we “densify” the manually collected ground points by interpolating them with triangles onto a raster of 2 meter resolution that we store as LAZ points with las2dem. You could consider other interpolation schemes to “densify” the ground points, here we use simple linear interpolation to prove the concept. Due to the varying distance between the manually surveyed ground points we allow interpolating triangles with edge lengths of up to 125 meters. These triangles then also cover narrow open areas next to the forest, so we clip the interpolated ground points against the forest stand delineating polygon with lasclip to classify those points that are really in the forest as “key points” (class 8) and all others as “noise” (class 7).

las2dem -i ManualGround.laz ^
        -step 2 ^
        -kill 125 ^
        -odix _2m -olaz

lasclip -i ManualGround_2m.laz ^
        -set_classification 7 ^ 
        -poly forest.shp ^
        -classify_as 8 -interior ^
        -odix _forest -olaz

Below we show the resulting densified ground points colored by elevation that survive the clipping against the forest stand delineating polygon and were classified as “key points” (class 8). The interpolated ground points in narrow open areas next to the forest that fall outside this polygon were classified as “noise” (class 7) and are shown in violet. They will be dropped in the next step.

We then merge the dense matching points with the densified manual ground points (while dropping all the violet points marked as noise) as input to lasthin and reclassify the lowest point per 1 meter by 1 meter with a temporary code (here we use class 9 that usually refers to “water”). Only the subset of lowest points that receives the temporary classification code 9 will be used for ground classification later.

lasthin -i DenseMatching.laz ^
        -i ManualGround_2m_forest.laz ^
        -drop_class 7 ^
        -merged ^
        -lowest -step 1 -classify_as 9 ^
        -o DenseMatchingAndDensifiedGround.laz

We use the GUI of lasview to pick several interesting areas for visual inspection. The selected points load much faster when the LAZ file is spatially indexed and therefore we first run lasindex. For better orientation we also load the forest stand delineating polygon as an overlay into the GUI.

lasindex -i DenseMatchingAndDensifiedGround.laz 

lasview -i DenseMatchingAndDensifiedGround.laz -gui

We pick the area shown below that contains the target forest with manually collected and densified ground points and a forested area with only dense matching points. The difference could not be more drastic as the visualizations show.

Now we run ground classification using lasground with option ‘-town’ using only the points with the temporary code 9 by ignoring all other classifications 0 and 8 in the file. We leave the temporary classification code 9 unchanged for all the points that were not classified with “ground” code 2 so we can visualize later which ones those are.

lasground -i DenseMatchingAndDensifiedGround.laz ^
          -ignore_class 0 8 ^
          -town ^
          -non_ground_unchanged ^
          -o GroundClassified.laz

We again use the GUI of lasview to pick several interesting areas after running lasindex and again load the forest stand delineating polygon as an overlay into the GUI.

lasindex -i GroundClassified.laz 

lasview -i GroundClassified.laz -gui

We pick the area shown below that contains all three scenarios: the target forest with manually collected and densified ground points, an open area with only dense matching points, and a forested area with only dense matching points. The result is as expected: in the target forest the manually collected ground points are used as ground and in the open area the dense-matching points are used as ground. But there is no useful ground in the other forested area.

Now we can compute the heights of the points above ground for our target forest with lasheight and either replace the z elevations in the file of store them separately as “extra bytes”. Then we can compute, for example, a Canopy Height Model (CHM) that color codes the height of the vegetation above the ground with lasgrid. Of course this will only be correct in the target forest where we have “good” ground but not in the other forested areas. We also compute a hillshaded DTM to be able to visually inspect the topography of the generated terrain model.

lasheight -i GroundClassified.laz ^
          -store_as_extra_bytes ^
          -o GroundClassifiedWithHeights.laz

lasgrid -i GroundClassifiedWithHeights.laz ^
        -step 2 ^
        -highest -attribute 0 ^
        -false -set_min_max 0 25 ^
        -o chm.png

las2dem -i GroundClassified.laz ^
        -keep_class 2 -extra_pass ^
        -step 2 ^ 
        -hillshade ^
        -o dtm.png

Here you can download the resulting color-coded CHM and the resulting hill-shaded DTM as Google Earth KMZ overlays. Clearly the resulting CHM is only meaningful in the target forest where we used the manually collected ground points to create a reasonable DTM. In the other forested areas the ground is only correct near the forest edges and gets worse with increasing distance from open areas. The resulting DTM exhibits some interesting looking  bumps in the middle of areas with manually collected ground point. Those are a result of using the dense-matching points as ground whenever their elevation is lower than that of the manually collected points (which is decided in the lasthin step). Whether those bumps represent true elevations of are artifacts of low erroneous elevation from dense-matching remains to be investigated.

For forests on complex and steep terrain the number of ground points that needs to be manually collected may make such an approach infeasible in practice. However, maybe you have another source of elevation, such as a low-resolution DTM of 10 or 25 meter provided by your local government. Or maybe even a high resolution DTM of 1 or 2 meter from a LiDAR survey you did several years ago. While the forest may have grown a lot in the past years, the ground under the forest will probably not have changed much …

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

LASmoons: Alen Berta

Alen Berta (recipient of three LASmoons)
Department of Terrestrial Ecosystems and Landscape, Faculty of Forestry
University of Zagreb and Oikon Ltd Institute for Applied Ecology, CROATIA

Background:
After becoming the EU member state, Croatia is obliged to fulfill the obligation risen from the Kyoto protocol: National Inventory Report (NIR) of the Green House Gasses according to UNFCCC. One of the most important things during the creation of the NIR is to know how many forested areas there are and their wood stock and increment. This is needed to calculate the size of the existing carbon pool and its potential for sequestration. Since in Croatia, according to legislative, it is not mandatory to calculate the wood stock and yield of the degraded forest areas (shrubbery and thickets) during the creation of the usual forest management plans, this data is missing. So far, only a rough approximation of the wood stock and increment is used during the creation of NIR. However, these areas are expanding every year due to depopulation of the rural areas and the cessation of traditional farming.

very diverse stand structure of degraded forest areas (shrubbery and thickets)

Goal:
This study will focus on two things: (1) Developing regression models for biomass volume estimation in continental shrubberies and thickets based on airborne LiDAR data. To correlate LiDAR data with biomass volume, over 70 field plots with a radius of 12 meters have been established in more than 550 ha of the hilly and lowland shrubberies in Central Croatia and all trees and shrubberies above 1 cm Diameter at Breast Height (DBH) were recorded with information about tree species, DBH and height. Precise locations of the field plots are measured with survey GNNS and biomass is calculated with parameters from literature. For regression modeling, various statistics from the point clouds matching the field plots will be used (i.e. height percentiles, standard deviation, skewness, kurtosis, …). 2) Testing the developed models for different laser pulse densities to find out if there is a significant deviation from results if the LiDAR point cloud is thinner. This will be helpful for planning of the later scanning for the change detection (increment or degradation).

Data:
+
641 square km of discrete returns LiDAR data around the City of Zagreb, the capitol of Croatia (but since it is highly populated area, only the outskirts of the area will be used)
+ raw geo-referenced LAS files with up to 3 returns and an average last return point density of 1 pts/m².

LAStools processing:
1)
extract area of interest [lasclip or las2las]
2) create differently dense versions (for goal no. 2) [lasthin]
3) remove isolated noise points [lasnoise]
4) classify point clouds into ground and non-ground [lasground]
5) create a Digital Terrain Model (DTM) [las2dem]
6) compute height of points above the ground [lasheight]
7) classify point clouds into vegetation and other [lasclassify]
8) normalize height of the vegetation points [lasheight]
9) extract the areas of the field plots [lasclip]
10) compute various metrics for each plot [lascanopy]
11) convert LAZ to TXT for regression modeling in R [las2txt]

Rapidlasso receives “Green Asia Award” at ACRS 2015

PRESS RELEASE (for immediate release)
November 16, 2015
rapidlasso GmbH, Gilching, Germany

At the Asian Conference on Remote Sensing 2015 (ACRS 2015) held in Manila, rapidlasso GmbH was honored with the “Green Asia Award” by the Chinese Society of Photogrammetry and Remote Sensing (CSPRS). This award is given to a paper that directs Asia towards a greener future using remote sensing technology. This year’s award commends rapidlasso GmbH on advancing the area of LiDAR processing through their PulseWaves effort. PulseWaves is a vendor-neutral full waveform LiDAR data exchange format and API that simplifies access to full waveform data and allows researchers to focus on algorithms and share results. In the future this technology may prove valuable to improve biomass estimates for carbon credit programs such as the TREEMAPS project of WWF.

Prof. Kohei Cho and Prof. Peter T. Y. Shih present the award

Prof. Kohei Cho and Prof. Peter T. Y. Shih present the Green Asia Award

The society communicated to Dr. Martin Isenburg, CEO of rapidlasso GmbH, that this award was also meant to honor his many years of teaching and capacity building across the Asian region. Since the beginning of 2013 rapidlasso GmbH has conducted well over 50 seminars, training events, and hands-on workshops at universities, research institutes, and government agencies in Thailand, Malaysia, Myanmar, Vietnam, Indonesia, Singapore, Taiwan, Japan, and the Philippines. The on-going LiDAR teaching efforts of rapidlasso GmbH in Asia and elsewhere can be followed via their event page.

Green Asia Award for CEO of rapidlasso GmbH

Green Asia Award given to the CEO of rapidlasso GmbH

The award certificate that was presented to Dr. Martin Isenburg by Prof Kohei Cho and Prof Peter Shih during the closing ceremony of ACRS 2015 came with a cash reward of USD 300. The award money was donated to the ISPRS summer school that followed the ACRS conference to top off the pre-existing “green sponsorship” by rapidlasso GmbH that was already supporting a “green catering” of summer school lunches and dinners to avoid single-use cups, plastic cutlery and styrofoam containers. The additional award money was used for hosting the main summer school dinner at a sustainable family-run restaurant serving “happy chickens” and “happy pigs” raised organically on a local farm.

during the closing ceremony of ACRS 2015

Award Ceremony held during Closing of ACRS 2015

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.

Two ASPRS awards for “pit-free” CHM algorithm

PRESS RELEASE (for immediate release)
July 29, 2015
rapidlasso GmbH, Gilching, Germany

The paper “Generating Pit-free Canopy Height Models from Airborne LiDAR” co-authored by rapidlasso GmbH and published in the September 2014 issue of PE&RS (the journal of the ASPRS) was awarded twice at the IGTF 2015 – ASPRS Annual Conference in Tampa, Florida last May. The paper took home the John I. Davidson President’s Award for Practical Papers (2nd Place) as well as the Talbert Abrams Award (2nd Honorable Mention).

The John I. Davidson President’s Award for Practical Papers (2nd Place).

The “pit-free” CHM paper wins the John I. Davidson President’s Award for Practical Papers (2nd Place) and the Talbert Abrams Award (Second Honorable Mention).

The “pit-free” CHM paper is joint work with Anahita Khosravipour, Andrew K. Skidmore, Tiejun Wang, and Yousif A. Hussin of ITC and University of Twente. It describes a technique that can create raster Canopy Height Models (CHMs) without the so called “pits” that tend to hamper subsequent extraction of individual tree attributes such as number, location, height, and crown diameter. The paper uses data measured in the field by ITC researchers to show that “pit-free” CHMs significantly lower the commission and omission errors in single tree detection.

Side-by-side comparison of a "standard" CHM and a "pit-free" CHM.

Visual side-by-side comparison of a “standard” versus a “pit-free” CHM.

The “pit-free” CHM algorithm can easily be implemented with LAStools either by modifying an available batch script or by executing the LAStools Pipelines distributed with the toolboxes for ArcGIS and QGIS. A detailed blog article that compares various different methods for creating CHMs is available via the Web pages of rapidlasso GmbH.

We at rapidlasso GmbH like to especially congratulate the main author, Ms. Anahita Khosravipour, who managed to get two awards with her very first academic publication. Those who like our “pit-free” CHM algorithm will probably also love the new technique that our team will introduce later this year at SilviLaser 2015 in France.

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

Andrew Smith (recipient of three LASmoons)
Advanced Diploma in GIS Applications
Vancouver Island University, Victoria, CANADA

Background:
Forestry and natural resource development are a large part of British Columbia’s economy and tourism. BC has some of the largest areas of forest-covered land in Canada, and some of the toughest terrain, which makes it a great place for towering trees and nature lovers. With proper resource management we can help ensure that the forestry sector meets its needs for wood while maintaining its obligations to ecosystem management and land use planning.

lasmoons_andrew_smith_0Goal:
The aim of this project is to assess the efficiency of utilising LiDAR data to develop a more spatially
detailed and accurate classification system of tree species. By identifying key tree species, and
tree populations, we can help forestry companies to be more efficient, and cost effective, allowing them
to have better ecosystem management and land use planning. This will help with harvesting profiles,
cutting schedules, opening size, as well as proper habitat analysis and development in line with
BC Forests and Range Practices Act, which spells out forestry obligations for wildlife management and
habitat creation.

Data:
+ 5 square kilometres of LiDAR data covering BC Coastal Forest on Vancouver island, Canada.
+ Average point density: 37 points per square metre.

LAStools processing:
1)
classify the ground points to create a DTM [lasground]
2) normalize the height the trees. the tallest tree in the sample plot is 70 meters tall [lasheight]
3) generate a Canopy Height Model (CHM) using CHM using the pit-free algorithm from (Khosravipour et. al, 2014) with the workflow described here [lasthin, las2dem, lasgrid]
4) overlay the CHM with high-resolution ortho images to first extract tree species, location and canopy diameter and to then derive tree diameter at breast height (DBH) for all identificatied single trees using Taper equations.

Reference:
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.
BC Forests and Range Practices Act, Forest and Range Evaluation Program Strategic Plan 2011 – 2013 November 8, 2011. http://www.for.gov.bc.ca/ftp/hfp/external/!publish/frep/library/FREP-Strategic-Plan-2011.pdf

Beautiful Full-Waveform LiDAR of a Tropical Rainforest

Scanning forests with LiDAR is done to predict timber yield, estimate biomass, and monitor ecosystems. A detailed Digital Terrain Model (DTM) is needed to accurately compute forestry metrics such as tree height. For generating a DTM it is important to have enough “ground returns” that measure the elevation of the bare-earth for which the laser pulse needs to penetrate through the vegetation down to the forest floor and reflect back to the plane. This can be difficult in a dense tropical forest with tall trees and many layers of canopy as the light of the laser pulse has many opportunities to be reflected or absorbed before reaching the ground. Often none of the light energy will reach the forest floor or the returning ground echo will be too weak to be registered back at the plane.

The "airborne sensors" of  Asian Aerospace Services is equpiied with a RIEGL Q680i airborne LiDAR scanner.

The “airborne sensors” aircraft of Asian Aerospace Services is equipped with a RIEGL Q680i airborne LiDAR scanner.

Asian Aerospace Services Ltd. and rapidlasso GmbH flew a RIEGL LMS Q680i LiDAR system out of a Diamond DA42 “Airborne Sensors” aircraft above dense tropical rainforest in Thailand. We did a short test flight to determine whether the LiDAR would be able to penetrate the canopy of Khao Yai National Forest. The original objective was to demonstrate that the scanner is able to “see through” the canopy in order to produce data for the TREEMAPS project of WWF Thailand funded with EUR 1,730,205 through the International Climate Initiative of the German government. Unfortunately Germany had to pull the funding for TREEMAPS after political tinkering started to interfere with technical objectives, and the LiDAR produced by our “pro-bono flight” has become the only tangible outcome produced by TREEMAPS. We are now making this data available under a Creative Commons Public Domain license (for download see below).

Above Khao Yao National Forest scanning for the pro-bono test flight to investigate the canopy penetration of our LiDAR.

Above Khao Yao National Forest during a test scan to investigate the canopy penetration of the RIEGL Q680i LiDAR scanner.

Low clouds disrupted our original plans but we had one cloud-free stretch during a ferry flight between target locations. We were scanning for 90 seconds at 1200 meters above ground with a speed of 220 km/h. The Q680i was set to operate at a laser pulse rate of 80 kHz with the polygon mirror rotating 10 times per second. The effective scan rate was 240 / 360 * 80,000 = 53,333 shots per second, as each of the 4 mirror facets covers 60 degrees – summing up to 240 degrees – and there are no laser shots leaving the aircraft for the remaining 120 degrees of between facets.

The 53,333 shots per second are reflected over the length of 10 * 4 = 40 mirror facets so that a single scan line contains a sequence of 53,333 / 40 = 1333 shots spaced about 12.4 microseconds apart. After each scan line comes a short pause of 8.3 milliseconds during which the mirror rotates 30 degrees to the next facet. Flying at 220 km/h and 1200 meters above ground means that each second we scanned 62 meters of a 1430 meters wide swath. This resulted in a pulse spacing of 62 / 40 = 1.55 meters between subsequent scan lines and a pulse spacing of 1430 / 1333 = 1.07 meters along each scan line. The 90 second scan resulted in a strip length of about 5.5 km.

The PulseWaves export options available in RIEGL's RiPROCESS.

The PulseWaves export options available in RIEGL’s RiPROCESS.

For each shot the outgoing as well as the returning waveform are digitized and stored to RIEGL’s proprietary SDF format. We used RIEGL’s RiPROCESS software (version 1.5.8) to export the 4,770,152 shots as full waveform LiDAR to PulseWaves and to extract a total of 7,956,587 returns as LASzip-compressed discrete LiDAR to the LAZ format.

The DSM, DTM, and CHM of the 5.5 km long strip scanned above Khao Yai National Forest.

The DSM, DTM, and CHM of the 5.5 km long strip scanned above Khao Yai National Forest.

The Q680i produced sufficiently many ground returns to compute a “plausible” DTM (using lasground and las2dem of LAStools) despite the dense 45 meter tropical canopy during leaf-on season and high humidity. Typical terrain features such as streams and slopes are clearly visible in the hillshading. The corresponding Canopy Height Model (CHM) uses a false coloring that maps vegetation heights of 0 meters to blue and heights of 45 meters to red. There is less vegetation (darker blue) in areas that appear to be in the steeper parts of the DTM.

Only 5 % of pulses emitted produced a discrete return on the ground … this is rather low. It means that the generated DTM will be of much coarser resolution than the number of pulses shot per square meter would typically suggest. Of the 250,473 points classified as ground only 30,931 are single returns where the laser had an unobstructed view of the bare earth. More than half of the ground returns come from laser shots that had already produced two or more vegetation returns and 7,665 of them had their light energy weakened by grazing four or more layers before hitting the forest floor. Scrutinizing the full waveform LiDAR might lead to more ground information to improve the accuracy of the DTM.

We are happy to provide both – the discrete return as well as the full-waveform LiDAR – under a Creative Commons Public Domain license to researchers world-wide. You can download the data here:

You can decompress the pair of compressed PulseWaves files PLZ + WVZ into an pair of uncompressed PulseWaves files PLS + WVS with pulsezip.exe – a recent addition to PulseTools. We would be happy to hear back about your experiences with the data. You could join the PulseWaves user forum at http://pulsewaves.org and share your insights with the group. Below a visualization of a small part of the full waveforms …

Full Waveform LiDAR from Khao Yai National forest.

Full Waveform LiDAR from Khao Yai National forest.