First Open LiDAR in Germany

UPDATE: (January 6th): Our new tutorial “downloading Bonn in LiDAR“.
UPDATE: (January 9th): Now a second state went open LiDAR as well.
UPDATE: (March 6th): New tutorial merging Bonn into proper LAS files

Kudos to OpenNRW for offering online download links for hundreds of Gigabytes of open LiDAR for the entire state of North Rhine-Westphalia (Nordrhein-Westfalen) as announced a few months ago:

More and more countries, states, and municipalities are deciding to make their LiDAR archives accessible to the general public. Some are doing entirely for free with instant online download and a generous open license that allows data sharing and commercial use. Others still charge a “small administrative fee” and require filling our actual paperwork with real signatures in ink and postal mailing of hard drives that can easily take half a year to complete. Some licences are also stricter in terms of data sharing and commercial use.

And then there was Germany where all the LiDAR data has traditionally been locked up in some cave by the 16 state survey departments and was sold for more than just a fee. Financial reasons would usually prohibit residents in Germany from making, for example, an elevation profile for their favorite mountain bike trail for hobby purposes. Or starting a small side business that (for 5 Euro a sheet) sells “Gassi Maps” with elevation-optimized dog walking paths of low incline that go past suitable potty spots and dog-friendly coffee shops.

The reason that many national and state mapping agencies have opened their LiDAR holdings for free and unencumbered access are manifold. A previous blog post had looked at the situation in England whose Environment Agency also used to sell LiDAR data and derivatives. The common argument that government agencies have been using to justify selling data (paid for with taxes) to those very same tax payers was that this would be used to finance future surveys.

It was the “freedom of information” request by Louise Huby on November 21st, 2014 that exposed this as a flawed argument. The total amount of revenue generated for all LiDAR and derivatives sales by the Environment Agency was just around £323,000 per year between 2007 and 2014. This figure was dwarfed by its annual operational budget of £1,025,000,000 in 2007/08. The revenue from LiDAR sales was merely 0.03 percent of the agencies’s budget. That can maybe pay for free coffee and tea in the office, but not for future airborne LiDAR flights. The reaction was swift. In September 2015 the Environment Agency made all their DTM and DSM rasters down to 0.25 meter resolution available online for open access and in March 2016 added the raw point clouds for download in LAZ format with a very permissible license. It has since been an incredible success story and the Environment Agency has been propelled into the role of a “champion for open data” as sketched in my ACRS 2016 keynote talk that is available on video here.

Frustrated with the situation in Germany and inspired by the change in geospatial data policy in England, we have been putting in similar “Frag den Staat” freedom of information requests with 11 of the 16 German state mapping agencies, asking about how much sales revenue was generated annually from their LiDAR and derivative sales. These four states denied the request:

We never heard back from Lower Saxony (Niedersachsen) and Thuringia (Thüringen) wanted more fees than we were willing to spend on this, but our information requests were answered by five states that are here listed with the average amount of reported revenue per year in EUR:

LiDAR acquisitions are expensive and while it would be interesting to also find out how much each state has actually spent on airborne surveys over the past years (another “Frag den Staat” request anyone?), it is obvious that the reported revenues are just a tiny fraction of those costs. Exact details of the reported revenue per year can be accessed via the links to the information requests above. The cost table of each answer letter is copied below.

However, it’s not all bad news in Germany. Some of you may have seen my happy announcement of the OpenNRW initiative that come January 2017 this would also include the raw LiDAR points. And did it happen? Yes it did! Although the raw LiDAR points are maybe a little tricky to find, they are available for free download and they come with a very permissible license.

The license is called “Datenlizenz Deutschland – Namensnennung – Version 2.0” or “dl-de/by-2-0” and allows data and derivative sharing as well as commercial use. It merely requires you to properly name the source. For the LiDAR you need to list the “Land NRW (2017)” with the year of the download in brackets as the source and specify the data set that was used via the respective Universal Resource Identification (URI) for the DOM and/or the DGM.

The OpenNRW portal now also offers the download of the LiDAR "punktwolke" (German for point cloud).

The OpenNRW portal now also offers the download of the LiDAR “punktwolke” (German for point cloud).

Follow this link to get to the open data download portal. Now type in “punktwolke” (German for “point cloud”) into the search field and on this January 3rd 2017 that gives me 2 “Ergebnisse” (German for “results”). The LiDAR point cloud representation is a bit unusual by international standards. One link is for the DGM (German for DTM) and the other for the DOM (German for DSM). Both links are eventually leading you to unstructured LiDAR point clouds that are describing these surfaces. But it’s still a little tricky to find them. First click on the “ATOM” links that get you to XML description with meta information and a lot of links for the DTM and the DSM. Somewhere hidden in there you find the actual download links for hundreds of Gigabytes of LiDAR for the entire state of North Rhine-Westphalia (Nordrhein-Westfalen):

We download the two smallest zipped files DGM and DOM for the municipality of Wickede (Ruhr) to have a look at the data. The point cloud is in EPSG 5555 which is a compound datum of horizontal EPSG 25832 aka ETRS89 / UTM zone 32N and vertical EPSG 5783 aka the “Deutches Haupthohennetz 1992”.

The contents of the DGM zip file contains multiple files per tile.

The contents of the DGM zip file contains multiple files per tile.

The DGM zip file has a total of 212 *.xyz files that list the x, y, and z coordinate for each point in ASCII format. We first compress them and add the EPSG 25832 code with laszip. The compressed LAZ files are less than half the size of the zipped XYZ files. Each file corresponds to a particular square kilometer. The name of each tile contains the lower left coordinate of this square kilometer but there can be multiple files for each square kilometer:

  • 14 files with “ab” in the name contain very few points. They look like additional points for under bridges. The “b” is likely for “Brücke” (German for “bridge”).
  • 38 files with “ag” in the name contain seem to contain only points in areas where buildings used to cover the terrain but with ground elevation. The “g” is likely for “Gebäude” (German for “building”).
  • 30 files with “aw” in the name contain seem to contain only points in areas where there are water bodies but with ground elevation. The “w” is likely for “Wasser” (German for “water”).
  • 14 files with “brk” in the name also contain few points. They look like the original bridge point that are replaced by the points in the files with “ab” in the name to flatten the bridges. The “brk” is also likely for “Brücke” (German for “bridge”).
  • 42 files with “lpb” in the name. They look like the last return LiDAR points that were classified as ground. The “lpb” is likely for “Letzter Pulse Boden” (German for “last return ground”).
  • 42 files with “lpnb” in the name. They look like those last return LiDAR points that were classified as non-ground. The “lpnb” is likely for “Letzter Pulse Nicht Boden” (German for “last return not ground”).
  • 32 files with “lpub” in the name contain very few points. They look like the last return points that are too low and were therefore excluded. The “lpub” is likely for “Letzter Pulse Unter Boden” (German for “last return under ground”).

It is left to an exercise to the user for figure out which of those above sets of files should be used for generating a raster DTM. (-: Give us your ideas in the comments. The DOM zip file has a total of 72 *.xyz files. We also compress them and add the EPSG 25832 code with laszip. The compressed LAZ files are less than half the size of the zipped XYZ files. Again there are multiple files for each square kilometer:

  • 30 files with “aw” in the name contain seem to contain only points in areas where there are water bodies but with ground elevation. The “w” is likely for “Wasser” (German for “water”).
  • 42 files with “fp” in the name. They look like the first return LiDAR points. The “fp” is likely for “Frühester Pulse” (German for “first return”).

If you use all these points from the DOM folder your get the nice DSM shown below … albeit not a spike-free one.

A triangulated first return DSM generated mainly from the file "dom1l-fp_32421_5705_1_nw.laz" with the points from "dom1l-aw_32421_5705_1_nw.laz" for areas with water bodies shown in yellow.

A triangulated first return DSM generated mainly from the file “dom1l-fp_32421_5705_1_nw.laz” with the points from “dom1l-aw_32421_5705_1_nw.laz” for areas with water bodies shown in yellow.

Kudos to OpenNRW for being the first German state to open their LiDAR holdings. Which one of the other 15 German state survey departments will be next to promote their LiDAR as open data. If you are not the last one to do so you can expect to get featured here too … (-;

UPDATE (January 5th): The folks at OpenNRW just tweeted us information about the organization of the zipped archives in the DTM (DGM) and DSM (DOM) folders. We guessed pretty okay which points are in which file but the graphic below (also available here) summarizes it much more nicely and also tells us that “a” was for “ausgefüllt” (German for “filled up”). Maximally two returns per pulse are available: either a single return or a first return plus a last return. There are no intermediate returns, which may be an issue for those interested in vegetation mapping.

Illustration of which LiDAR point is in which file.

Nice illustration of which LiDAR point is in which file. All files with ‘ab’, ‘ag’, or ‘aw’ in the name contain synthetic points that fill up ground areas not properly reached by the laser.

Fixing Intensity Differences between Flightlines (“quick and dirty”)

Visiting our users on-site, such as last week at Mariano Marcos State University in Ilocos Norte in the Philippines, we sometimes come across situations as pictured below where the intensity values of the returns of one flightline are drastically different from that of other flightlines.

The intensity of returns in the left most flightline is different from that of other flightlines.

The intensity of returns in the left most flightline is different from that of other flightlines.

Using intensity rasters with such dark strips as an additional input for land cover classification may likely make things worse. Radiometrically “correct” intensity calibration is a complex topic and may not always be possible to do using only the LAZ files without meta information such as the internals of the scanning system and the aircraft trajectory. However, we now describe a “quick and dirty” fix to the situation shown above so that the intensity grids (that were computed as averages of first return intensities) at least “look” as sensible as for the one square tile (shown below) that was corrected by a simple multiplication with 5 for all intensities of the dark strip.

Simply multiplying all intensities of the dark flightline with 5 seems to "fix" the issue.

Simply multiplying all intensities of the dark flightline with 5 seems to “fix” the issue for our test tile.

The number 5 was determined by a quick glance at the intensity histograms that we can generate with lasinfo. We decide to only look at single returns as we expect them to have a higher correlation: Their locations are more likely to be “seen similarly” from and their energy is more likely “reflected similarly” to different flightlines compared to that of multiple returns.

lasinfo -i strip1.laz strip2.laz strip3.laz ^
        -keep_single ^
        -histo intensity 1 ^
        -nmm -nh -nv ^
        -odix _histo_int -otxt

The resulting histograms for the dark ‘strip1.laz’ is quite different from that of the much brighter ‘strip2.laz’ and ‘strip3.laz’. The average single return intensity for the dark ‘strip1.laz’ is a meager 5.13 whereas the  brighter ‘strip2.laz’ and ‘strip3.laz’ have similar averages of 24.15 and 24.50 respectively.

Draw your own conclusion about which scale factor to use. We have the choice to match either the peak of the histograms or their averages. Scaling the peak of 3 for ‘strip1.laz’ to match the 25 of the other two strips is too much of an upscaling. But the average 24.15 divided by 5.13 gives a potential scale of 4.71 and the average 24.50 divided by 5.13 gives a potential scale of 4.77 and we already saw a multiplication by 5 giving reasonable results. So this is how we can fix the intensity:

las2las -i strip1.laz ^
        -scale_intensity 4.75 ^
        -odix _corr_int -olaz

But what if your data is already in tiles? How can you adjust only the intensity of those returns that are from the flightline 1? Assuming that your flightline information is properly stored in the point source ID field of every point this is easily done with the new ‘-filtered_transform’ in LAStools using las2las on as many cores as you have as follows:

las2las -i tiles/*.laz ^
        -keep_point_source 1 ^
        -filtered_transform ^
        -scale_intensity 4.75 ^
        -odir tiles_corr -olaz ^
        -cores 8

This is not currently exposed in the GUI of las2las but you can simply add it by typing it into the ‘RUN’ pop-up window as shown below.

Scaling only the intensities of flightline 1 by 4.9 using the new '-filtered_transform'.

Scaling only the intensities of flightline 1 by 4.9 using the new ‘-filtered_transform’.

After this “quick and dirty” intensity correction we again ran lasgrid as follows:

lasgrid -i tiles_corr/*.laz ^
        -gray -set_min_max 0 60 ^
        -odir tiles_int_rasters -opng ^
        -cores 8

And the result is shown below. The obvious flightline-induced discontinuity in the intensities has pretty much disappeared. Do you have similar flightline-related intensity issues? We like to hear from you whether this technique works or if we need to implement something more clever in the future …

lasgrid_intensity_differences_3

New ‘laspublish’ creates Web Portals for 3D Viewing and Downloading of LiDAR

PRESS RELEASE (for immediate release)
February 18, 2016
rapidlasso GmbH, Gilching, Germany

Just in time for ILMF 2016, the makers of the open LASzip LiDAR compressor annouce the latest addition to their LiDAR processing software LAStools. The new ‘laspublish‘ from rapidlasso GmbH creates stand-alone Web Portals for interactive 3D viewing of LiDAR points and for selective downloading of LAZ or LAS files. The new tool is based on the cutting-edge streaming point cloud viewing technology of Potree that optimizes large LiDAR point clouds for streaming via the Web such that anyone can visualize, explore, and (optionally) download them with any modern browser.

3D viewer and download portal created with 'laspublish'

3D viewer and download portal created with ‘laspublish’

The interactive 3D viewer streams on demand only the relevant parts of the point clouds. It not only visualizes the LiDAR in many useful and intuitive ways, but is also equipped with measurement tools to calculate distances or areas and profile or clipping tools for close-up inspections. With its integrated LASzip compression, options to color LiDAR by classification, return type, intensity, and point source ID, stunning visuals via Eye Dome Lighting (EDL), and the optional 2D download map, the new ‘laspublish‘ empowers professional and novice users alike to create stand-alone LiDAR Webportals with just a few button clicks.

a few clicks and 'laspublish' creates a professional LiDAR portal

a few clicks and ‘laspublish’ creates a professional LiDAR portal

The new ‘laspublish‘ is now an integral part of the LAStools software and is bundled together with all necessary components of the Potree software. It provides an instant and cost-effective solution for generating a set of Web pages that realize a self-contained LiDAR portal offering interactive online visualization and exploration as well as easy and intuitive distribution of large LiDAR data sets via the optional download map.

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. As the only Diamond sponsor, rapidlasso GmbH has been the main financial supporter of the open source Potree package by Markus Schütz over the past two years.

About Potree:
Potree is a WebGL based viewer for large point clouds. The project evolved as a Web based viewer from the Scanopy desktop point cloud renderer by TU Wien, Institute of Computer Graphics and Algorithms. It will continue to be free and open source with a FreeBSD license to enable anyone to view, analyze and publicly share their large datasets. Visit http://potree.org for more information.

Create proper LAS 1.4 files with LAStools (for free)

So your LiDAR processing workflow produces beautiful LAS or LAZ output. You think the files are nothing short of perfect. But they are in LAS 1.2 format and the tender document explicitly requests delivery in the latest LAS 1.4 format. The free and open source las2las module of LAStools – also known as the Swiss Army knife of LiDAR processing – can easily up-convert them.

There are two possible scenarios: (1) The client wants the data as LAS 1.4 files but does not care about which point type is used. (2) The client wants LAS 1.4 and specifically asks for one of the new point types such as point type 6.

Let us assume you have LAS 1.2 content with point type 1 like these LAZ files from Martin county in Minnesota. Use these four tiles 5126-05-42.laz5126-05-43.laz5126-06-42.laz, and 5126-06-43.laz to repeat the LAS 1.4 up-conversion we are doing in the following. You will need to use version 160110 of LAStools (or newer) that you can download here.

A report generated with lasinfo shows geo-referencing information stored as GeoTIFF keys in the first VLR. Then there are 10 useless numbers stored in the second VLR that seem left-over from an earlier version of geo-referencing without EPSG codes. The two vendor-specific ‘NIIRS10’ VLRs should probably be removed unless they are meaningful to you. The LAS header is followed by 1564 user-defined bytes (sometimes used as “padding”) that we can also delete before delivery.

D:\LAStools\bin>lasinfo -i martin\5126-05-42.laz
reporting all LAS header entries:
 file signature: 'LASF'
 file source ID: 0
 global_encoding: 0
 project ID GUID data 1-4: E603549E-B74A-4696-3BAD-EA49F643C221
 version major.minor: 1.2
 system identifier: 'LAStools (c) by Martin Isenburg'
 generating software: 'lassort (110815) unlicensed'
 file creation day/year: 119/2011
 header size: 227
 offset to point data: 2163
 number var. length records: 4
 point data format: 1
 point data record length: 28
 number of point records: 13772409
 number of points by return: 13440460 240598 78743 12608 0
 scale factor x y z: 0.01 0.01 0.01
 offset x y z: 0 0 0
 min x y z: 361818.33 4855898.31 349.89
 max x y z: 364400.98 4859420.79 404.22
variable length header record 1 of 4:
 reserved 43707
 user ID 'LASF_Projection'
 record ID 34735
 length after header 40
 description 'by LAStools of Martin Isenburg'
 GeoKeyDirectoryTag version 1.1.0 number of keys 4
 key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
 key 3072 tiff_tag_location 0 count 1 value_offset 26915 - ProjectedCSTypeGeoKey: NAD83 / UTM 15N
 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
variable length header record 2 of 4:
 reserved 43707
 user ID 'LASF_Projection'
 record ID 34736
 length after header 80
 description 'GeoTiff double parameters'
 GeoDoubleParamsTag (number of doubles 10)
 500000 0 -93 0.9996 0 1 6.37814e+006 298.257 0 0.0174533
variable length header record 3 of 4:
 reserved 43707
 user ID 'NIIRS10'
 record ID 1
 length after header 26
 description 'NIIRS10 Tile Index'
variable length header record 4 of 4:
 reserved 43707
 user ID 'NIIRS10'
 record ID 4
 length after header 10
 description 'NIIRS10 Timestamp'
the header is followed by 1564 user-defined bytes
LASzip compression (version 2.1r0 c2 50000): POINT10 2 GPSTIME11 2
reporting minimum and maximum for all LAS point record entries ...
 X 36181833 36440098
 Y 485589831 485942079
 Z 34989 40422
 intensity 1 4780
 return_number 1 4
 number_of_returns 1 4
 edge_of_flight_line 0 0
 scan_direction_flag 0 1
 classification 1 10
 scan_angle_rank -24 24
 user_data 0 32
 point_source_ID 5401 8015
 gps_time 243852.663596 403514.323142
number of first returns: 13440460
number of intermediate returns: 91408
number of last returns: 13440349
number of single returns: 13199808
overview over number of returns of given pulse: 13199808 323647 198449 50505 0 0 0
histogram of classification of points:
 6560507 unclassified (1)
 6744481 ground (2)
 401464 medium vegetation (4)
 55433 building (6)
 100 noise (7)
 10157 water (9)
 267 rail (10)

(1) Create LAS 1.4 files with point type 1

The las2las command shown below turns a folder of LAS 1.2 files into a folder of  LAS 1.4 files with option ‘-set_version 1.4’ without changing the point type. It removes the last three VLRs with option ‘-remove_vlrs_from_to 1 3’ and the user-defined bytes in the LAS header with option ‘-remove_padding’ (*). The result is as LAZ with option ‘-olaz’ and the task is run on up to 4 CPUs in parallel with option ‘-cores 4’.

(*) option ‘-remove_padding’ is called  ‘-remove_extra’ in older versions of LAStools.

las2las -i martin/*.laz ^
        -remove_vlrs_from_to 1 3 ^
        -remove_padding ^
        -set_version 1.4 ^
        -odir martin_LAS14_pt1 -olaz ^
        -cores 4

(2) Create LAS 1.4 files with point type 6

The las2las command shown below turns a folder of LAS 1.2 files into a folder of  LAS 1.4 files with option ‘-set_version 1.4’ and changes the point type from 1 to 6 with option ‘-set_point_type 6’. As before some VLRs and the header padding is removed. It also adds a second description of the georeferencing information by rewriting the existing information stored as GeoTIFF tags into a properly formatted OGC WKT string with option ‘-set_ogc_wkt’. Thanks to ESRI’s assault on the LAZ format the output still needs to be in LAS. It can be compressed with laszip in a subsequent step using the ‘LAS 1.4 compatibility mode‘ sponsored by NOAA.

las2las -i martin/*.laz ^
        -remove_vlrs_from_to 1 3 ^
        -remove_padding ^
        -set_version 1.4 ^
        -set_point_type 6 ^
        -set_ogc_wkt ^
        -odir martin_LAS14_pt6 -olas ^
        -cores 4

laszip -i martin_LAS14_pt6/*.las -cores 4

Below the output of lasinfo for the “upgraded” LAS 1.4 file with the OGC WKT string. Some of the key changes have been marked in red.

D:\LAStools\bin>lasinfo -i martin_LAS14_pt6\5126-05-42.las
reporting all LAS header entries:
 file signature: 'LASF'
 file source ID: 0
 global_encoding: 16
 project ID GUID data 1-4: E603549E-B74A-4696-3BAD-EA49F643C221
 version major.minor: 1.4
 system identifier: 'LAStools (c) by rapidlasso GmbH'
 generating software: 'las2las (version 160110)'
 file creation day/year: 119/2011
 header size: 375
 offset to point data: 1135
 number var. length records: 2
 point data format: 6
 point data record length: 30
 number of point records: 0
 number of points by return: 0 0 0 0 0
 scale factor x y z: 0.01 0.01 0.01
 offset x y z: 0 0 0
 min x y z: 361818.33 4855898.31 349.89
 max x y z: 364400.98 4859420.79 404.22
 start of waveform data packet record: 0
 start of first extended variable length record: 0
 number of extended_variable length records: 0
 extended number of point records: 13772409
 extended number of points by return: 13440460 240598 78743 12608 0 0 0 0 0 0 0 0 0 0 0
variable length header record 1 of 2:
 reserved 43707
 user ID 'LASF_Projection'
 record ID 34735
 length after header 40
 description 'by LAStools of Martin Isenburg'
 GeoKeyDirectoryTag version 1.1.0 number of keys 4
 key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
 key 3072 tiff_tag_location 0 count 1 value_offset 26915 - ProjectedCSTypeGeoKey: NAD83 / UTM 15N
 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
variable length header record 2 of 2:
 reserved 43707
 user ID 'LASF_Projection'
 record ID 2112
 length after header 612
 description 'by LAStools of rapidlasso GmbH'
 WKT OGC COORDINATE SYSTEM:
 PROJCS["NAD83 / UTM 15N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-93],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26915"]]
reporting minimum and maximum for all LAS point record entries ...
 X 36181833 36440098
 Y 485589831 485942079
 Z 34989 40422
 intensity 1 4780
 return_number 1 4
 number_of_returns 1 4
 edge_of_flight_line 0 0
 scan_direction_flag 0 1
 classification 1 10
 scan_angle_rank -24 24
 user_data 0 32
 point_source_ID 5401 8015
 gps_time 243852.663596 403514.323142
 extended_return_number 1 4
 extended_number_of_returns 1 4
 extended_classification 1 10
 extended_scan_angle -4000 4000
 extended_scanner_channel 0 0
number of first returns: 13440460
number of intermediate returns: 91408
number of last returns: 13440349
number of single returns: 13199808
overview over extended number of returns of given pulse: 13199808 323647 198449 50505 0 0 0 0 0 0 0 0 0 0 0
histogram of classification of points:
 6560507 unclassified (1)
 6744481 ground (2)
 401464 medium vegetation (4)
 55433 building (6)
 100 noise (7)
 10157 water (9)
 267 rail (10)

Potree puts Big and Beautiful LiDAR in Your Browser

PRESS RELEASE (for immediate release)
September 14, 2015
rapidlasso GmbH, Gilching, Germany

Just in time for INTERGEO 2015, the Potree software was released in its latest 1.3 version. Potree is a WebGL based point cloud viewer for very large datasets. The Potree software allows to publish large LiDAR point clouds on the Web such that anyone can explore the data with nothing more but a modern browser. The interactive 3D viewer not only visualizes the LiDAR in many useful and intuitive ways but also comes with tools to perform various measurements. As its only Gold Sponsor, rapidlasso GmbH is the main supporter of this powerful open source package by Markus Schütz.

17.7 billion points around San Simeon,, CA courtesy of Open Topography

17.7 billion points from San Simeon, CA courtesy of Open Topography

The long-term sponsorship of rapidlasso GmbH has directly supported a number of useful features such as the integration of our award-winning LASzip compressor using the pure javascript version contributed by Hobu Inc, optimization for massive airborne LiDAR data, profile selection, tools for distance and area measurements, options to color by classification, return type, and point source ID, and a clipping tool. The particular features sponsored in the most recent 1.3 release of Potree are the incredible Eye Dome Lighting (EDL) and faster data conversion for large data sets. A number of interesting showcases (including the CA13 example shown here) are available on the Potree page.

In the near future the Potree software will be distributed together with the LAStools package to offer a one-click solution for generating Webportals that host and distribute large LiDAR data sets and offer interactive online visualization and exploration. Potree is open source software that is free for anyone to acquire and to deploy. Please remember that using open source software is not the same as supporting open source software. Given the positive experience that rapidlasso GmbH has had with Potree we can only encourage other geospatial companies to support with time or money those open source projects that help your business.

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 Potree:
Potree is a WebGL based viewer for large point clouds. The project evolved as a Web based viewer from the Scanopy desktop point cloud renderer by TU Wien, Institute of Computer Graphics and Algorithms. It will continue to be free and open source with a FreeBSD license to enable anyone to view, analyze and publicly share their large datasets. Visit http://potree.org for more information.

England Releases National LiDAR DEM with Insane (!) Vertical Resolution

This article could also be titled “How not to implement a national open data policy for massive geospatial data sets” or “Forget single-photon LiDAR, England already has single-quantum LiDAR” … (-:

You may have heard about the amazing open data release by the Environment Agency. So far LiDAR-derived DTM and DSM rasters have been released for 72% of the entire English territory at horizontal resolutions of 50 cm, 1 m, and 2 m. They can be downloaded here. The rasters are distributed as zipped archives of tiles in textual ASC format (*.asc). While easy to parse it would not be our first format of choice for such a large release as it loads slower than a comparable binary format like GeoTIFF or BIL … but so far so good.
Open data download portal for DSM and DTM rasters

Open data download portal for DSM and DTM rasters of England

But here comes the shocker and I would to make this a learning experience for those planning similar download portals. Again, the horizontal resolutions of the DTM and DSM rasters is 50 cm, 1 m, and 2 m. But what vertical resolution was chosen? I can still not quite believe it. It is more than micrometer, more than nanometers, and even more than picometers. I had to look up the name. The vertical resolution ranges from femtometers to attometers. This means that the ASCII numbers that specify the elevation for each grid cell are written down with 15 to 17 digits after the decimal point. Here an overview of units and the corresponding number of digits after the decimal point:

 0 - meters:      1.0
 1 - decimeters:  0.1
 2 - centimeters: 0.01
 3 - millimeters: 0.001
 6 - micrometers: 0.000001
 9 - nanometers:  0.000000001
12 - picometers:  0.000000000001
15 - femtometers: 0.000000000000001
18 - attometers:  0.000000000000000001
Wikipedia states that “The picometre’s length is of an order such that its application is almost entirely confined to particle physics, quantum physics, chemistry and acoustics. Atoms are between 62 and 520 pm in diameter, and the typical length of a carbon-carbon single bond is 154 pm.” and the “femtometer […] was so named in honour of physicist Enrico Fermi, as it is a typical length-scale of nuclear physics. […] For example, the charge radius of a proton is approximately 0.84–0.87 femtometres while the radius of a gold nucleus is approximately 8.45 femtometres.” There is no individual Wikipedia entry for attometers because it’s just too small for most practical use … except for specifying the elevations in the DSM and DTM rasters across England … (-; … this interactive animation gives you a sense of those scales.
a Helium atom has a diameter of about 62 picometers.

diameter of Helium atom =  62 picometers

No seriously. This is a gigantic waste of network bandwidth, storage, and – more importantly – people’s time. Please fix this as soon as possible. Here an example: I downloaded LIDAR-DSM-1M-SP37.zip (237.96 MB compressed) and a quick look at one DSM after unzipping the 100 tiles (1891.13 MB uncompressed) was reason enough for this article:

D:\LAStools\bin>more LIDAR-DSM-1M-SP37\sp3070_DSM_1m.asc
ncols        1000
nrows        1000
xllcorner    430000.000000000000
yllcorner    270000.000000000000
cellsize     1.000000000000
NODATA_value  -9999
 79.9499969482421875 80.23999786376953125 80.95999908447265625 80.9199981689453125 80.90000152587890625 81.44000244140625 80.3300018310546875 79.68000030517578125 79.76000213623046875 79.69000244140625 79.56999969482421875 [...]

If you look at these numbers more carefully you see that they really only ought to have centimeter resolution. I quickly changed the resolution to centimeter with a run of lasgrid on 4 cores:

D:\LAStools\bin>lasgrid -i LIDAR-DSM-1M-SP37\*.asc ^
                      -step 1 -use_bb ^
                      -odir LIDAR-DSM-1M-SP37-NO-FLUFF -oasc ^
                      -cores 4

The result is a DSM that is identical for all practical purposes … just compare the first ten elevations below with those ones above.

D:\LAStools\bin>more LIDAR-DSM-1M-SP37-NO-FLUFF\sp3070_DSM_1m.asc
ncols 1000
nrows 1000
xllcorner 430000.000000
yllcorner 270000.000000
cellsize 1.000000
NODATA_value -9999.0
79.95 80.24 80.96 80.92 80.90 81.44 80.33 79.68 79.76 79.69 79.57 [...]

The resulting 100 *.asc tiles use only 580.45 MB uncompressed on disk: an instant storage saving of nearly 70 percent over those tiles with the insanely high resolution. After compressing them back into a single zipped archive I get a compressed file of size 161.99 MB – still a whopping 32 percent less than the zipped archive that I had originally downloaded.

Environment Agency, please lower the vertical resolution of all your DSM and DTM rasters to centimeters. This will directly translate into enourmous storage and bandwidth savings for you over the coming years with each download being around 30 percent smaller and faster. It will also allow your users to work more efficient with the rasters as decompressing and parsing the files will be quicker. In the future I will happily work with you to pick the perfect format for distributing your soon-to-be-open raw LiDAR points and with all the money you will safe for the storage and tranmission of the rasters you could easily become the third Gold Sponsor of the LASzip LiDAR compressor … (-;

PS: Just curious … which software did you use to generate those insanely high vertical resolutions in the first place?

Realizing Open LiDAR with ArcGIS Online and LASzip

Putting aside for a moment our controversy with ESRI over their “uncool” move to create yet another proprietary format whose only benefit is vendor lock-in (read this, this, or this for more on their “LAZ clone“), they do have some neat Web technology. And maybe not all hope is lost and we can rejoin our efforts for a single “native LAS 1.4 compressor” as suggested a few days ago.

But on to the actual story: The Sonoma County Vegetation Mapping and LiDAR Program is using ESRI’s ArcGIS online Web map technology to implement an interactive download portal for open LiDAR and derived products such as contours, building footprints, and elevation rasters.

The five year SonomaVegMap program – a joint program of Sonoma County Open Space and the Sonoma County Water Agency – aims to survey Sonoma County’s topography, physical and biotic features, and diverse plant communities and habitats. The program has already produced several products including county-wide LiDAR data which is now freely available for download as LASzip compressed LAZ files that are linked in increments of watersheds via an ArcGIS online Web map as shown below.

open LiDAR download portal of SonomaVegMap

download portal of SonomaVegMap project using ArcGIS online

Clicking the download icon opens another window (as shown below) where the user can select either the raw LiDAR or the various derivatives (or all) and start downloading.

selecting the ZIP archive

selecting ZIP file with LASzip-compressed LiDAR in LAZ format for download

After about an hour of downloading and a few more minutes of unzipping the 1.4 GB ‘Bodega Head LAZ Files.zip’ file we finally get to see its contents: 103 LAZ files totalling 1.4 GB that decompress to 296,897,792 LiDAR points.

sonomavegmap_lidar_1

A small personal plea: “Please do not put spaces in file names of folder names but use underscores instead!“. Spaces have traditionally been used to separate commands and having them in file or directory names can break batch processing or python scripting. Using the GUI of lasview we load the 103 files and see how this watershed is tiled:

sonomavegmap_lidar_2

The information bar on the right tells us that the LasMonkey software of Watershed Science Inc. has processed points of type 3 (i.e. with GPS time stamp and RGB color) into California State Plane II coordinates and stored them with a resolution of “centi survey feet” (i.e. scale x y z are all 0.01). We pick a horizontal strip and view the points with different visualization options.

(1) RGB (2) intensity (3) classification (4) elevation (5) flightline (6) user data (7) hillshaded ground

(1) RGB (2) intensity (3) classification (4) elevation (5) flightline (6) user data (7) hillshaded TIN of ground points

Everything looks correct although we are not quite sure which values are stored in the ‘user data’ field of each point. Apparently these values more intense for water returns and maybe someone from Watershed Science Inc. will be able to enlighten us. Below a single tile as well as a close-up of one particular area of this tile followed by a lasinfo report. All 103 tiles pass the lasvalidate test with flying colors, nice job guys!

a closer look at tile 'SOCO_0024_121.laz'

a closer look at tile ‘SOCO_0024_121.laz’ of the Bodega Head watershed

close-up with ground and building points triangulated

close-up with ground and building points triangulated

lasinfo report for tile '   '

lasinfo report for tile ‘SOCO_0024_121.laz’

Let me close this article with a 1 minute and 30 second video that summarizes the SonomaVegMap project.