Removal of Cloud Returns With a Coarse DTM

Flying LiDAR in regions with frequent cloud cover presents a significant challenge. If flight plan constraints do not allow to stay below all of the clouds then some of them will be scanned from above. For denser clouds this often means that all of the laser’s energy gets reflected or absorbed by the cloud and no returns on the terrain are generated. Clouds of points that are all cloud points can spell trouble in subsequent processing steps like automated classification with lasground. This is especially true for large dense clouds that start to look like features.

lasheight_masbate_lidar_clouds

scanned clouds in four neighboring LiDAR flight lines

Airborne LiDAR surveys are often carried out to create an improved Digital Terrain Model (DTM) with higher resolution than previous elevation products. Hence, usually there is already some lower resolution model that – as we will see in the following – can be used to robustly remove or mark all cloud points, at least those sufficiently high above this older ground approximation. We use data from the DREAM LiDAR Project in the Philippines who often acquire LiDAR in areas with a lot of cloud cover.

Our input are the 4 very short LiDAR strips in LAZ format shown above with lots of cloud returns and a coarse Aster DTM in ASC format at 50 m resolution. We will classify all LiDAR points far above the Aster DTM because they correspond to cloud returns. We need to be very conservative (because the Aster DTM is coarse and inaccurate) and only remove points that are really far above the Aster DTM whose ASC header is shown below.

ncols         2109
nrows         2162
xllcorner     512546.427859
yllcorner     1290961.682335
cellsize      50
NODATA_value  -9999
-9999 -9999 -9999 -9999 -9999 [...]

First we convert the Aster DTM from the inefficient ASC format to the efficient LAZ format using lassort. Why lassort? Because that puts the rasters that are on-the-fly converted to points on a grid into a spatially coherent order. This will allow efficient area-of-interest (AOI) queries once we have LAXxed the file with lasindex.

>> lassort -i masbate.asc -o masbate.laz
>> dir masbate*
    10,993,886 masbate.asc
     2,120,442 masbate.laz

Next we create a spatial indexing for ‘masbate.laz’ with lasindex. The default granularity of lasindex is 100 by 100 meter because it is set for airborne LiDAR. Given that there is only one point every 50 by 50 meters in the point cloud grid of the Aster DTM we increase the granularity to 1000 by 1000 meters.

>> lasindex -i masbate.laz -tile_size 1000 -append
>> dir masbate*
    10,993,886 masbate.asc
     2,132,006 masbate.laz

The ‘masbate.laz’ file is a little bit bigger than before because the tiny ‘masbate.lax’ file (that can also be stored separately) was appended to the end of ‘masbate.laz’ via the ‘-append’ option. Spatial indexing is realized with an underlying quadtree that can be visualized with lasview by pressing ‘Q’ (or with the pop-up menu) and the points corresponding to each quadtree cell can be turned blue by hovering above a cell and pressing ‘q’.

The points of the 50m Aster DTM as a LAZ file with the LAX spatial indexing quadtree  generated by lasindex.

Points of the 50m Aster DTM as a LAZ file and its spatial indexing quadtree generated by lasindex.

Before we use the sorted points from the Aster DTM to classify, flag, or delete the LiDAR returns far above the ground with lasheight we do a visual sanity check with lasview using the GUI settings shown below.

>> lasview -i raw_strips\*.laz ^
           -i masbate.laz ^
           -gui
GUI settings in lasview for picking a small area

GUI settings in lasview for picking a small area

By triangulating only the Aster DTM points we can visually confirm that the LiDAR points on the terrain are also close to the Aster DTM whereas the cloud returns are far up with a clear seperation between. The pink “drop lines” of length 50 meters that are dangling off each point give us a sense of scale (enabled via the pop-up menu that appears with a right click). Note how many last returns get stuck in the clouds. Those would likely cause troubles in subsequent ground classification with lasground.

Finally we run lasheight on 4 cores to classify all points that are 150 meter or more above the Aster DTM as noise (class 7).

>> lasheight -i raw_strips\*.laz ^
             -ground_points masbate.laz ^
             -classify_above 150 7 ^
             -odix _cloud -olaz ^
             -cores 4

Below the result with the cloud returns in pink (i.e. the points classified as 7). This process scales well even for larger ground points files, because lasheight uses an area-of-index (AOI) query to load only those ground points from the spatially indexed ‘masbate.laz’ file, that fall into the (slightly enlarged) bounding box of the raw LiDAR strip that is being processed.

lasheight_masbate_lidar_aster_classified_clouds

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.

Tutorial: editing LAS or LAZ files “by hand” with lasview

This tutorial describes how to manually edit LiDAR using the new inspection and editing functionality available in ‘lasview.exe’ with the latest release of LAStools (version 140301). We will work with the familiar ‘fusa.laz’ sample LiDAR data set from the LAStools distribution that was recently reported to have shown strange symptoms assumed to be side-effects of the LAZ cloning experiments in the ESRI labs … (-;

Inspecting LiDAR files with cross sections

Copy ‘fusa.laz’ from the folder ‘lastools\data’ to the folder ‘lastools\bin’. Run ‘lasview.exe’ so that it loads ‘fusa.laz’. Either do this via the GUI by double-clicking ‘lasview.exe’, loading ‘fusa.laz’ via the ‘browse …’ menu, and then clicking the ‘VIEW’ button or by entering the command below:

C:\lastools\bin>lasview -i fusa.laz

Press the <x> key to toggle to “select cross” where you can pick a rectangular “cross” section. The default cross section is a profile extending across the bottom of the bounding box.

tutorial4_lasview_01_select_cross

By pressing the <x> key again you toggle back to actually view the cross section. Holding down <ALT> you can rotate the view to look at the cross section from the side. Holding down <CTRL> you can zoom in and out. Holding down <SHIFT> you can translate up and down or left and right. Increase or decrease the size of the points pressing <=> or <->. Hover with the mouse over a point and press <i> to inspect its coordinates and attributes.

tutorial4_lasview_02_cross_inspect_point

Traverse the LiDAR file visually by moving the cross section with the arrow keys <UP> <DOWN> <LEFT> and <RIGHT>. You can move either in the “select cross” view and see the picked rectangle move or in the “cross” view and “walk” through the LiDAR. Hold down the <SHIFT> key simultaneously to take bigger steps or the <ALT> key to take smaller steps. Inspect other points by hovering over them with the cursor and pressing <i>. The point information disappears when pressing <i> with the cursor over the background.

tutorial4_lasview_03_traverse_with_arrow_keysToggle back to the “select cross” view with <x> and pick approximately the same rectangle as shown below:

tutorial4_lasview_04_pick_missclassified_roof

Changing Classifications and Deleting Points

Continuing the steps above, toggle back to the “cross” view by pressing <x>. Note that part of the roof of the house has been miss-classified as vegetation while others are left unclassified. Press <e> to turn on the “EDIT” mode and right-click to select “reclassify points as building (6)” via the pop-up menu.

tutorial4_lasview_05_reclassify_menu

Now use the cross-hair cursor to draw a polygonal fence around all points that should be reclassified. Press <ESC> to remove the last vertex of the polygon if you miss-placed it by a mistake.

tutorial4_lasview_06_reclassify_fenceOnce you are happy with your polygon press <r> to register the edit. A note appears informing you how many points had their classification changed. In the top right corner an “undo” counter appears informing you how many changes you can undo by pressing <CTRL-u>. Try it. Immediataly the changes disappear and a “redo” counter appears instead. Press <CTRL-o> to redo the change you have just undone.

tutorial4_lasview_07_edit_result

Press <CTRL-s> to save this edit as a tiny LAY file using the recently introduced LASlayers concept. In case there was already an existing LAY file (that was not applied with ‘-ilay’ when starting lasview) you will be warned and have to press <CTRL-f> to force overwriting it as shown below.

tutorial4_lasview_07a_first_save

Press the key sequence <SHIFT-b>, <t>, and <a> to get the same visuals above.

tutorial4_lasview_07a_second_save

Press <SHIFT-t> to remove the triangulation again. After saving an edit it can no longer be undone via <CTRL-u>. Instead you will have to strip off this particular layer with the layer management available through “laslayers.exe” as described here. Now press <x> to toggle to the “cross select” view.

tutorial4_lasview_08_move_select_cross

Use the <DOWN> arrow to move the selected cross section to the area shown above that has a few unclassified points in the middle of the roof. Press <x> to go back to the “cross” view and try to understand why these points are not part of the roof. Looks like they are from the top of a chimney, and antenna, or a satellite dish as they do not fit the otherwise planar roof.

Assume we need to remove them for some reason. Pan, translate, and zoom the view such that these points can be easily surrounded by a polygon. Now press <d> to enter the “DELETE” mode, fence in these points, and press <r> to register the deletion.

tutorial4_lasview_09_delete_points

It can be tricky to place a clearly seperating polygon and you may be worried about deleting a few orange building points as well. Press <u> to only display the unclassified points before pressing <r> to register the deletion.

tutorial4_lasview_10_delete_points_unclassified_only

Press <a> to see all points again, then delete the other two points by finding a good view point, pressing <d>, and drawing a polygon.

tutorial4_lasview_11_delete_other_points

After registering this deletion of two points your “undo” counter should be at two. Press <CTRL-u> twice to undo this and the last deletion, then press <CTRL-o> twice to redo them both.

tutorial4_lasview_12_delete_other_points_undo_count_2

Now press <CTRL-s> to save this deletion as another layer. It will be appended to the LAY file that already contains one layer with the roof re-classification edit we did first. Press <t> to triangulate the points in the “cross” view. See how nicly flat the triangulated roofs are now that we deleted these 6 chimney points.

tutorial4_lasview_13_second_save

Look at the size of the tiny LAY file called ‘fusa.lay’ that is in the same folder as the ‘fusa.laz’ file. It contains all the edits we have done so far and mine is only 681 bytes in size. The original LAZ file has not changed. Maybe this is all you want for now. You could send only this tiny LAY file to a colleague elsewhere and he or she could apply those changes locally when needed using the ‘-ilay’ switches. For more on this see the LASlayers page.

C:\lastools\bin>lasview -i fusa.laz
C:\lastools\bin>lasview -i fusa.laz -ilay 1
C:\lastools\bin>lasview -i fusa.laz -ilay 2

However, you may want to eventually apply the changes and produce a new LAZ file. This will be a lot slower as it requires rewriting the entire file. It will also make changes permanent. Press <CTRL-a> and a new file is produced called ‘fusa_1.laz’ that has 6 points less than ‘fusa.laz’ and 69 points with a different classification as “building”. One more thing, press <CTRL-x> if you want to toggle between the “cross” section view and the default view.

tutorial4_lasview_14_applied_laslayers You need to have a license to LAStools to save edits for file that contain 1 million points or more.

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.

Call for Input on Compression of LAS 1.4

PRESS RELEASE
January 21, 2014
rapidlasso GmbH, Gilching, Germany

The creators of the widely-used open source LiDAR compressor LASzip are issuing a “Call for Input” for extending the popular LAZ format to the new point types 6 to 10 that were introduced with the LAS 1.4 specification. Their LiDAR technology company rapidlasso GmbH invites all interested parties to “The LAS room” for an open discussion: “As LAZ has gained significant traction in the community we want to involve stakeholders across industry, academia, and government,” so Dr. Martin Isenburg, founder of rapidlasso GmbH.

Recent conversations with large agencies suggest that LAS 1.4 is not yet being tendered and that there are no significant LAS 1.4 holdings. “So we have some time,” says Dr. Isenburg, “No need to rush and put out a half-baked solution. I have the feeling LAS 1.4 will stay with us for a while, so we better get it right.”

The LASzip compression scheme will not change for LAS 1.X files that only contain point types 0 to 5 to guarantee “forward-compatibility” with older, already existing LAZ readers. “The new point types,” so Dr. Isenburg, “provide a unique opportunity to enhance the compression scheme further.”

Apart from the improvements mentioned in the LASzip journal paper, one change worth discussing is whether to encode GPS time stamps, RGBI colors, and WavePackets separately from other point attributes to allow selective reading of compressed files. Plain visualization, for example, does not need GPS time information. Another consideration is (optional) spatial indexing using LASindex. The new concept of “Appended VLRs” in LAS 1.4 allows appending the tiny LAX files of LASindex to the end of a LAZ file.

Dr. Isenburg asks especially those who do not yet use LAZ (for technical reasons) to join “The LAS room” and contribute suggestions: “Let’s make LASzip with LAS 1.4 support not only a better LAZ, but also your LAZ.”

LAX files contain a quadtree to describe which parts of a file have to be loaded to get all points within an area of interest

The LAX files of LASindex contain a quadtree that describes which parts of the file need to be loaded to get all points within an area of interest.

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? (-;

LASindex – spatial indexing of LiDAR data

Salzburg is a beautiful city in December. The European LiDAR Mapping Forum coincided with the days when the “Krampus” (= “Christmas monsters”) are roaming the Christmas markets in the old town to scare children and adults alike. One gave me a painful whipping in the legs with its leathery tail when I tried to protect a LAStools user … (-;

More to the point, here is my talk at ELMF 2012 on “LASindex – simple spatial indexing of LiDAR data”. I first give a little update on LASzip, then talk about spatial indexing with LAX, before sneak-previewing PulseWaves – our new and open LiDAR format for storing full waveform data.

Some more detail:
Airborne LiDAR surveys collect large amounts of elevations samples, often resulting in Terabytes of data. The acquired LiDAR points are typically stored and distributed in the LAS format or – its lossless compressed twin – the LAZ format. However, managing a folder of LAS or LAZ files is not a trivial task when a survey consists, for example, of 500 flight strips containing around 200 million points each. Even a simple area-of-interest (AOI) query requires opening all files and loading all those whose bounding box overlaps the queried AOI. One solution is to copy the survey into a dedicated data base such as Oracle Spatial or PostgreSQL. We present a much simpler alternative that works directly on the original LAS or LAZ files.

Our minimal-effort spatial indexing scheme has very small setup costs, avoids creating a second copy of the data, and is already in use in the LAStools software suite. For each LiDAR file we generate a tiny LAX file that resides in the same folder as the *.las or *.laz file and has the same name but with a *.lax extension. The LAX files are generally as small as 0.01 percent (for a LAS file) or 0.1 percent (for a LAZ file) of the file containing the LiDAR data and they can be generated as fast as the points can be read off disk.

The LAX files describe an adaptive quadtree over the x and y coordinates of all points. Each occupied quadtree cell stores a list of point index intervals that together reference all points falling into this cell. By merging all intervals of a cell that are less than 1000 apart in point index space we significantly reduce the number of intervals, the size of the LAX files, and the number of file seek operations.

Although individual cells typically reference too many points this is usually amortized as a typical AOI query will require returning a union of all intervals from many quadtree cells. However, our in-place spatial indexing relies on a certain degree of spatial coherency to be present in the point order. A simple measure of the efficiency of the existing order is to calculate the overhead factor when loading each quadtree cell individually from disk.

The source code for LASindex is part of the open source library LASlib of LAStools. It has been extensively field-tested in the LiDAR delivery pipeline of Open Topography (OT) where it is used to efficiently gather data from folders of LAZ files in accordance to area-of-interest queries that are generated by users via OT’s popular web-based LiDAR download interface. Another important use is on-the-fly point buffering. When batch processing, for example, 2km by 2km LiDAR tiles to create DTMs via rasterization of a temporary TIN, it is beneficial to load a 100 meter point buffer around each tile to avoid tile boundary artifacts. The presence of LAX files allows doing so efficiently on-the-fly.