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

 

Tutorial: quality checking

This is part one 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 the 7 example flight strips. For simplicity we will work inside the ‘.\LAStools\bin’ directory. Therefor move the folder ‘strips_raw’ that the ZIP file contains such that it its full and final path is ‘.\LAStools\bin\strips_raw’. Then 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:’.

I usually first run ‘lasview’ on newly received LiDAR to see if the data makes any sense. The command below will on-the-fly down-sample the data and display only around 10 million points:

C:\LAStools\bin>lasview -i strips_raw\LDR*.laz -points 10000000

I usually inspect a few area-of-interests (AOIs) to get a few close-up looks at the data. If we do this several times for different areas it is best to first create a spatial indexing of the strips with ‘lasindex’ which it significantly speeds-up our spatial queries.

C:\LAStools\bin>lasindex -i strips_raw\LDR*.laz -cores 3

Only add the ‘-cores 3’ option if you have a computer with at least 4 cores. If you have an 8 core machine you may even add ‘-cores 7’ and so on. Now start ‘lasview’ in the GUI mode either by double-clicking the executable and using the ‘browse …’ roll-out on the left side to load the seven flight strips or simply by adding ‘-gui’ to the command line.

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

Now you can maneuver up and down the file list to highlight the different strips (close the ‘browse …’ roll-out to see the header contents listed on the bottom left). You can again look at all the data by pressing the ‘VIEW’ button or you can inspect a smaller area by first picking a rectangular region as shown in the screenshot below.

lasview_pick_aoi_in_GUIIn the GUI above you can see that the files have no projection information. We will add this later. The x, y, and z scaling factors are set to 0.001 which means that points are stored with millimeter resolution. Because airborne LiDAR is far from being that accurate we will later change the scaling factor to a more appropriate centimeter resolution. Finally, the creation day and year is missing, which – if we know when the data was flown – will be an easy fix.

After starting ‘lasview’ you can use the ‘c’ hot-key to toggle through the different ways to color the points. Another way to change the point coloring is by presetting it in the GUI as shown above or via the ‘color by …’ menu that opens when you right click.

lasview visualize flight-lines (aka point_source_ID)The coloring by flight-line which is supposed to be stored in the ‘point_source_ID’ field of each point seems rather odd. This field was probably used for some other purpose and not cleaned up before the data was distributed. Not nice and we will fix that later. There are a few other visualization shown below. Coloring points by return turns single returns yellow, first of many returns red, and last of many returns blue. Usually there are also a few green points corresponding to intermediate (neither first nor last) of many returns but there are none in this data set as the ‘lasinfo’ report will confirm later. Toggle through the colorings with the ‘c’ key until it displays the intensities. You may need to press the ‘=’ key a few times to increase the point size. Next, press the ‘t’ key to triangulate the points and then decrease the point size with the ‘-‘ key until they have disappeared and you can see a hillside shading of the TIN. Use the ‘h’ key to toggle through different shadings for the TIN triangles.

Before attempting to improve these LAZ files and prepare them for processing we want to make sure that our work will be worthwhile by running a quick visualization based of how well the flight strips fit together. We do this with the ‘lasoverlap’ tool.

C:\LAStools\bin>lasoverlap -i strips_raw\LDR*.laz ^
                           -files_are_flightlines ^
                           -utm 19N ^
                           -step 1.0 -max_diff 1.0 ^
                           -o strip_overlap.png

I happen to know that the missing projection in the LAZ files is the UTM zone 19N. I can simply add this to the command line with ‘-utm 19N’. Then ‘lasoverlap’ will produce in addition to the PNG output also a KML file that allows visualizing the overlap and difference rasters within the geo-spatial context of Google Earth.

Serious troubles in the data would be evidenced by void areas in the overlap raster that are obviously caused by poor flight planning (not water bodies) and by deeply blue or deeply red saturated areas in open (=> non-forested) terrain in the difference raster. If either was the case you should send the data back to the vendor to have him fix it … (-:

lasoverlap_strip_overlap_overlasoverlap_strip_overlap_diffNow we compare the LiDAR elevations with a set of 29 exactly known control points that were obtained through a ground survey that are listed below. If you want to follow this exercise along you should make sure that you already have the file ‘strips_raw\cps.csv’. Or simply copy the text below and store it to a file called ‘strips_raw\cps.csv’.

Type,Easting,Northing,Z
Open/Paved,273299.68,4715133.88,74.17
Open/Paved,273477.61,4714979.29,73.85
Open/Paved,274001.2,4714540.29,72.5
Open/Paved,273670.13,4714817.4,73.34
Open/Paved,273677.42,4715018.66,74.26
Open/Paved,273400.1,4714528.98,73.15
Open/Paved,274511.59,4714905.97,95.7
Open/Paved,275074.66,4714841.98,120.18
Open/Paved,275409.65,4714994.76,113.41
Field,273321.18,4714946.83,73.46
Field,273601.49,4715101.78,74.35
Field,273646.76,4714972.94,73.97
Field,273890.5,4714457.59,71.13
Field,274650.24,4714903.44,105.13
Field,274522.36,4714829.74,98.59
Field,275474.47,4714780.03,127.63
Field,275636.39,4714868.85,120.72
Field,274747.37,4714932.57,116.11
Field,272795.36,4714503.86,126.9
Forested,272547.02,4714623.09,127.71
Forested,273205.33,4714900.27,79.36
Forested,272530.52,4715045.46,113.48
Forested,275237.48,4715049.57,120.31
Forested,275268.37,4714543.82,104.99
Forested,274666.09,4714497.49,108.86
Forested,274403.56,4715053.43,93.17
Forested,274901.6,4714493.63,114.64
Forested,274658.37,4715072.74,104.49
Forested,274121.73,4714524.51,72.06

The first entry describes the location of the control point and the following three entries are its x, y, and z coordinate. We now compare the seven LiDAR strips against these 29 control points. You can also do thus via the GUI as shown below but you need to make sure to that the ‘keep ground (2) ‘ and the ‘keep keypoint (8)’ check-buttons are not checked since our LiDAR data is not yet classified.

lascontrol_strips_GUI

C:\LAStools\bin>lascontrol -i strips_raw\*.laz ^
                           -cp strips_raw\cps.csv ^
                           -parse sxyz
 diff,lidar_z,Type,Easting,Northing,Z
 0.1451,74.3151,Open/Paved,273299.68,4715133.88,74.17
 0.0190857,73.8691,Open/Paved,273477.61,4714979.29,73.85
 0.08872,72.5887,Open/Paved,274001.2,4714540.29,72.5
 -0.00419681,73.3358,Open/Paved,273670.13,4714817.4,73.34
 0.054868,74.3149,Open/Paved,273677.42,4715018.66,74.26
 0.0318439,73.1818,Open/Paved,273400.1,4714528.98,73.15
 -0.0720668,95.6279,Open/Paved,274511.59,4714905.97,95.7
 -0.0670072,120.113,Open/Paved,275074.66,4714841.98,120.18
 0.0756029,113.486,Open/Paved,275409.65,4714994.76,113.41
 0.00132683,73.4613,Field,273321.18,4714946.83,73.46
 0.0750162,74.425,Field,273601.49,4715101.78,74.35
 0.00228472,73.9723,Field,273646.76,4714972.94,73.97
 0.166194,71.2962,Field,273890.5,4714457.59,71.13
 -0.392266,104.738,Field,274650.24,4714903.44,105.13
 -0.0705893,98.5194,Field,274522.36,4714829.74,98.59
 -0.423622,127.206,Field,275474.47,4714780.03,127.63
 0.217528,120.938,Field,275636.39,4714868.85,120.72
 -0.12489,115.985,Field,274747.37,4714932.57,116.11
 0.0195411,126.92,Field,272795.36,4714503.86,126.9
 17.0342,144.744,Forested,272547.02,4714623.09,127.71
 9.20693,88.5669,Forested,273205.33,4714900.27,79.36
 26.355,139.835,Forested,272530.52,4715045.46,113.48
 20.2143,140.524,Forested,275237.48,4715049.57,120.31
 14.5363,119.526,Forested,275268.37,4714543.82,104.99
 0.173861,109.034,Forested,274666.09,4714497.49,108.86
 20.5742,113.744,Forested,274403.56,4715053.43,93.17
 0.00432622,114.644,Forested,274901.6,4714493.63,114.64
 23.8816,128.372,Forested,274658.37,4715072.74,104.49
 0.3645,72.4245,Forested,274121.73,4714524.51,72.06
 sampled TIN at 29 of 29 control points.
 avgabs/rms/stddev/avg of elevation errors are
 4.63438/9.61987/8.62324/4.55475 meter. skew is 1.47593.

The output of ‘lascontrol’ could be directed into a file using the ‘-o control.txt’ option. The tool prefixes two numbers to the beginning of every line: the elevation value that was computed from the LiDAR data at the xy location of the control point and the difference between the computed elevation and the elevation of the control point. At the end a summary reports the average absolute error, the relative mean-square error, the standard deviation, and the average error across all the differences. Of course, for control points of type ‘Forested’ we have terrible results because we use all LiDAR points in this computation and not just those that correspond to ground returns. Therefore we will have to repeat this exercise at a later point once we have ground-classified the LiDAR. Via the GUI of ‘lascontrol’ you can zoom in on different areas-of-interests and visually inspect how well the data matches the control points (see below).

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

C:\LAStools\bin>lasinfo -i strips_raw\LDR030828_211804_0.laz -cd
reporting all LAS header entries:
 file signature:             'LASF'
 file source ID:             0
 global_encoding:            0
 project ID GUID data 1-4:   0 0 0 ''
 version major.minor:        1.0
 system identifier:          'ALSXX'
 generating software:        'ALSXX_PP V2.56c'
 file creation day/year:     0/0
 header size:                227
 offset to point data:       5587
 number var. length records: 3
 point data format:          0
 point data record length:   20
 number of point records:    2404613
 number of points by return: 2210130 0 194483 0 0
 scale factor x y z:         0.001 0.001 0.001
 offset x y z:               0 4000000 0
 min x y z:                  272254.830 4714389.375 65.050
 max x y z:                  275391.197 4714711.445 169.208
 variable length header record 1 of 3:
 reserved             43707
 user ID              'LeicaGeo'
 record ID            1001
 length after header  5120
 description          ''
 variable length header record 2 of 3:
 reserved             43707
 user ID              'LeicaGeo'
 record ID            1002
 length after header  22
 description          'MissionInfo'
 variable length header record 3 of 3:
 reserved             43707
 user ID              'LeicaGeo'
 record ID            1003
 length after header  54
 description          'UserInputs'
 the header is followed by 2 user-defined bytes
 LASzip compression (version 2.0r2 c2 50000): POINT10 2
 reporting minimum and maximum for all LAS point record entries ...
 X  272254812  275391197
 Y  714389375  714711445
 Z      65050     169208
 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     1
 scan_angle_rank   -13    13
 user_data           0     0
 point_source_ID     0   255
 WARNING: 2 points outside of header bounding box
 number of last returns: 2209098
 covered area in square units/kilounits: 728996/0.73
 point density: all returns 3.30 last only 3.03 (per square units)
 spacing: all returns 0.55 last only 0.57 (in units)
 overview over number of returns of given pulse: 2014615 389998 0 0 0 0 0
 histogram of classification of points:
 2404613  Unclassified (1)
 real max x larger than header max x by 0.000422
 real min x smaller than header min x by 0.017715
Any ‘WARNING’ messages in the output of ‘lasinfo’ spells potential trouble. In addition to the things already mentioned (e.g. wrong scaling factor, no creation date, missing projection VLR, random point_source_IDs) we also notice three VLRs that mean nothing to us as well as a slightly inaccurate bounding box. Another strange thing is that although all returns are reported to be from pulses with maximally 2 returns there are many returns with number 3. At the same time all return numbers are either 1 or 3 and there is no return with number 2.  Was this an older scanner with maximally two return per pulse? This does not seem LAS specification conform. In the next part of this tutorial we will turn these imperfect strips into a beautiful tiling.

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.