Category Archives: spatial indexing

National Open LiDAR Strategy of Latvia humiliates Germany, Austria, and other European “Closed Data” States

Posted on by

Latvia, officially the Republic of Latvia, is a country in the Baltic region of Northern Europe has around 2 million inhabitants, a territory of 65 thousand square kilometers and – since recently – also a fabulous open LiDAR policy. Here is a list of 65939 tiles in LAS format available for free download that cover the entire country with airborne LiDAR with a density from 4 to 6 pulses per square meters. The data is classified into ground, building, vegetation, water, low noise, and a few other classifications. It is licensed Creative Commons CC0 1.0 – meaning that you can copy, modify, and distribute the data, even for commercial purposes, all without asking permission. And there is a simple and  functional interactive download portal where you can easily download individual tiles.

latvia_open_data_portal_01

Interactive open LiDAR download portal of Latvia.

We downloaded the 5 by 5 block of square kilometer tiles matching “4311-32-XX.las” for checking the quality and creating a 1m DTM and a 1m DSM raster. You can follow along after downloading the latest version of LAStools.

Quality Checking

We first run lasvalidate and lasinfo on the downloaded LAS files and then immediately compress them with laszip because multi-core processing of uncompressed LAS files will quickly overwhelm our file system, make processing I/O bound, and result in overall longer processing times with CPUs waiting idly for data to be loaded from the drives.

lasinfo -i 00_tiles_raw\*.las ^
        -compute_density ^
        -histo z 5 ^
        -histo intensity 256 ^
        -histo user_data 1 ^
        -histo scan_angle 1 ^
        -histo point_source 1 ^
        -histo gps_time 10 ^
        -odir 01_quality -odix _info -otxt ^
        -cores 3
lasvalidate -i 00_tiles_raw\*.las ^
            -no_CRS_fail ^
            -o 01_quality\report.xml

Despite already excluding a missing Coordinate Reference System (CRS) from being a reason to fail (the lasinfo reports show that the downloaded LAS files do not have any geo-referencing information) lasvalidate still reports a few failing files, but scrutinizing the resulting XML file ‘report.xml’ shows only minor issues.

Usually during laszip compression we do not alter the contents of a file, but here we also add the EPSG code 3059 for CRS “LKS92 / Latvia TM” as we turn bulky LAS files into slim LAZ files so we don’t have to specify it in all future processing steps.

laszip -i 00_tiles_raw\*.las ^
       -epsg 3059 ^
       -cores 2

Compression reduces the total size of the 25 tiles from over 4.1 GB to below 0.6 GB.

Next we use lasgrid to visualize the last return density which corresponds to the pulse density of the LiDAR survey. We map each 2 by 2 meter pixel where the last return density is 2 or less to blue and each 2 by 2 meter pixel it is 8 or more to red.

lasgrid -i 00_tiles_raw\*.laz ^
        -keep_last ^
        -step 2 ^
        -density_16bit ^
        -false -set_min_max 2 8 ^
        -odir 01_quality -odix _d_2_8 -opng ^
        -cores 3

This we follow by the mandatory lasoverlap check for flight line overlap and alignment where we map the number of overlapping swaths as well as the worst vertical difference between overlapping swaths to a color that allows for quick visual quality checking.

lasoverlap -i 00_tiles_raw\*.laz ^
           -step 2 ^
           -min_diff 0.1 -max_diff 0.2 ^
           -odir 01_quality -opng ^
           -cores 3

The results of the quality checks with lasgrid and lasoverlap are shown below.

last return density: blue = 2 per square meter, red = 8 per square meter
flight line overlap: blue – 1 flight line, aqua – 2 flight lines, yellow – 3 flight lines, orange – 4 flight lines, red – 5 or more flight lines
flight line alignment: white – error less than 10 cm, red or blue – error more than 20 cm

Raster Derivative Generation

Now we use first las2dem to create a Digital Terrain Model (DTM) and a Digital Surface Model (DSM) in RasterLAZ format and then use blast2dem to create merged and hill-shaded versions of both. Because we will use on-the-fly buffering to avoid edge effects along tile boundaries we first spatially index the data using lasindex for more efficient access to the points from neighboring tiles.

Resulting hill-shaded DTM.
Resulting hill-shaded DSM. 
lasindex -i 00_tiles_raw\*.laz ^
         -cores 3

las2dem -i 00_tiles_raw\*.laz ^
        -keep_class 2 9 ^
        -buffered 25 ^
        -step 1 ^
        -use_orig_bb ^
        -odir Latvia\02_dtm_1m -olaz ^
        -cores 3

blast2dem -i 02_dtm_1m\*.laz ^
          -merged ^
          -hillshade ^
          -step 1 ^
          -o dtm_1m.png

las2dem -i 00_tiles_raw\*.laz ^
        -drop_class 1 7 ^
        -buffered 10 ^
        -spike_free 1.5 ^
        -step 1 ^
        -use_orig_bb ^
        -odir 03_dsm_1m -olaz ^
        -cores 3

blast2dem -i 03_dsm_1m\*.laz ^
          -merged ^
          -hillshade ^
          -step 1 ^
          -o dsm_1m.png

Because the overlaid imagery does not look as nice in our new Google Earth installation, below are the DTM and DSM at versions down-sampled to 25% of their original size.

hill-shaded DTM down-sampled to 25%
hill-shaded DSM down-sampled to 25%

Many thanks to SunGIS from Latvia who tweeted us about the Open LiDAR after we chatted about it during the Foss4G 2019 gala dinner. Kudos to the Latvian Geospatial Information Agency (LGIA) for implementing a modern national geospatial policy that created opportunity for maximal return of investment by opening the expensive tax-payer funded LiDAR data for re-purposing and innovation without barriers. Kudos!

Plots to Stands: Producing LiDAR Vegetation Metrics for Imputation Calculations

Posted on by

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.

flightline overlap: blue=1, cyan=2, yellow=3, orange=4, red=5
flightline differences: vertical up to 20 cm = white, above 40 cm = saturated red or blue 
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?

pulse density per 2 by 2 meter: blue=20
pulse density per 5 by 5 meter: blue=20

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 lastile, lasnoise, 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
vendor-delivered ground classification
new ground classification from LAStools

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
maximum metric is sensitive to high noise
95th percentile is more robust against noise
50th percentile (of points with heights above cutoff)
25th percentile (of points with heights above cutoff)
percentage of points with heights between 2 and 10 meters (of those above cutoff)
percentage of points with heights between 10 and 100 meters (of those above cutoff)
30 percent bincentile
80 percent bincentile
canopy cover above 3 meters

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

NRW Open LiDAR: Merging Points into Proper LAS Files

Posted on by

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.

Tile with classified returns and point flagged as synthetic.
Recovered return types: single, first of two, last of two.
Synthetic points visible under the Kennedy bridge.
Synthetic points also visible under the Opera house.

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

RIEGL Becomes LASzip Sponsor for LAS 1.4 Extension

Posted on by 
PRESS RELEASE (for immediate release)
August 31, 2015
rapidlasso GmbH, Gilching, Germany

We are happy to announce that RIEGL Laser Measurement Systems, Austria has become a sponsor of the award-winning LASzip compressor. Their contribution at the Silver level will kick-off the actual development phase of the “native LAS 1.4 extension” that had been discussed with the LiDAR community over the past two years. This “native extension” for LAS 1.4 complements the existing “compatibility mode” for LAS 1.4 that was supported by Gold sponsor NOAA and Bronze sponsors Quantum Spatial and Trimble Geospatial. The original sponsor who initiated and financed the open sourcing of the LASzip compressor was USACE – the US Army Corps of Engineers (see http://laszip.org).

The existing “LAS 1.4 compatibility mode” in LASzip was created to provide immediate support for compressing the new LAS 1.4 point types by rewriting them as old point types and storing their new information as “Extra Bytes”. As an added side-benefit this has allowed legacy software without LAS 1.4 support to readily read these newer LAS files as most of the important fields of the new point types 6 to 10 can be mapped to fields of the older point types 1, 3, or 5.

In contrast, the new “native LAS 1.4 extension” of LASzip that is now sponsored in part by RIEGL will utilize the “natural break” in the format due to the new point types of LAS 1.4 to introduce entirely new features such as “selective decompression”, “rewritable classifications and flags”, “integrated spatial indexing”, … and other functionality that has been brain-stormed with the community since rapidlasso GmbH had issued the open “call for input” on native LASzip compression for LAS 1.4 in January 2014. We invite you to follow the progress or contribute to the development via the discussions in the “LAS room“.

silverLASzip_m60_512_275

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.

About RIEGL:
Austrian based RIEGL Laser Measurement Systems is a performance leader in research, development and production of terrestrial, industrial, mobile, bathymetric, airborne and UAS-based laser scanning systems. RIEGL’s innovative hard- and software provides powerful solutions for nearly all imaginable fields of application. Worldwide sales, training, support and services are delivered from RIEGL‘s Austrian headquarters and its offices in Vienna, Salzburg, and Styria, main offices in the USA, Japan, and in China, and by a worldwide network of representatives covering Europe, North and South America, Asia, Australia and Africa. Visit http://riegl.com for more information.

Use Buffers when Processing LiDAR in Tiles !!!

Posted on by

We often process LiDAR in tiles for two reasons: first, to keep the number of points per file low and use main memory efficient, and second, to speed up the computation with parallel tile processing and keep all cores of a modern CPU busy. However, it is very (!!!) important to take the necessary precautions to avoid “edge artifacts” when processing LiDAR in tiles. We have to include points from neighboring tiles during certain LAStools processing steps to avoid edge artifacts. Why? Here is an illustration from our PHIL LiDAR tour earlier this year:

Buffers are important to avoid edge artifacts along tile boundaries during DTM creation.

Buffers are important to avoid edge artifacts along tile boundaries.

What you see is the temporary TIN of ground points created (internally) by las2dem or blast2dem that is then rastered at the user-specified step size onto a grid. Without a buffer (right side) there will not always be a triangle to cover every pixel. Especially in the corners of the tile you will often find empty pixels. Furthermore the poorly shaped “sliver triangles” along the boundary of the TIN do not interpolate the ground elevations properly. In contrast, with buffer (left side) the TIN generously covers the entire area that is to be rastered with nicely shaped triangles.

The

Christmas cookie analogy: buffers are like generously rolling out the dough

Here the christmas cookies analogy: You need to roll out the dough larger than the cookies you will cut to make sure your cookies will have nice edges. Think of the TIN as the dough and the square tile as your cookie cutter. You need to use a sufficiently large buffer when you roll out your TIN to assure an edge without crumbles when you cut out the tile … (-: … otherwise you are pretty much guaranteed to get results that – upon closer inspection – have these kind of artifacts:

Without buffers processing artifacts also happen when classifying points with lasground or lasclassify, when calculating height above ground or height-normalizing LiDAR tiles with lasheight, when removing noise with lasnoise, when creationg contours with las2iso or blast2iso, or any other operation where an incomplete neighborhood of points can affect the results. Hence, we need to surround each tile with a temporary buffer of points. Currently there are two ways of working with buffers with LAStools:

  1. creating buffered tiles during the initial tiling step with the ‘-buffer 25’ option of lastile, maintaining buffered tiles throughout processing and finally using the ‘-use_tile_bb’ option of lasgrid, las2dem, blast2dem, or lascanopy to raster the tiles without the temporary buffer.
  2. creating buffered tiles from non-overlapping (= unbuffered) tiles with “on-the-fly” buffering using the ‘-buffered 25’ option of most LAStools such as lasground, lasheight, or las2dem. For some workflows it is useful to also add ‘-remain_buffered’ if buffers are needed again in the next step. Finally, we use the ‘-use_orig_bb’ option of lasgrid, las2dem, blast2dem, or lascanopy to raster the tiles without the temporary buffer.

In the following three (tiny) examples using the venerable ‘fusa.laz’ sample that is distributed with LAStools to illustrate the two types of buffering as well as to show what happens when no buffers are used. In each example we will first cut the small ‘fusa.laz’ sample into nine smaller tiles and then process these separately on 4 cores in parallel.

1. Initial buffer creation with lastile

This is what most of my tutorials teach. It assumes you are the one creating the tiling in the first place. If you do it with lastile and add a buffer right from the start things are pretty easy.

lastile -i ..\data\fusa.laz ^
        -set_classification 0 -set_user_data 0 ^
        -tile_size 100 -buffer 20 ^
        -odir 1_raw -o futi.laz

We cut the input into 100 meter by 100 meter tiles but add a 20 meter buffer around each tile. That means that each tile on disk will contain the points for an area of up to 140 meter by 140 meter. The GUI for LAStools shows the overlap and if you scrutinize the bounding box values that the cursor points to you notice the extra 20 meters in each direction.

tiles_buffered_with_lastile

Now we can forget about the buffers and run the standard workflow consiting of lasground, lasheight, and lasclassify to distinguish ground, vegetation, and building points in the LiDAR tiles.

lasground -i 1_raw\futi*.laz ^
          -city ^
          -odir 1_ground -olaz ^
          -cores 4
lasheight -i 1_ground\futi*.laz ^
          -drop_above 50 ^
          -odir 1_height -olaz ^
          -cores 4
lasclassify -i 1_height\futi*.laz ^
            -odir 1_classify -olaz ^
            -cores 4

At the end – when we generate raster products – we have to remember that the tiles were buffered by lastile and cut off the buffers when we raster the TIN with option ‘-use_tile_bb’ of las2dem.

las2dem -i 1_classify\futi*.laz ^
        -keep_class 2 6 ^
        -step 0.25 -use_tile_bb ^
        -odir 1_dbm -obil ^
        -cores 4

We created a digital terrain model with buildings (DBM) by keeping the points with classification 2 (ground) and 6 (building). After loading the resulting 9 tiles into QGIS and generating a virtual raster we see a nice seamless DBM without any edge artifacts.

The DEM of the 9 tiles computed with buffers created by lastile has no edge artifacts acoss tile boundaries.

The DBM of the 9 tiles computed with buffers created by lastile has no edge artifacts acoss tile boundaries.

If you need to deliver the LiDAR files you should remove the buffers with lastile and option ‘-remove_buffer’.

lastile -i 1_classify\futi*.laz ^
        -remove_buffer ^
        -odir 1_final -olaz ^
        -cores 4

2. On-the-fly buffering

Now assume you are given LiDAR tiles without buffers. We generate them here with lastile.

lastile -i ..\data\fusa.laz ^
        -set_classification 0 -set_user_data 0 ^
        -tile_size 100 ^
        -odir 2_raw -o futi.laz

The only difference is that we do not request the 20 meter buffer and the result is a typical tiling as you may receive it from a vendor or download it from a LiDAR portal. The GUI for LAStools shows that there is no overlap and if you scrutinize the bounding box values that the cursor points to, you see that the tiles is exactly 100 meters by 100 meters.

tiles_without_buffer

Now we have to think about buffers a lot. When using on-the-fly buffering we should first spatially index the tiles with lasindex for faster access to the points from neighbouring tiles.

lasindex -i 1_raw\futi*.laz -cores 4

Below in red are the modifications for on-the-fly buffering to the standard workflow of lasground, lasheight, and lasclassify. The first lasground run uses ‘-buffered 20’ to add buffers to each tile and ‘-remain_buffered’ to write those buffers to disk. This way they do not have to created again by lasheight and lasclassify.

lasground -i 2_raw\futi*.laz ^
          -buffered 20 -remain_buffered ^
          -city ^
          -odir 2_ground -olaz ^
          -cores 4
lasheight -i 2_ground\futi*.laz ^
          -remain_buffered ^
          -drop_above 50 ^
          -odir 2_height -olaz ^
          -cores 4
lasclassify -i 2_height\futi*.laz ^
            -remain_buffered ^
            -odir 2_classify -olaz ^
            -cores 4

At the end we have to remember that the tiles still have on-the-fly buffers and them cut off with option ‘-use_orig_bb’ of las2dem.

las2dem -i 2_classify\futi*.laz ^
        -keep_class 2 6 ^
        -step 0.25 -use_orig_bb ^
        -odir 2_dbm -obil ^
        -cores 4

Again, we created a digital terrain model with buildings (DBM) by keeping the points with classification 2 (ground) and 6 (building). The resulting hillshade computed from a virtual raster that combines the 9 BIL rastera into one looks perfectly smooth in QGIS.

The hillshaded DEM of the 9 tiles computed with on-the-fly buffering has no edge artifacts acoss tile boundaries.

The hillshaded DBM of 9 tiles computed with on-the-fly buffering has no edge artifacts acoss tile boundaries.

If you need to deliver the LiDAR files you should probably remove the buffers first … (-:

lastile -i 2_classify\futi*.laz ^
        -remove_buffer ^
        -odir 2_final -olaz ^
        -cores 4

3. Bad: No buffering

Here what you are *not* supposed to do. Assuming you get unbuffered tiles.

lastile -i ..\data\fusa.laz ^
        -set_classification 0 -set_user_data 0 ^
        -tile_size 100 ^
        -odir 3_raw -o futi.laz

Bad. You do not take care about buffering when processing the tiles.

lasground -i 3_raw\futi*.laz ^
          -city ^
          -odir 3_ground -olaz ^
          -cores 4
lasheight -i 3_ground\futi*.laz ^
          -drop_above 50 ^
          -odir 3_height -olaz ^
          -cores 4
lasclassify -i 3_height\futi*.laz ^
            -odir 3_classify -olaz ^
            -cores 4

Bad. You do not take care about buffering when generating the DBM.

las2dem -i 3_classify\futi*.laz ^
        -keep_class 2 6 ^
        -step 0.25 ^
        -odir 3_dbm -obil ^
        -cores 4

Bad. You get crappy results with edge artifacts clearly visible in the hillshade.

The hillshaded DBM of 9 tiles computed WITHOUT using buffers has severe edge artifacts acoss tile boundaries.

The hillshaded DBM of 9 tiles computed WITHOUT using buffers has severe edge artifacts acoss tile boundaries.

Bad. If you zoom in on a corner where 4 tiles meet you find missing pixels and incorrect elevation values. Bad. Bad. Bad. So please folks. Try this on your own data. Notice the horrible edge artifacts. Then always use buffers … (-:

PS: Usually no buffers are needed for running lasgrid, lasoverlap, or lascanopy as they perform simple binning operations that do not make use of neighbour information.

Preparing raw LiDAR for efficient (online) distribution

Posted on by

On August 14th, Prof. David Pyle tweeted about raw LiDAR being publically available for the vulcanic island Nea Kameni and its little sibling Palea Kameni that are part of the Santorini caldera in Greece. A big “kudos” to all those like Prof. Pyle who share ‎‎raw LiDAR‬ data online for all to download – be it for transparent research, as an open data policy, or to enable innovation. Although fancy download portals like OpenTopography are great, already a folder full of files accessible via simple FTP or HTTP is an incredible resource. For the latter, here are some tips from rapidlasso how to prepare your raw LiDAR flightlines with LAStools so that they are as good as they can be for those that download them … (-:

The twelve LiDAR flight lines for the Kameni islands are a perfectly-sized example on how to prepare raw LiDAR for efficient (online) distribution. They are provided both in the LAS format as well as simple ASCII files. Stored in the LAS format, the twelve strips (1,2,3,4,5,6,7,8,9,10,11,12)  are about 1 Gigabyte:

lidar_for_download_santorini_size_raw

A lasinfo report tells us the following:
(1) The exporting software used the wrong GeoTIFF tag to specify the UTM 35 (north) projection via EPSG code 32635.
(2) There are four proprietary “LeicaGeo” VLRs totalling 22000 bytes stored in each header. Does someone know what they contain and how to read them?
(3) There are two legacy bytes following the header that are part of the (now somewhat dated) LAS 1.0 specification.
(4) The coordinates are stored with millimeter resolution (i.e. the scale factors are 0.001). This is an overkill for airborne LiDAR. Those millimeters are just scanning noise, are miss-leading, and negatively affect compression.
(5) The file does not store flight line numbers in the “file source ID” field.

lidar_for_download_santorini_lasinfo_raw

Visualizing the twelve flight lines with lasview further illustrates:
(1) A cruise ship and a cloud was captured. the latter is classified as noise (7).
(2) An intensity value is stored for each return.
(3) Almost all laser shots resulted in a single return, which is not surprising for a vulcanic island without vegetation. The multi-returns from the cloud are a colorful exception.
(4) There are some strange (=> useless) numbers stored in the “point source ID” field of each point that should really store the flight line number of each point.

We suggest to fix up this raw LiDAR to make it more efficient and useful for those receiving it as follows. We run las2las to do the following:
(1) Remove the legacy two bytes following the header.
(2) Switch to the 1.2 version of the LAS format.
(3) Remove all existing VLRs.
(4) Set the horizontal projection to EPSG code 32635 and the vertical datum to WGS84.
(5) Set the file source ID of each file and the point source ID of each point to the same value. This values starts at one and is incremented with each file.
(6) Rescale coordinate resolution to centimeters.
(7) Append the string ‘_cm’ to the original file name to output LASzip-compressed files that end in ‘*.laz’.

las2las -i *.LAS ^
        -remove_extra ^
        -set_version 1.2 ^
        -remove_all_vlrs ^
        -epsg 32635 -vertical_wgs84 ^
        -files_are_flightlines ^
        -rescale 0.01 0.01 0.01 ^
        -odix _cm -olaz

This result in much much smaller files that are easier to host and faster to download. We achieve almost a factor 10 in compression as the new files are only 11.2 % of the size of the original files.

lidar_for_download_santorini_size_fixed

We also get a much cleaner raw LiDAR file that has only meaningful VLRs, more suitable coordinate resolution, and properly populated flight line information:

lidar_for_download_santorini_lasinfo_fixed

And we now have flightline information that tells us, for example, that the (blue) cloud was captured by the third flight line.

lidar_for_download_santorini_05

We noticed something odd when loading the files with lasview: Although most flightlines were still in acquisition order (i.e. the points in the file are in the order acquired by the scanner), flightline 8 and 10 were not. Perfectionists like us employ lassort with option ‘-gps_time’ and reorder all twelve files by GPS time on 4 cores in parallel.

lassort -i LDR*cm.laz ^
        -gps_time ^
        -odix _sort -olaz ^
        -cores 4

Several calls to lasdiff – a tool that reports any content and order difference between two files – confirm that only points of flightline 8 and 10 are actually reordered by the above call. Here two example calls:

lidar_for_download_santorini_lasdiff_sorted

One – rather new – thing that we are now recommend doing is to also add spatial indexing information and do this by storing the little LAX files directly into the compressed LAZ files. This can be done in-place with

lasindex -i LDR*cm.laz ^
         -append ^
         -cores 8

Such spatial indexing information (see the video of our ELMF 2012 talk for more details) allows faster spatial queries that only read and decompress the actually queried area-of-interest. This is already used by all LAStools and can be exploited via the LASlib application programming interface. Soon there will also be support for spatially-indexed queries in the LASzip compression DLL.

lidar_for_download_santorini_06

Finally, should you have ortho photo imagery in TIF format (like it was kindly provided here by Prof. David Pyle) then you can embellish your LiDAR points with RGB using lascolor as shown here:

lascolor -i LDR*_cm.laz ^
         -image KameniOrtho.tif ^
         -odix _col -olaz ^
         -cores 4

Now you have raw LiDAR strips for online distribution that are as good as they get … (-:

lidar_for_download_santorini_07

Esri and rapidlasso develop joint LiDAR compressor

Posted on by
PRESS RELEASE (April Fools’ Day)
April 1, 2014
rapidlasso GmbH, Gilching, Germany

In a positive spin of events, Esri and rapidlasso are announcing to join forces and together develop a LiDAR compressor for LAS 1.4 in open source avoiding unnecessary format fragmentation. Their new “LASeasy” tool not only compresses but also optimizes LAS files for efficient area-of-interest queries. LASeasy extends the popular LASzip compressor to handle LAS 1.4 content and includes the tiny spatial indexing *.lax files into the *.laz file via Extended Variable Length Records (EVLRs). More importantly, LASeasy provides new features such as optional spatial sorting and precomputed statistics – motivated by Esri – that allow exploiting LiDAR in the cloud.

To minimize disruption in existing workflows, their joint effort uses a clever strategy that capitalizes on the natural “break” in the ASPRS LAS format from version 1.3 to 1.4. LAS files compressed by Esri will automatically be upgraded to the new point types introduced with LAS 1.4 (and be losslessly downgraded on decompression). LiDAR software already supporting LAZ will instantly be able to read all LiDAR produced by Esri with the same DLL update that will be needed to access future compressed LAS 1.4 content – achieving maximum compatibility with minimal disruption for users of ArcGIS, LASzip, and the larger LiDAR community,

Martin Isenburg, chief scientist and CEO of rapidlasso GmbH, was all smiles during the announement. “Yes, I had some hard feelings when hearing about their ‘LAZ clone‘ because our presumed open dialogue suddenly felt so very one-sided,” he said, “So over Martin Luther King weekend I proposed this LAS 1.4 trick as a joint development quoting MLK’s ‘We must accept finite disappointment, but never lose infinite hope’ and that seemed to resonate with them.” Speaking on the condition of anonymity an executive of Esri’s management added “For a global geospatial player like us it can happen that we do something ‘evil-by-accident’. We occasionally need someone like Martin to poke some good-natured fun at Esri to remind us of our values.”

LASeasy optimizes LAS files by reordering points along an adaptive space-filling curve for efficient LiDAR queries in the cloud. To access the corner of the LiDAR tile only the points shown in blue need to be loaded and decompressed.

LASeasy optimizes LAS files by reordering points along an adaptive space-filling curve for efficient LiDAR queries in the cloud. To access the corner of the LiDAR tile only the points shown in blue need to be loaded and decompressed.

Warming up for ILMF 2014, rapidlasso puts lean yet plush “LASlayers” on LiDAR

Posted on by
PRESS RELEASE
February 14, 2014
rapidlasso GmbH, Gilching, Germany

As a sweet foretaste to ILMF 2014, the creators of LAStools, LASzip, and PulseWaves are announcing “LASlayers” already on Valentine’s Day. The new functionality nicely complements their popular and widely-used LiDAR compressor making the compressed LAZ files editable for most practical purposes. LASlayers significantly reduce I/O load for writing modification to LAS or LAZ files, especially when batch-processing LiDAR tiles on many cores in parallel or when sending changes to LiDAR files across bandwidth-limited networks.

Conceptually LASlayers add

LASlayers store modifications and additional  attributes to raw LAS or LAZ files in small LAY files avoiding to replicate data that has not changed.

Most point attributes (e.g. coordinates, intensities, scan angles, GPS times, etc.) are not modified when processing LiDAR. LASlayers avoids re-writing the unchanged portions of LAS or LAZ by storing only the actual differences layer by layer to a new “LAY” file. Changing the point classifications or deleting a few points, for example, can be done with LAY files that are just a tiny fraction of the size of a new LAS or LAZ file. Adding new attributes such as RGB colors or the height-above-ground with LASlayers means only this new information needs to be written.

This also provides simultaneous access to different versions of the data: a LiDAR server or a Web portal can store only a single copy of the raw LiDAR and apply LASlayers as needed on-the-fly, for example, to replace ellipsoidal with orthometric elevations or to add RGB colors.

Even users of other LiDAR processing software can readily take advantage of LASlayers with the new “laslayers.exe” tool that computes the difference between a raw and a modified LAS or LAZ file and expresses it as a LAY file (assuming the point order has not changed). A typical use case is the exchange of modifications to LiDAR files between remote locations such as a vendor in Australia or Canada and a data processing center in China or India. Instead of up- and downloading the entire modified LAS or LAZ files, only the much smaller LAY files need to be send back and forth.

A fully featured prototype of LASlayers is available (10 MB including data) together with three simple exercises that illustrate the concept and allow interested parties to test it already today on their own data.

new compressed LAS format by ESRI

Posted on by

NEWSFLASH: update on Jan 7th, 12th, 19th, and Feb 7th (see end of article)

Today I got an email from a LAStools user at NOAA pointing out a new entry in the ArcGIS 10.2 documentation of ESRI that mentions a *.zlas format for the first time. This may have been an oversight at ESRI since there was no press release, blog post, etc preceding this documentation update (that happened 11 days ago). A screenshot of the entry can be found below.

I have heard about LAS compression by ESRI since Gene mentioned it in a blog entry after ESRI’s 3D Mapping and LiDAR Forum. Back then I throught they were talking about LAZ and that our 1.5 years of talking about including support for LASzip-compressed LiDAR into ArcGIS were finally getting somewhere. But turns out they have been doing their own thing. Here some rumors I have heard about ESRI’s new *.zlas format:

  • similar compression rates as LAZ
  • includes spatial indexing
  • (maybe) re-orders points during compression
  • performance is like laszip.exe or better
  • will be available in ArcGIS 10.2.1
  • can be used without the LAS Dataset
  • “free” Windows executable will be available soon
  • development libraries with API will follow
  • ESRI has been giving data providers heads up that clients may soon demand this format

My first thought was that this might be a reengineered version of the LizardTech LiDAR CompressorTM but it is not. This seems to be ESRI’s own development. Does anyone have more details on this?

What was their motivation? Is LAZ too slow for them? I would have happily adressed whatever LASzip was lacking as (compatible) extensions to the LAZ format – which has become de-facto standard for LiDAR compression and is open source. But instead they invested serious money and man-power into creating an entirely new format. Anyone want to speculate why …?

screenshot of ArcGIS 10.2 documentation

screenshot of the ArcGIS 10.2 documentation mentioning *.zlas

UPDATE (January 7th): It is official now. Apparently, Gene who has mentioned our rumors on his blog just received heads up from Clayton Crawford at ESRI that the LAS Optimizer is available for download here. The EzLAS.zip file contains a PDF with interesting details. Below a snapshot of the GUI and what the website says:

“This executable is used to optimize and compress LAS format lidar. It creates *.zlas files, an optimized version of LAS that’s useful for archiving, sharing, and direct use. zLAS files are much smaller and more efficient to use, especially on the cloud and over networks, than regular LAS.
  •  The standalone executable does not require an ArcGIS install or license.
  •  The same executable is used for both compression and decompression.
  •  The download zip file contains more information and help in an included pdf document.”
GUI of the LAS Optimizer

GUI of the LAS Optimizer

Apparently, ArcGIS 10.2.1 is available for general release on January 7th and will support direct read of optimized LAS (*.zlas) via the LAS dataset. Now it’s your turn to try out what ESRI has cooked up and comment … (-:

UPDATE (January 9th): I was told that ESRI will be releasing an “official statement” soon explaining why they have developed their own LiDAR compression format. And they really should do so. I have received (and continue to receive) a fair number of “off-the-record” emails from people across the industry expressing feelings that range between “disappointment”, “anger”, and “disgust” over what is seen as an attempt to sabotage our multi-year effort of creating an open, free, and efficient compressed LiDAR exchange format … time for some xkcd humor.

UPDATE (January 12th): This might just be that “official statement“. Seems someone is really trying hard to avoid using the word LASzip … (-;

UPDATE (January 19th): The story has been picked up by a number of blogs like Paul Ramsey’s “LiDAR format wars“, James Fee’s “LAS, LAZ, LASzip, zLAS and You“, and Randal Hale’s “LiDAR and your software“.

UPDATE (February 7th): The front-lines harden as an unlikely coalition of open source knights, laser guardians, imperial agencies, and competing thugs forms a rebel movement against the approaching Desri Star who threatens the free world announcing that the dArcG is going to spawn “parallelized LAZ clones“. Encrypted instructions to Jedis spread like “Point Clouds on the Horizon” “Towards an Open Future” as the insurgency prepares a better future for compression, free of proprietary oppression. The clone wars might be starting soon. Will the FOSS be with LAZ? (-;

Tutorial: LiDAR preparation

Posted on by

This is part two of a three-part tutorial on how to use LAStools to implement a pipeline that (1) quality checks a newly received set of raw LiDAR flight strips, (2) tiles and prepares the LiDAR for subsequent exploitation, and (3) generates raster and vector derivatives such as DTMs, DSMs, building footprints, and iso-contours with multi-core batch processing.

To get started, download the latest version of LAStools and these 7 example flight strips and put them into a folder called ‘.\LAStools\bin\strips_raw’. For simplicity we will work inside the ‘.\LAStools\bin’ directory, so open a DOS command line window and change directory to where this folder is on your computer, for example, with the command ‘cd c:\software\LAStools\bin’ or with ‘cd d:\LAStools\bin’ followed by ‘d:’. It may be helpful (not mandatory) to follow the tutorial on quality checking first.

Create a new folder called ‘.\LAStools\bin\tiles_raw’ for storing the LiDAR tiles. Then run ‘lastile.exe’ in GUI mode and load the 7 example flight strips. You can either do this by double-clicking ‘lastile.exe’ and using the ‘browse …’ menu or by entering the command below:

C:\LAStools\bin>lastile -i strips_raw\LDR*.laz -gui

Set the GUI options as shown below to create a buffered tiling from the original flight strips. Check the button ‘files are flightlines’ so that points from different flight lines get a unique flight line ID stored in their ‘point source ID’ attribute. This makes it possible to identify later which points of a tile came from the same flight strip. Use the default 1000 to specify the tile size and request a buffer of 50 meters around every tile. This buffer helps to reduce edge artifacts at the tile boundaries in a tile-based processing pipeline.  Shift the tiling grid by 920 in x direction and by 320 in y direction so the flight strips fit in exactly 4 tiles. This is not mandatory but results in fewer tiles and therefore fewer cuts (experimentally discovered). Requests compressed output to overcome the I/O bottleneck. Set the projection to UTM zone 19 north and the vertical datum to NAVD88. When projection information is missing in the original flight strips it is a good idea to add this as early as possible so it is available in subsequent processing steps.

tutorial2 lastile GUI settings

Press ‘RUN’ and the ‘RUN’ window will pop up. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Alternatively, if you prefer to work in the command-line, you can use this almost equivalent command here. It has been enhanced by one additional parameter shown in red. It quantizes the input on-the-fly to a more appropriate centimeter [cm] resolution because airborne LiDAR does not have the millimeter [mm] resolution that the raw flight strips are stored with. See the lasinfo report shown at the end of the tutorial on quality checking to see the excessive millimeter resolution that the additional parameter in red is fixing.

C:\LAStools\bin>lastile -i strips_raw\LDR*.laz ^
                        -files_are_flightlines ^
                        -rescale 0.01 0.01 0.01 ^
                        -utm 19N -vertical_navd88 ^
                        -tile_size 1000 -buffer 50 ^
                        -tile_ll 920 320 ^
                        -odir tiles_raw -o fitch.laz

Create a new folder called ‘.\LAStools\bin\tiles_ground’ for storing the ground-classified LiDAR tiles. Then run ‘lasground.exe’ in GUI mode and load the 4 tiles. You can do this by double-clicking ‘lasground.exe’ and using the ‘browse …’ menu or by entering the command below:

C:\LAStools\bin>lasground -i tiles_raw\fitch*.laz -gui

Set the GUI options as shown below to classify the bare-earth points in all tiles. Select the ‘metro’ setting (which uses a step size of 50m) because there are very large hangars in the point cloud and the ‘ultra fine’ setting for the initial ground estimate. For more details on these parameters read the corresponding README.txt file. Specify the output directory as ’tiles_ground’ and select compressed LAZ output. Select 4 cores only if you have a computer with that many cores. This will assign tiles to different cores to process multiple tiles in parallel.

tutorial2 lasground GUI settings

Press ‘RUN’ and the ‘RUN’ window will pop up. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Alternatively, if you prefer to work in the command-line, you can use this equivalent command here:

C:\LAStools\bin>lasground -i tiles_raw\fitch*.laz ^
                          -metro -ultra_fine ^
                          -odir tiles_ground -olaz ^
                          -cores 4

Run ‘lasview.exe’ on ‘fitch_273920_4714320.laz’ in the ’tiles_ground’ folder to inspect the result of your bare-earth classification. Press ‘g’ on the keyboard to show only the ground points. After all points have loaded press ‘t’ on the keyboard to triangulate only the points shown. Press ‘-‘ to make the points disappear and press ‘h’ multiple times to iterate over the possible ways to shade the triangulation. Now press ‘=’ to make the points reappear, press ‘u’ to show the unclassified (non-ground) points, and press ‘c’ a few times to iterate over the possible ways to color the points.

tutorial2 lasview ground and cloud points

There are “dirt” points below the ground and “air” points far above the ground such as the one pointed at by the mouse cursor. We remove them with ‘lasheight.exe’. This tool computes for each point the vertical distance to the triangulation of the ground points and stores it in the ‘user_data’ field. We also need these heights for the following step of finding buildings and vegetation. Create a new folder called ‘.\LAStools\bin\tiles_height’ for storing the cleaned LiDAR tiles with height above ground information. Run ‘lasheight.exe’ in GUI mode and load the 4 ground-classified tiles by double-clicking ‘lasheight.exe’ and using the ‘browse …’ menu or by entering the command below:

C:\LAStools\bin>lasheight -i tiles_ground\fitch*.laz -gui

Configure the GUI with the settings shown below. By default ‘lasheight.exe’ uses the points classified as ground to construct a TIN and then calculates the height of all other points in respect to this ground surface. With ‘drop_above 30’ and  ‘drop_below -2′ all points that are 30 meters above or 2 meters below the ground are removed from the compressed LAZ tiles output to the ’tiles_height’ folder. Run the process on multiple cores if possible.

tutorial2 lasheight GUI settings

Press ‘RUN’. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Or use this in the DOS command-line:

C:\LAStools\bin>lasheight -i tiles_ground\fitch*.laz ^
                          -drop_below -2 -drop_above 30 ^
                          -odir tiles_height -olaz ^
                          -cores 4

Next, create a new folder called ‘.\LAStools\bin\tiles_classified’ for storing classified LiDAR tiles with building and vegetation points. Run ‘lasclassify.exe’ in GUI mode and load the 4 height-cleaned tiles by double-clicking or by entering the command below:

C:\LAStools\bin>lasclassify -i tiles_height\fitch*.laz -gui

Set the GUI settings as shown below to identify buildings and trees. Mileage may vary on this step because automatic LiDAR classification is a hard problem. The default settings require a density of at least 2 pulses per square. Setting the step size for computing planarity (or ruggedness) to 3 meters has fewer false positives but misses some smaller buildings. For more details on tuning these parameters see the corresponding README.txt file.

tutorial2 lasclassify GUI settings

Press ‘RUN’. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Or use this command-line:

C:\LAStools\bin>lasclassify -i tiles_height\fitch*.laz ^
                            -step 3 ^
                            -odir tiles_classified -olaz ^
                            -cores 4

Run ‘lasview.exe’ on ‘fitch_273920_4714320.laz’ in the ’tiles_classified’ folder and be amazed about the result of your building and vegetation classification. Although not 100 percent correct, ‘lasclassify.exe’ identifies all man-made structures as class 6 and marks all vegetation above two meters as class 5.

tutorial2 lasview classified points

Create the last folder called ‘.\LAStools\bin\tiles_final’ for storing the final LiDAR tiles without the 50 meter buffers that are ready to be shipped to the customer. Run ‘lastile.exe’ in GUI mode and load the 4 classified tiles with:

C:\LAStools\bin>lastile -i tiles_classified\fitch*.laz -gui

Set the GUI settings as shown below to remove the buffers that were added in the very beginning and carried through all processing steps to avoid artifacts along the boundaries.

tutorial2 lastile GUI settings to remove buffer

Press ‘RUN’. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Or use this command-line that has one additional parameter shown in red. This sets the user data field of each point back to zero that was set by ‘lasheight.exe’ to the height above ground (in decimeters). This number is meaningless to the customer and setting it to zero improves compression.

C:\LAStools\bin>lastile -i tiles_classified\fitch*.laz ^
                        -set_user_data 0 ^
                        -remove_buffer ^
                        -odir tiles_final -olaz

Finally let us run ‘lasinfo’ on one of the generated tiles and also compute the density:

C:\LAStools\bin>lasinfo -i tiles_final\fitch_271920_4714320.laz -cd
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            0
  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000
  version major.minor:        1.0
  system identifier:          'LAStools (c) by Martin Isenburg'
  generating software:        'lastile (131011) commercial'
  file creation day/year:     0/0
  header size:                227
  offset to point data:       5689
  number var. length records: 4
  point data format:          0
  point data record length:   20
  number of point records:    2571486
  number of points by return: 2196212 0 375274 0 0
  scale factor x y z:         0.01 0.01 0.01
  offset x y z:               0 4000000 0
  min x y z:                  272044.80 4714466.57 82.34
  max x y z:                  272919.99 4715243.91 169.21
variable length header record 1 of 4:
  reserved             43707
  user ID              'LeicaGeo'
  record ID            1001
  length after header  5120
  description          ''
variable length header record 2 of 4:
  reserved             43707
  user ID              'LeicaGeo'
  record ID            1002
  length after header  22
  description          'MissionInfo'
variable length header record 3 of 4:
  reserved             43707
  user ID              'LeicaGeo'
  record ID            1003
  length after header  54
  description          'UserInputs'
variable length header record 4 of 4:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34735
  length after header  48
  description          'by LAStools of Martin Isenburg'
    GeoKeyDirectoryTag version 1.1.0 number of keys 5
      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
      key 3072 tiff_tag_location 0 count 1 value_offset 32619 - ProjectedCSTypeGeoKey: PCS_WGS84_UTM_zone_19N
      key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
      key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter
      key 4096 tiff_tag_location 0 count 1 value_offset 5103 - VerticalCSTypeGeoKey: VertCS_North_American_Vertical_Datum_1988
the header is followed by 2 user-defined bytes
LASzip compression (version 2.2r0 c2 50000): POINT10 2
LAStiling (idx 8, lvl 2, sub 0, bbox 271920 4.71232e+006 275920 4.71632e+006)
reporting minimum and maximum for all LAS point record entries ...
  X   27204480   27291999
  Y   71446657   71524391
  Z       8234      16921
  intensity 0 255
  edge_of_flight_line 0 1
  scan_direction_flag 0 1
  number_of_returns_of_given_pulse 1 2
  return_number                    1 3
  classification      1     6
  scan_angle_rank   -13    21
  user_data           0     0
  point_source_ID     1     7
number of last returns: 2196301
covered area in square meters/kilometers: 408700/0.41
point density: all returns 6.29 last only 5.37 (per square meter)
      spacing: all returns 0.40 last only 0.43 (in meters)
overview over number of returns of given pulse: 1821027 750459 0 0 0 0 0
histogram of classification of points:
          286793  Unclassified (1)
          589019  Ground (2)
         1586840  High Vegetation (5)
          108834  Building (6)