RIEGL Becomes LASzip Sponsor for LAS 1.4 Extension

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 !!!

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 bty 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 … but that is not yet implemented. (-:

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

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

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

PRESS RELEASE
February 14, 2014
rapidlasso GmbH, Gilching, Germany

As a sweet foretaste to ILMF 2014, the creators of LAStoolsLASzip, 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

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

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: strip1, strip2, strip3, strip4, strip5, strip6, strip7 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)