LASmoons: Leonidas Alagialoglou

Posted on December 17, 2020 by rapidlasso

Leonidas Alagialoglou (recipient of three LASmoons)
Multimedia Understanding Group, Aristotle University of Thessaloniki
Thessaloniki, GREECE

Background:
Canopy height is a fundamental geometric tree parameter in supporting sustainable forest management. Apart from the standard height measurement method using LiDAR instruments, other airborne measurement techniques, such as very high-resolution passive airborne imaging, have also shown to provide accurate estimations. However, both methods suffer from high cost and cannot be regularly repeated.

Preliminary results of predicted CHE based on multi-temporal satellite images against ground-truth LiDAR measurements. The 3rd column depicts pixel-wise absolute error of prediction. Last column depicts pixel-wise uncertainty estimation of the prediction (in means of 3 standard deviations).

Goal:
In our study, we attempt to substitute airborne measurements with widely available satellite imagery. In addition to spatial and spectral correlations of a single-shot image, we seek to exploit temporal correlations of sequential lower resolution imagery. For this we use a convolutional variant of a recurrent neural network based model for estimating canopy height, based on a temporal sequence of Sentinel-2 images. Our model’s performance using sequential space borne imagery is shown to outperform the compared state-of-the-art methods based on costly airborne single-shot images as well as satellite images.

Digital Terrain Model of a part of the study area

Data:
The experimental study area of approximately 940 squared km is includes two national parks, Bavarian Forest National Park and Šumava National Park, which are located at the border between Germany and Czech Republic. LiDAR measurements of the area from 2017 and 2019 will be used as ground truth height measurements that have been provided by the national park’s authorities. Temporal sequences of Sentinel-2 imagery will be acquired from the Copernicus hub for canopy height estimation.

LAStools processing:
Accurate conversion of LAS files into DEM and DSM in order to acquire ground truth canopy height model.
1) Remove noise [lasthin, lasnoise]
2) Classify points into ground and non-ground [lasground, lasground_new]
3) Create DTMs and DSMs [lasthin, las2dem]

LASmoons: Zak Kus

Posted on June 20, 2020 by rapidlasso

Zak Kus (recipient of three LASmoons)
Topology Enthusiast
San Francisco, USA

Background:
While LiDAR data enables a lot of research and innovation in a lot of fields, it can also be used to create unique and visceral art. Using the high resolution data available, a 3D printer, and a long tool chain, we can create a physical, 3D topological map of the San Francisco bay area that shows off both the city’s hilly geology, and its unique skyline.

Test print of San Francisco’s Golden Gate Park.

Test print of San Francisco’s Golden Gate Park.

Goal:
The ultimate goal of this project is to create an accurate, unique physical map of San Francisco, and the surrounding areas, which will be given to a loved one as a birthday gift. Using the data from the 2010 ARRA-CA GoldenGate survey, we can filter and process the raw lidar data into a DEM format using LAStools, which can be converted using a python script into a “water tight” 3D printable STL file.

While the data works fairly well out of the box, it does require a lot of manual editing, to remove noise spikes, and to delineate the coast line from the water in low lng areas. Interestingly, while many sophisticated tools exist to edit STLs that could in theory be used to clean up and prepare the files at the STL stage, few are capable of even opening files with so much detailed data. Using LAStools to manually classify, and remove unwanted data is the only way to achieve the desired level of detail in the final piece.

Data:
+ LiDAR data provided through USGS OpenTopography, using the ARRA-CA GoldenGate 2010 survey
+ Average point density of 3.33 pts/m^2 (though denser around SF)
+ Covers 2638 km^2 in total (only a ~100 km^2 subset is used)

LAStools processing:
1) Remove noise [lasnoise]
2) Manually clean up shorelines and problematic structures [lasview, laslayers]
3) Combine multiple tiles (to fit 3d printer) [lasmerge]
4) Create DEMs (asc format) for external tool to process [las2dem]

LASmoons: Martin Romain

Posted on June 13, 2020 by rapidlasso

Martin Romain (recipient of three LASmoons)
Marshall Islands Conservation Society
Majuro, Republic of the MARSHALL ISLANDS

Background:
As a low-lying coastal nation, the Republic of the Marshall Islands (RMI) is at the forefront of exposure to climate change impacts. RMI has a strong dependence on natural resources and biodiversity not only for food and income but also for culture and livelihood. However, these resources are threatened by rising sea levels and associated coastal hazards (king tides, storm surges, wave run-up, saltwater intrusion, erosion). This project aims at addressing the lack of technical capacity and available data to implement effective risk reduction and adaptation measures, with a particular focus on inundation mapping and local evacuation planning in population centers.

Typical low-lying coastal area of the Republic of the Marshall

Goal:
This project intends to use LAStools to generate a DEM of the inhabited sections of 3 remote atolls (Aur, Ebon, Likiep) and 1 island (Mejit). The resulting DEM will be used to produce an inundation exposure model (and map) under variable sea level rise projections for each site. The ultimate goal is to integrate the results into each site’s disaster risk reduction strategy (long-term outcome) and present it through community consultations in schools, community centers, and council houses.

Data:
+ Aerial imagery of 11.5 square kilometers of land (6.3% of total national landmass) using DJI Matrice 200 V2 & DJI Zenmuse X5S with a minimum overlap of 75/75 and maximum altitude of 120m.

LAStools processing:
1) tile large point cloud into tiles with buffer [lastile]
2) remove noise points [lasthin, lasnoise]
3) classify points into ground and non-ground [lasground]
4) create Digital Terrain Models and Digital Surface Models [lasthin, las2dem]

Potential LAStools pipelines:
1) Removing Excessive Low Noise from Dense-Matching Point Clouds
2) Digital Pothole Removal: Clean Road Surface from Noisy Pix4D Point Cloud
3) Creating DTMs from dense-matched points of UAV imagery from SenseFly’s eBee

Coco Loco for LiDAR: Open Data and 3D GIS Carnival Cruise into Sao Paulo

Posted on May 1, 2020 by rapidlasso

The municipality of São Paulo just dropped a MASSIVE amount of open LiDAR upon the geospatial community – a whopping 33 billion points. São Paulo city now has their own Entwine Point Tiles (EPT) service available for everyone interested in downloading, visualizing, or otherwise accessing this amazing data set. The EPT point cloud service can be accessed here and an interactive 3D portal that allows you to explore this massive data with nothing but a Web browser can be accessed here. It simple serves these EPT tiles via Potree.

sao_paulo_online_open_data_01_rgb

What you see above is the beautiful Matarazzo Building (Edifício Matarazzo in Portuguese) where the City Hall of Sao Paulo is located. On May 27 2019 we had a one day workshop in the similarly beautiful Martinelli building (Edifício Martinelli in Portuguese) to brainstorm how to process the data and I seem to remember being rather VOCAL about repeating the success of Open LiDAR in Guadalajara and also providing easy access to the data for Sao Paolo. KUDOS to those who made it happen.

This data was originally collected for the Urban Agriculture project “São Paulo: Growing Farmers’ Income, Shrinking Urban Sprawl” that was supported after winning the 2016 MAYORS CHALLENGE by Bloomburg Philanthropies. At this point I would like to thank former New York mayor Michael Rubens Bloomburg (“Hey Mike, do they also constantly miss-spell your family name? So annoying!”) for ending his 2020 bid for the Democratic presidential nomination.

This project stands on the shoulders of many many giants such as Howard Butler and his team who worked for years to create Entwine, Markus Schuetz who came out of nowhere and suddenly dropped Potree upon our community, and the many contributors in person at OSGeo hackfests and codesprints or remotely via patches and pull requests on the software infrastructure behind this effort. Don’t forget to support these folks!

point colored by RGB

point colored by elevation

point colored by classification

points colored by intensity

point colored by return type

Above other visualizations of the the Matarazzo Building where the City Hall of Sao Paulo is located. Below an overview of the city after zooming out a bit.

point colored by RGB

point colored by elevation

point colored by classification

Surprise Release of Airborne LiDAR in Germany’s “Closed Data State” Bavaria

Posted on February 17, 2020 by rapidlasso

You have guessed correctly. This is mostly fake news as our “Freistaat” (read “Free State“) of Bavaria continues to tightly guard all of its tax-payer funded geospatial basis data for no good reason. Our other “Free State” – that of Thuringia – has become one poster child of open data in Germany with North Rhine-Westphalia being the original one. But there is “some” open LiDAR in Bavaria now.

The authors of a recent paper on change detection in urban areas have published two interesting airborne LiDAR data sets from 2008 and 2009 for the town of Abenberg in [Hebel, Arens, Stilla, 2013]. What is interesting about these data sets is that they (a) were flown with a forward looking laser scanner with flight trajectories from 4 different directions (as illustrated in the image below) and (b) were surveyed again the same way in the following year for temporal change detection.

Our lasview is able to visualize how the four different flight lines scan this house from four different directions, once we have reconstructed a properly populated LAZ file with flight line information, return numbers, and unique GPS time stamps.

The data is provided for download here as zipped ASCII files that have one line per returns containing 11 comma-separated values. Below is a sample of the first 10 lines of the 2008 data:

1, 290.243, 28.663, -11.787, 0.060, -0.052, 0.997, 517.3170, -58.6934, 313.0817, 52
1, 290.208, 28.203, -11.825, 0.062, -0.056, 0.996, 517.3167, -58.6934, 313.0817, 49
1, 290.182, 27.739, -11.852, 0.063, -0.055, 0.997, 517.3164, -58.6935, 313.0817, 53
1, 290.165, 27.272, -11.866, 0.061, -0.058, 0.996, 517.3161, -58.6935, 313.0817, 53
1, 290.163, 26.800, -11.858, 0.061, -0.053, 0.997, 517.3157, -58.6935, 313.0817, 68
1, 290.152, 26.334, -11.864, 0.059, -0.054, 0.997, 517.3154, -58.6936, 313.0817, 57
1, 290.092, 25.882, -11.938, 0.050, -0.057, 0.997, 517.3151, -58.6936, 313.0817, 57
1, 290.103, 25.406, -11.911, 0.043, -0.058, 0.997, 517.3147, -58.6937, 313.0817, 63
1, 290.067, 24.947, -11.952, 0.043, -0.061, 0.997, 517.3144, -58.6937, 313.0817, 63
1, 290.034, 24.488, -11.989, 0.044, -0.063, 0.997, 517.3141, -58.6937, 313.0817, 56

The first number is either a classification into ground, vegetation, or other surface, or represents an identifier for a planar shape that the return is part of. The next three numbers in red are the x, y, and z coordinate of the LiDAR point in a local coordinate system. The next three numbers in green are the x, y, and z coordinates of an estimated surface normal. The next three numbers in blue are the x, y, and z coordinates of the sensor position. The last number is the intensity of the LiDAR return.

This textual representation makes it difficult to efficiently load the data into most LiDAR processing software. Also several attributes such as return number, number of returns, flight line ID, and GPS time stamps are missing.

We use this as an opportunity for a little exercise in converting ASCII to LAZ while preserving any “additional attributes” using the “extra bytes” functionality available since the LAS 1.4 specification. This is a timely experiment as the LAS Working Group of the ASPRS is currently contemplating how to standardize some useful “additional attributes”. Here you can download the resulting files:

In order to replicate these steps, please get your hands of the very latest version of LAStools. First we use txt2las to convert the TXT file to LAZ format as follows:

txt2las -i abenberg_data_2008.txt ^
        -set_point_type 1 ^
        -parse 0xyz123456i ^
        -set_scale 0.001 0.001 0.001 ^
        -add_attribute 3 "planar shape ID" "preliminary classification" ^
        -add_attribute 4 "normal x coord" "local normal direction estimate" 0.001 ^
        -add_attribute 4 "normal y coord" "local normal direction estimate" 0.001 ^
        -add_attribute 4 "normal z coord" "local normal direction estimate" 0.001 ^
        -add_attribute 6 "sensor x coord" "sensor position" 0.0001 ^
        -add_attribute 6 "sensor y coord" "sensor position" 0.0001 ^
        -add_attribute 6 "sensor z coord" "sensor position" 0.0001 ^
        -odix _temp1 -olaz

We use a back and forth of lasview and las2las with option ‘-subseq 12345 67890’ to interactively find the exact index of the return where each flightline is ending and the next one is starting. The command below allows you to visualize the the trajectories.

lasview -i abenberg_data_2008_temp1.laz ^
        -copy_attribute_into_x 4 ^
        -copy_attribute_into_y 5 ^
        -copy_attribute_into_z 6 ^
        -points 10000000 ^
        -point_size 5

abenberg_open_data_lasview_trajectory

By hovering with the mouse over a point where the trajectory end and pressing ‘i’ like info we can query the coordinates and attributes of a point. The console output also lists the index of the point in the file. We use this index as the start index for the manual search of the exact index where the flight lines really ends by piping the coordinates as text to stdout and looking for a “jump” in coordinates that indicates the start of a new flightline as shown below.

las2las -i abenberg_data_2008_temp1.laz ^
        -subseq 1213485 1213505 ^
        -oparse xyz -otxt -stdout

-279.846 -98.442 -11.973
-279.984 -98.900 -12.123
-280.150 -99.357 -12.311
-280.195 -99.825 -12.332
-280.229 -100.294 -12.340
-280.275 -100.763 -12.363
-280.302 -101.233 -12.361
-280.344 -101.700 -12.379
-280.371 -102.171 -12.376
-280.156 -102.661 -12.042
110.259 304.077 3.873
110.646 304.118 3.925
111.036 304.154 3.970
111.424 304.167 3.982
111.811 304.211 4.037
112.201 304.252 4.088
112.588 304.278 4.118
112.976 304.315 4.164
113.364 304.336 4.186
113.752 304.344 4.192

Then we use the found index to seperate the first flight line form the rest with with two more runs of las2las:

las2las -i abenberg_data_2008_temp1.laz ^
         -subseq 0 1213495 ^
         -odix _strip1 -olaz

las2las -i abenberg_data_2008_temp1.laz ^
        -subseq 1213495 100000000 ^
        -odix _rest234 -olaz

We then repeat this procedure for the rest until we have four individual strips. Next we use las2las to change the point type from 0 to 1 to have a GPS time attribute that we will then populate with “fake” but useful GPS time stamps:

las2las -i abenberg_data_2008_temp1_strip*.laz ^
        -set_point_type 1 ^
        -odix _pt1 -olaz

We noticed in the original text file that subsequent groups of points often have the exact same value for the three numbers in blue that are the x, y, and z coordinates of the sensor position. This suggests that those points are multiple returns from the same laser shot. We wrote a small tool using LASlib that exploits this observed pattern to recover the “return number” and the “number of returns” attribute for each point and to set a “fake” but unique GPS time for each such group of returns. You can download the source code for this tool here. We run this tool for each strip with a different start GPS time:

lasrecover -i abenberg_data_2008_temp1_strip1_pt1.laz ^
           -gpstime_start 1000000 ^
           -odix _rec -olaz

lasrecover -i abenberg_data_2008_temp1_strip2_pt1.laz ^
           -gpstime_start 2000000 ^
           -odix _rec -olaz

lasrecover -i abenberg_data_2008_temp1_strip3_pt1.laz ^
           -gpstime_start 3000000 ^
           -odix _rec -olaz

lasrecover -i abenberg_data_2008_temp1_strip4_pt1.laz ^
           -gpstime_start 4000000 ^
           -odix _rec -olaz

Finally we merged the four strips back into one file using lasmerge:

lasmerge -i abenberg_data_2008_temp1_strip1_pt1_rec.laz ^
         -i abenberg_data_2008_temp1_strip2_pt1_rec.laz ^
         -i abenberg_data_2008_temp1_strip3_pt1_rec.laz ^
         -i abenberg_data_2008_temp1_strip4_pt1_rec.laz ^
         -faf ^
         -o abenberg_data_2008_temp2.laz

With lasview we are now able to visualize not only the multiple returns per shot but also the different angles from which the laser scanner has observed the scene.

flight line 4 scan direction

flight line 2 scan direction

flight line 3 scan direction

Finally we use las2las and a “filtered transform” with the custom classfication code provided in the first “additional attribute” to populate the official LAS classification codes as ratified by the ASPRS. First we turn code 1 (“Ground level”) into ground points.

las2las -i abenberg_data_2008_temp2.laz ^
        -keep_attribute_between 0 1 1 ^
        -filtered_transform ^
        -set_classification 2 ^
        -o abenberg_data_2008_temp3.laz

Then we turn code 5 (“Vegetation”) into high vegetation points, code 6 (“Other Surface”) into keypoints, and code 9 or higher (“Planar shape #”) into building points.

las2las -i abenberg_data_2008_temp3.laz ^
        -keep_attribute_between 0 5 5 ^
        -filtered_transform ^
        -set_classification 5 ^
        -o abenberg_data_2008_temp4.laz

las2las -i abenberg_data_2008_temp4.laz ^
        -keep_attribute_between 0 6 6 ^
        -filtered_transform ^
        -set_classification 8 ^
        -o abenberg_data_2008_temp5.laz

las2las -i abenberg_data_2008_temp5.laz ^
         -keep_attribute_above 0 8 ^
         -filtered_transform ^
         -set_classification 6 ^
         -o abenberg_data_2008.laz

Then we are done. Here you can download the resulting files:

colored by classification

colored by return type

colored by flight line

References

Hebel, M., Arens, M., Stilla, U., 2013. Change detection in urban areas by object-based analysis and on-the-fly comparison of multi-view ALS data. ISPRS Journal of Photogrammetry and Remote Sensing 86, pp. 52-64. [DOI: 10.1016/j.isprsjprs.2013.09.005] [PDF]

Philippines use Taal Vulcano Eruption as Opportunity to become Very First Asian Country with Open LiDAR

Posted on January 28, 2020 by rapidlasso

UPDATE: As of January 30th also orthophotos and classified LAZ tiles are available for download.

It took just a few years of nagging, a vulcanic eruption, and then a few more weeks of nagging but now it has happened. The Philippines have become the first country in Asia to offering LiDAR as open data for free and unencumbered download. The portal created by the UP Training Center for Applied Geodesy and Photogrammetry (UP TCAGP) and their DREAM and PHIL LiDAR program already offers LiDAR-derived 1 meter DTM and DSM data flown between 2013 and 2017 as part of a national mission to aquire flood mapping data for a certain area around the Taal Vulcano. In the coming days orthophotos and the classified LiDAR point cloud will be added (at the moment the data is still undergoing another quality assurance review process).

As a quick test we went to the new online portal and downloaded the 34 DTM raster tiles that cover the Taal Vulcano Lake as seen in the screenshot below.

Downloading the area-of-interest is easy with LiPAD’s nice download portal.

The downloaded 1 meter DTM tiles are in TIF format and each cover an area of 1000 by 1000 meter. However, they are overlapping because they have a 50 meter buffer, so that each raster contains elevation samples organized in 1100 columns by 1100 rows plus “no data” values. We use two LAStools commands to remove the buffers. First we use our new demzip to turn the TIF to RasterLAZ format. Use demzip from version 200131 of LAStools (or newer) as older releases did not handle “no data” values correctly.

demzip -i Taal\DTM\*.tif ^
       -olaz

The conversion from TIF to RasterLAZ also reduces the total file size for the 34 files from 157 MB to 27 MB. Next we remove the buffers using a new functionality in lasgrid (make sure you have the latest LAStools version 200112 or newer).

lasgrid -i Taal\DTM\*.laz ^
        -step 1 ^
        -use_tile_size 1000 ^
        -odir Taal\DTM_unbuffered ^
        -olaz

Without buffers the total file size in RasterLAZ format shrinks to 22 MB. Now we have the data in a format that can either be treated as a raster or as a point cloud. Hence we can use laspublish and quickly create a visualization of the Taal Vulcano Island with Potree which we then copied onto our university Web space for you to play with. This was he are able to instantly create an 3D visualization portal that lets anyone do various simple and also more complex measurements.

laspublish -i Taal\DTM_unbuffered ^
           -elevation ^
           -odir Taal\DTM_portal ^
           -o TaalVulcanoIsland.html ^
           -title "DTM of Taal Vulcano Island" ^
           -description "DTM of Taal Vulcano Island" ^
           -olaz -overwrite

Below we see the result visualized with the Desktop version of Potree. You can access the interactive portal we have created here with any Web browser.

Visualizing the 1 meter DTM of Taal Vulcano Island as RasterLAZ point cloud with Potree to instantly create interactive portal allowing simple measurements that give an intuition about the height and the size of the vulcanic formation that makes up Taal Vulcano Island.

We would like to acknowledge the UP Training Center for Applied Geodesy and Photogrammetry (UP TCAGP) and their DREAM and PHIL LiDAR program for providing easy and unencumbered open access to this data with a license that encourages data reuse and repurposing. Kudos for being first in Asia to make open LiDAR happen!!!

Converting Rasters from inefficient ASCII XYZ to more compact LAZ or TIF Formats

Posted on January 13, 2020 by rapidlasso

The German state of Brandenburg has recently started to provide many of their basic geospatial data as open data, such as digital ortophotos in TIF and JPG formats, vertical and horizontal control points in gzipped XML format, LOD1 and LOD2 building models in zipped GML format, topographic maps from 1:10000 to 1:100000 in zipped TIF and PDF formats, cadastral data in zipped XML and TIF formats, as well as LiDAR-derived 1m DTM rasters and image-derived 1m DSM rasters both in zipped XYZ ASCII format. All this data is provided with the user-friendly license called “Datenlizenz Deutschland Namensnennung 2.0“. In this article we show how to convert the 1m DTM rasters and the 1m DSM rasters from verbose XYZ ASCII to more compact LAZ or TIF rasters.

Four 2000 by 2000 meter tiles of the Brandenburg 1m DTM.

One particularity about most official German and Austrian rasters (anywhere else?) is that they sample the elevations in the corners rather than in the center of each raster cell. Here a one square kilometer raster tile of 1 meter resolution will have 1001 columns by 1001 rows instead of the more familiar 1000 by 1000 layout. While this corner-based representation does have some benefits, we convert these rasters in to the more common area-based representation using new functionality recently added to lasgrid.

After downloading one sample DTM tile such as dgm_33250-5886.zip we find three files in the zip folder. Two files with meta data and license information and the actual data file, which is a 2 km by 2km corner-based raster tile called “dgm_33250-5886.xyz” with 2001 columns by 2001 rows. Here is how the 4004001 lines looks:

more DGM_33250-5886.xyz
250000.0 5886000.0 15.284
250001.0 5886000.0 15.277
250002.0 5886000.0 15.273
250003.0 5886000.0 15.275
250004.0 5886000.0 15.289
250005.0 5886000.0 15.314
[...]
251994.0 5888000.0 13.565
251995.0 5888000.0 13.567
251996.0 5888000.0 13.565
251997.0 5888000.0 13.565
251998.0 5888000.0 13.564
251999.0 5888000.0 13.564
252000.0 5888000.0 13.565

The first step is to convert these XYZ rasters to LAZ format. We do this with txt2las as shown below. In case the vertical datum is the “Deutsches Haupthoehennetz 2016” we should also add ‘-vertical_dhhn2016’ but not sure at the moment:

txt2las -i dgm\*.xyz ^
        -set_scale 1.0 1.0 0.001 ^
        -epsg 25833 ^
        -odir temp -olaz ^
        -cores 4

For 84 files this reduces the size by a factor of 31 or compresses it down to 3.2 percent of the original, namely from 8.45 GB for raw XYZ to 277 MB for LAZ. So far we have really just converted a list of x, y and z coordinates from verbose ASCII to more compact LAZ. We can easily go back to ASCII with las2txt whenever needed:

txt2las -i temp\*.laz ^
        -odir ascii -otxt ^
        -cores 4

Next we use lasgrid to convert from a corner-based raster to an area-based raster using the new option ‘-subsquare 0.2’ which replaces each input point by four points that are displaced by all possibilities of adding +/- 0.2 in x and y. We then average the exactly four points that fall into each relevant raster cell with option ‘-average’ and clip the output to the meaningful 2000 columns by 2000 rows with ‘-use_tile_size 2000’. You need to get the most recent version of LAStools to have these options.

lasgrid -i temp\*.laz ^
        -subsquare 0.25 ^
        -step 1 -average ^
        -use_tile_size 2000 ^
        -odir dgm -olaz ^
        -cores 4

Instead of RasterLAZ you can also choose the TIF, BIL, IMG, or ASC format here. The final result are standard 1 meter elevation products with 2000 columns by 2000 rows with the averaged elevation sample being associated with the center of the raster cell. The lasinforeport for a sample tile is shown at the end of this article.

You may proceed to optimize the RasterLAZ for area-of-interest queries by reordering the raster into a space-filling curve with lassort or lasoptimize and compute a spatial index. You may also classify the RasterLAZ elevation samples, for example, into building, high, medium, and low vegetation, ground, and other common classifications with lasclip or lascolor. You may also add RGB or intensity values to the RasterLAZ elevation samples using the orthophotos that are also available as open data with lascolor. These are some of the benefits of RasterLAZ beyond efficient storage and access.

We like to acknowledge the LGB (Landesvermessung und Geobasisinformation Brandenburg) for providing state-wide coverage of their geospatial data holdings as easily downloadable open data with the user-friendly Deutschland Namensnennung 2.0 license. But we also would like to ask to please add the raw LiDAR point clouds to the open data portal. The storage savings in going from ASCII XYZ to LAZ for the DTM and DSM rasters should free enough space to host the LiDAR … (-;

lasinfo (200112) report for 'dgm_33\DGM_33250-5886.laz'
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.2
  system identifier:          'raster compressed as LAZ points'
  generating software:        'LAStools (c) by rapidlasso GmbH'
  file creation day/year:     13/20
  header size:                227
  offset to point data:       455
  number var. length records: 2
  point data format:          0
  point data record length:   20
  number of point records:    4000000
  number of points by return: 4000000 0 0 0 0
  scale factor x y z:         0.5 0.5 0.001
  offset x y z:               200000 5800000 0
  min x y z:                  250000.5 5886000.5 13.419
  max x y z:                  251999.5 5887999.5 33.848
variable length header record 1 of 2:
  reserved             0
  user ID              'Raster LAZ'
  record ID            7113
  length after header  80
  description          'by LAStools of rapidlasso GmbH'
    ncols   2000
    nrows   2000
    llx   250000
    lly   5886000
    stepx    1
    stepy    1
    sigmaxy <not set>
variable length header record 2 of 2:
  reserved             0
  user ID              'LASF_Projection'
  record ID            34735
  length after header  40
  description          'by LAStools of rapidlasso GmbH'
    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 25833 - ProjectedCSTypeGeoKey: ETRS89 / UTM 33N
      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
LASzip compression (version 3.4r3 c2 50000): POINT10 2
reporting minimum and maximum for all LAS point record entries ...
  X              100001     103999
  Y              172001     175999
  Z               13419      33848
  intensity           0          0
  return_number       1          1
  number_of_returns   1          1
  edge_of_flight_line 0          0
  scan_direction_flag 0          0
  classification      0          0
  scan_angle_rank     0          0
  user_data           0          0
  point_source_ID     0          0
number of first returns:        4000000
number of intermediate returns: 0
number of last returns:         4000000
number of single returns:       4000000
overview over number of returns of given pulse: 4000000 0 0 0 0 0 0
histogram of classification of points:
         4000000  never classified (0)

LASmoons: Volga Lipwoni

Posted on December 18, 2019 by rapidlasso

Volga Lipwoni (recipient of three LASmoons)
Department of Geography, School of Earth and Environment
University of Canterbury, NEW ZEALAND

Background:
Structure from motion (SfM) photogrammetry, has emerged as an effective tool to accurately extract three-dimensional (3D) structures from a series of overlapping two-dimensional (2D) Unmanned aerial vehicles (UAVs) images. The bid to switch from the current labour-intensive, and time consuming forestry inventory practices has seen a lot of interest geared towards understanding the use of SfM photogrammetry to derive forest metrics (Iglhaut et al., 2019). There are a range of commercial, free and open source SfM photogrammetric software packages that can be used to process UAV images into 3D point clouds. Selection of the most appropriate package has become an important issue for most projects (Turner, Lucieer, & Wallace, 2013). A comparison of software performance in terms of accuracy, processing times and related costs would help foresters in deciding the best tool for the job.

Typical point cloud derived with SfM software from UAV imagery.

Goal:
The study will generate 3D point clouds of images of a young forest trial and LAStools will be used to derive canopy height models (CHM) for computing tree heights. Tree heights from LiDAR data will serve as a baseline for accuracy assessment of heights derived from the point clouds.

Data:
+ 422 UAV images processed into 3D point clouds using ten (10) different commercial and open source SfM software packages

References:
Iglhaut, J., Cabo, C., Puliti, S., Piermattei, L., O’Connor, J., & Rosette, J. (2019). Structure from motion photogrammetry in forestry: A review. Current Forestry Reports, 5(3), 155-168. doi:https://doi.org/10.1007/s40725-019-00094-3
Turner, D., Lucieer, A., & Wallace, L. (2013). Direct georeferencing of ultrahigh-resolution UAV imagery. EEE Transactions on Geoscience and Remote Sensing, 52(5), 2738-2745. doi:10.1109/TGRS.2013.2265295

Removing Noise from Single Photon LiDAR to Generate a Smooth DTM

Posted on December 16, 2019 by rapidlasso

A while back we had a first look at the Single Photon LiDAR from Leica’s SPL100 sensor (that eventually turned out just to be an SPL99 because one beamlet or one receiver in the 10 by 10 array was broken and did not produce any returns). Today we are taking a closer look at a strategy to remove the excessive noise in the raw Single Photon LiDAR data from a “proper” SPL100 sensor (where all of the 100 beamlets are firing) that was flown in 2017 in Navarra, Spain.

Profile through original points on top of generated DTM.

The data was provided as open data by the cartography section of Navarra’s Government and is available via a simple download FTP portal. We describe the LAStools processing steps that were used to eliminate the excessive noise and to generate a smooth DTM. In the following we are using the originally released version of the data, that we obtained shortly after the portal went online that seems to be a bit more “raw” than the current files available now. One starndard quality check with lasinfo was done with:

lasinfo -i 0_raw\*.laz ^
        -cd ^
        -histo intensity 1 ^
        -histo user_data 1 ^
        -histo point_source 1 ^
        -histo gps_time 10 ^
        -odir 1_quality -odix _info -otxt

Upon inspecting the lasinfo report we suggest a few changes in how to store this Single Photon LiDAR data for more efficient hosting via an online portal. We perform these changes here before starting the actual processing. First we use the las2las call shown below to fix an error in the global encoding bits, remove an irrelevant VLR, re-scale the coordinates from millimeter to centimeters, re-offset the coordinates to nice numbers, and – what is by far the most crucial change for better compression – remap the beamlet ID stored in the ‘user data’ field as described in an earlier article.

las2las -i 0_raw\*.laz ^
        -rescale 0.01 0.01 0.01 ^
        -auto_reoffset ^
        -set_global_encoding_gps_bit 1 ^
        -remove_vlr 1 ^
        -map_user_data beamlet_ID_map.txt ^
        -odir 2_fix_rescale_reoffset_remap -olaz ^
        -cores 3

Then we use two lassort calls, one to maximize compression and one to improve spatial coherence. One lassort call rearranges the points in increasing order first based on the GPS time stamps, then breaks ties based on the user data field (that stores the beamlet ID), and finally stores the returns of every beamlet ordered by return number. We also add spatial reference information in this step. The other lassort call rearranges the points into a spatially coherent layout. It uses a Z-order sort with the granularity of 50 meter by 50 meter buckets of points. Within each bucket the point order from the prior sort is kept.

lassort -i 2_fix_rescale_reoffset_remap\*.laz ^
        -epsg 25830 ^
        -gps_time ^
        -user_data ^
        -return_number ^
        -odir 2_maximum_compression -olaz ^
        -cores 3

lassort -i 2_maximum_compression\*.laz ^
        -bucket_size 50 ^
        -odir 2_spatial_coherence -olaz ^
        -cores 3

The resulting optimized nine tiles are around 200 MB each and can be downloaded as one file here or as individual tiles here:

Now we start the usual processing workflow by tiling the data with lastile into smaller 500 meter by 500 meter tiles with a 25 meter buffer. We also set the pre-existing point classification in the data to zero as we will compute our own later.

lastile -i 2_spatial_coherence\*.laz ^
        -set_classification 0 ^
        -tile_size 500 -buffer 25 -flag_as_withheld ^
        -odir 3_buffered -o yecora.laz

We notice that a large amount of the noise has intensity values below 1000. We are still a bit puzzled where those intensity values come from and what exactly they mean in a Single Photon LiDAR system. But it works. We run las2las with a “filtered transform” to set classification of all points whose intensity value is 1000 or less to the classification code 7 (aka “noise”).

las2las -i 3_buffered\*.laz ^
        -keep_intensity_below 1000 ^
        -filtered_transform ^
        -set_classification 7 ^
        -odir 4_intensity_denoised -olaz ^
        -cores 3

We then ignore this “easy-to-identify” noise and go after the remaining one with lasnoise by ignoring classification code 7 and setting the newly identified noise to classification code 9 – not because it’s “water” (the usual meaning of class 9) but because these points are drawn with a distinct blue color when checking the result with lasview.

 lasnoise -i 4_intensity_denoised\*.laz ^
         -ignore_class 7 ^
         -step_xy 1.0 -step_z 0.2 ^
         -isolated 5 ^
         -classify_as 9 ^
         -odir 4_isolation_denoised -olaz ^
         -cores 3

Of the surviving non-noise points we then use lasthin to reclassify the point closest to the 20th elevation percentile per 50 cm by 50 cm area with classification code 8 (for all areas that have more than 5 non-noise points per 50 cm by 50 cm area. We repeat the same for every 1 meter by 1 meter area.

lasthin -i 4_isolation_denoised\*.laz ^
        -ignore_class 7 9 ^
        -step 0.5 -percentile 20 5 ^
        -classify_as 8 ^
        -odir 5_thinned_p20_050cm -olaz ^
        -cores 3

lasthin -i 5_thinned_p20_050cm\*.laz ^
        -ignore_class 7 9 ^
        -step 1.0 -percentile 20 5 ^
        -classify_as 8 ^
        -odir 5_thinned_p20_100cm -olaz ^
        -cores 3

We then perform a more agressive second noise removal step one with lasnoise using only those points with classification code 8, namely those non-noise points that were the 20th elevation percentile in either a 50 cm by 50 cm cell or a 1 meter by 1 meter cell. This can be done by ignoring classification code 0, 7, and 9. We mark those noise points as 6 so they appear orange in the point cloud with lasview.

lasnoise -i 5_thinned_p20_100cm\*.laz ^
         -ignore_class 0 7 9 ^
         -step_xy 2.0 -step_z 0.2 ^
         -isolated 1 ^
         -classify_as 6 ^
         -odir 5_thinned_p20_100cm_denoised -olaz ^
         -cores 3

The 20th elevation percentile points that survive the last noise removal are then classified into ground (2) and non-ground (1) points with lasground_new by ignoring all other points, namely those with classification codes 0, 6, 7, and 9.

lasground_new -i 5_thinned_p20_100cm_denoised\*.laz ^
              -ignore_class 0 6 7 9 ^
              -town ^
              -odir 5_tiles_ground_050cm -olaz ^
              -cores 3

These images below illustrate the steps we took. They also show that not all data was used and might give you ideas where to tweak our workflow for even better results.

First returns are red, intermediate returns are green, last returns are blue and single returns are yellow.

Colored by intensity. Most noise has low intensity values.

Final classification.

Easy-to-remove noise (all intensity values below 1000).

Remaining noise points removed with lasnoise.

Unused points. Another workflow could try to use more of this data.

Final ground points used for generating the smooth DTM.

Finally we raster the ground points into 1 meter Digital Terrain Model (DTM) rasters with las2dem and store the result (without buffers) to the RasterLAZ format.

las2dem -i 5_tiles_ground_050cm\*.laz ^
        -keep_class 2 ^
        -step 1.0 ^
        -use_tile_bb ^
        -odir 6_tiles_dtm_100cm -olaz ^
        -cores 3

Finally we merged all RasterLAZ tiles into one and compute the final hillshaded DTM with blast2dem.

blast2dem -i 6_tiles_dtm_100cm\*.laz -merged ^
          -step 1.0 ^
          -hillshade ^
          -o yecora_dtm_100cm.png

The hillshaded DTM that is result of the entire sequence of processing steps described above is shown below.

DTM from ground classification created with LAStools

For comparison we generate the same DTM using the originally provided classification. According to the README file the original ground points are classified with code 22 in areas of flight line overlap and as the usual code 2 elsewhere. Hence we must use both classification codes to construct the DTM. We do this analogue to the earlier processing steps with the three LAStools commands lastile, las2dem, and blast2dem below.

lastile -i 2_spatial_coherence\*.laz ^
        -tile_size 500 -buffer 25 -flag_as_withheld ^
        -odir 3_tiles_buffered_orig -o yecora.laz

las2dem -i 3_tiles_buffered_orig\*.laz ^
        -keep_class 2 22 ^
        -step 1.0 ^
        -use_tile_bb ^
        -odir 6_tiles_dtm_100cm_orig -olaz ^
        -cores 3

blast2dem -i 6_tiles_dtm_100cm_orig\*.laz -merged ^
          -step 1.0 ^
          -hillshade ^
          -o yecora_dtm_100cm_orig.png

Below the hillshaded DTM generated from the ground classification that was provided with the LiDAR when it was originally released as open data.

DTM from ground classification of originally released data.

In the meantime Andorra’s SPL data have been updated with a newer version in the open data portal. The new version of the data contains a much better ground classification that might have been improved manually as the new files now have the the string ‘cam’ instead of ‘ca’ in the file name, which probably means ‘classified automatically and manually’ instead of the original ‘classified automatically’. We decided not to switch to the new data release as it seemed less “raw” than the original release. For example there are suddenly points with GPS times and returns counts and numbers of zero in the file that seem synthetic. But we also computed the hillshaded DTM for the new release which is shown below.

DTM from ground classification of newly released data.

We thank the cartography section of Navarra’s Government for providing their LiDAR as open data. This not only allows re-purposing expensive data paid for by public taxes but also generates additional value, encourages citizen science, and provides educational opportunity and insights such as this blog article.

Another European Country Opens LiDAR: Welcome to the Party, Slovakia!

Posted on December 11, 2019 by rapidlasso

We got a little note from Vítězslav Moudrý from CULS pointing out that the Geodesy, Cartography and Cadastre Authority of the Slovak Republic has started releasing LiDAR as open data on their interactive Web portal. Congratulations, Slovakia!!! Welcome to the Open Data Party!!! We managed to download some data starting from this Web portal link and describe the process of obtaining LiDAR data from the Low Tatras mountain range in central Slovakia with pictures below.

(1) click the new “data export” link

(2) change the export selection to “Shape”

(3) change the file format to “LAZ”

(4) zoom to a colored area-of-interest

(5) zoom further and draw a nice polygon

(6) edit polygon into nice shape and realize heart is red because area is too big

(7) zoom further and draw polygon smaller than 2 square kilometer

(8) when polygon turns green, accept license, enter email address and export

(9) short wait and you get download link to such an archive

(10) license conditions: PDF auto-translated from Slovak to English

(11) LiDAR are spatially indexed flight lines clipped to area-of-interest

open_data_portal_slovakia_12_density_all_returns_20_50

(12) all return density: blue = 20 and red = 50 returns per square meter

lasgrid -i LowTatras\*.laz -merged ^
        -step 2 -point_density_16bit ^
        -false -set_min_max 20 50 ^
        -o LowTatras\density_all_returns_20_50.png

open_data_portal_slovakia_13_density_last_returns_4_40

(13) last return density: blue = 4 and red = 40 last returns per square meter

lasgrid -i LowTatras\*.laz -merged ^
        -keep_last ^
        -step 2 -point_density_16bit ^
        -false -set_min_max 4 40 ^
        -o LowTatras\density_last_returns_4_40.png

open_data_portal_slovakia_14_density_ground_returns_4_40

(14) ground return density: blue = 4 and red = 40 ground returns per square meter

lasgrid -i LowTatras\*.laz -merged ^
        -keep_classification 2 ^
        -step 2 -point_density_16bit ^
        -false -set_min_max 4 40 ^
        -o LowTatras\density_ground_returns_4_40.png

open_data_portal_slovakia_14_overlap_10cm_20cm_diff

(15) flight line difference image: white <= +/- 10 cm and red/blue >= +/- 20 cm

lasoverlap -i LowTatras\*.laz -faf ^
           -drop_classification 7 18 ^
           -min_diff 0.1 -max_diff 0.2 ^
           -o LowTatras\overlap_10cm_20cm.png

Finally we compute a DSM and a corresponding DTM using the already existing ground classification with BLAST using the command sequence shown below.

hillshaded 1m DSM

hillshaded 1m DTM

lasthin -i LowTatras\*.laz -merged ^
        -drop_classification 7 18 ^
        -step 0.5 -highest ^
        -o LowTatras\highest_50cm.laz

blast2dem -i LowTatras\highest_50cm.laz ^
          -hillshade ^
          -o LowTatras -o dsm_1m_hillshaded.png

blast2dem -i LowTatras\*.laz -merged ^
          -keep_classification 2 ^
          -thin_with_grid 0.5 ^
          -hillshade ^
          -o LowTatras\dtm_1m_hillshaded.png

We thank the Geodesy, Cartography and Cadastre Authority of the Slovak Republic for providing their LiDAR as open data for both commercial and non-commercial purposes and name the source of the data used above (as the license requires) as the Office of Geodesy, Cartography and Cadastre of the Slovak Republic (GCCA SR) or – in Slovak – the Úrad geodézie, kartografie a katastra Slovenskej republiky (ÚGKK SR).

Which European country goes next? Czech Republic? Poland? Hungary? Switzerland?

fast tools to catch reality

References