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

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.

taal_vulcano_open_lidar_download_portal

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.

taal_vulcano_open_lidar_dtm

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

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.

brandenburg_dgm_258_5888_4000

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)

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

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.

open_data_portal_slovakia_01

(1) click the new “data export” link

open_data_portal_slovakia_02

(2) change the export selection to “Shape”

open_data_portal_slovakia_03

(3) change the file format to “LAZ”

open_data_portal_slovakia_04

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

open_data_portal_slovakia_05

(5) zoom further and draw a nice polygon

open_data_portal_slovakia_06

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

open_data_portal_slovakia_07

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

open_data_portal_slovakia_08

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

open_data_portal_slovakia_09

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

open_data_portal_slovakia_10

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

 

open_data_portal_slovakia_11

(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.

 

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?

 

 

Another German State Goes Open LiDAR: Saxony

Finally some really good news out of Saxony. 😊 After North Rhine-Westphalia and Thuringia released the first significant amounts of open geospatial data in Germany in a one-two punch in January 2017, we now have a third German state opening their entire tax-payer-funded geospatial data holdings to the tax-paying public via a simple and very easy-to-use online download portal. Welcome to the open data party, Saxony!!!

Currently available via the online portal are the LiDAR-derived raster Digital Terrain Model (DTM) at 1 meter resolution (DGM 1m) for everything flown since 2015 and and at 2 meter resolution (DGM 2m) or 20 meter resolution (DGM 20m) for the entire state. The horizontal coordinates use UTM zone 33 with ETRS89 (aka EPSG code 25833) and the vertical coordinate uses the “Deutsche Haupthöhennetz 2016” or “DHHN2016” (aka EPSG code 7837). Also available are orthophotos at 20 cm (!!!) resolution (DOP 20cm).

dgm_1000_rdax_87

Overview of current LiDAR holdings. Areas flown 2015 or later have LAS files and 1 meter rasters. Others have LiDAR as ASCII files and lower resolution rasters.

Offline – by ordering through either this online form or that online form – you can also get the 5 meter DTM and the 10 meter DTM, the raw LiDAR point clouds, LiDAR intensity rasters, hill-shaded DTM rasters, as well as the 1 meter and the 2 meter Digital Surface Model (DSM) for a small administrative fee that ranges between 25 EUR and 500 EUR depending on the effort involved.

Our immediate thought is to get a copy on the entire raw LiDAR points clouds (available as LAS 1.2 files for all  data acquired since 2015 and as ASCII text for earlier acquisitions) and find some portal willing to hosts this data online. We are already in contact with the land survey of Saxony to discuss this option and/or alternate plans.

Let’s have a look at the data. First we download four 2 km by 2 km tiles of the 1 meter DTM raster for an area surrounding the so called “Greifensteine” using the interactive map of the download portal, which are provided as simple XYZ text. Here a look at the contents of one ot these tiles:

more Greifensteine\333525612_dgm1.xyz
352000 5613999 636.26
352001 5613999 636.27
352002 5613999 636.28
352003 5613999 636.27
352004 5613999 636.24
[...]

Note that the elevation are not sampled in the center of every 1 meter by 1 meter cell but exactly on the full meter coordinate pair, which seems especially common  in German-speaking countries. Using txt2las we convert these XYZ rasters to LAZ format and add geo-referencing information for more efficient subsequent processing.

txt2las -i greifensteine\333*_dgm1.xyz ^
        -set_scale 1 1 0.01 ^
        -epsg 25833 ^
        -olaz

Below you see that going from XYZ to LAZ reduces the amount of  data from 366 MB to 10.4 MB, meaning that the data on disk becomes over 35 times smaller. The ability of LASzip to compress elevation rasters was first noted during the search for missing airliner MH370 and resulted in our new LAZ-based compressor for height grid called DEMzip.  The resulting LAZ files now also include geo-referencing information.

96,000,000 333525610_dgm1.xyz
96,000,000 333525612_dgm1.xyz
96,000,000 333545610_dgm1.xyz
96,000,000 333545612_dgm1.xyz
384,000,000 bytes

2,684,820 333525610_dgm1.laz
2,590,516 333525612_dgm1.laz
2,853,851 333545610_dgm1.laz
2,795,430 333545612_dgm1.laz
10,924,617 bytes

Using blast2dem we then create a hill-shaded version of the 1 meter DTM in order to overlay a visual representation of the DTM onto Google Earth.

blast2dem -i greifensteine\333*_dgm1.laz ^
          -merged ^
          -step 1 ^
          -hillshade ^
          -o greifensteine.png

Below the result that nicely shows how the penetrating laser of the LiDAR allows us to strip away the forest to see interesting geological features in the bare-earth terrain.

In a second exercise we use the available RGB orthophoto images to color one of the DTM tiles and explore it using lasview. For this we download the image for the top left of the four tiles that covers the area containing the “Greifensteine” from the interactive download portal for orthophotos. As the resolution of the TIF image is 20 cm and that of the DTM is only 1 meter, we first down-sample the TIF using gdalwarp of GDAL.

gdalwarp -tr 1 1 ^
         -r cubic ^
         greifensteine\dop20c_33352_5612.tif ^
         greifensteine\dop1m_33352_5612.tif

If you are not yet using GDAL today is a good day to start. It nicely complements the point cloud processing functionality of LAStools for raster inputs. Next we use lascolor to give each elevation pixel of the DTM stored in LAZ format its corresponding color from the orthophoto.

lascolor -i greifensteine\333525612_dgm1.laz ^
         -image greifensteine\dop1m_33352_5612.tif ^
         -odix _rgb -olaz

Now we can view the colored DTM in LAZ format interactively with lasview or any other LiDAR viewing software and turn on the RGB colors from the orthophoto as needed to understand the scene.

lasview -i greifensteine\333525612_dgm1_rgb.laz

We thank the “Staatsbetrieb Geobasisinformation und Vermessung Sachsen (GeoSN)” for giving us easy access to the 1 meter DTM and the 20 cm orthophoto that we have used in this article through their new open geodata portal as open data under the user-friendly license “Datenlizenz Deutschland – Namensnennung – Version 2.0.

National Open LiDAR Strategy of Latvia humiliates Germany, Austria, and other European “Closed Data” States

Latvia, officially the Republic of Latvia, is a country in the Baltic region of Northern Europe has around 2 million inhabitants, a territory of 65 thousand square kilometers and – since recently – also a fabulous open LiDAR policy. Here is a list of 65939 tiles in LAS format available for free download that cover the entire country with airborne LiDAR with a density from 4 to 6 pulses per square meters. The data is classified into ground, building, vegetation, water, low noise, and a few other classifications. It is licensed Creative Commons CC0 1.0 – meaning that you can copy, modify, and distribute the data, even for commercial purposes, all without asking permission. And there is a simple and  functional interactive download portal where you can easily download individual tiles.

latvia_open_data_portal_01

Interactive open LiDAR download portal of Latvia.

We downloaded the 5 by 5 block of square kilometer tiles matching “4311-32-XX.las” for checking the quality and creating a 1m DTM and a 1m DSM raster. You can follow along after downloading the latest version of LAStools.

Quality Checking

We first run lasvalidate and lasinfo on the downloaded LAS files and then immediately compress them with laszip because multi-core processing of uncompressed LAS files will quickly overwhelm our file system, make processing I/O bound, and result in overall longer processing times with CPUs waiting idly for data to be loaded from the drives.

lasinfo -i 00_tiles_raw\*.las ^
        -compute_density ^
        -histo z 5 ^
        -histo intensity 256 ^
        -histo user_data 1 ^
        -histo scan_angle 1 ^
        -histo point_source 1 ^
        -histo gps_time 10 ^
        -odir 01_quality -odix _info -otxt ^
        -cores 3
lasvalidate -i 00_tiles_raw\*.las ^
            -no_CRS_fail ^
            -o 01_quality\report.xml

Despite already excluding a missing Coordinate Reference System (CRS) from being a reason to fail (the lasinfo reports show that the downloaded LAS files do not have any geo-referencing information) lasvalidate still reports a few failing files, but scrutinizing the resulting XML file ‘report.xml’ shows only minor issues.

Usually during laszip compression we do not alter the contents of a file, but here we also add the EPSG code 3059 for CRS “LKS92 / Latvia TM” as we turn bulky LAS files into slim LAZ files so we don’t have to specify it in all future processing steps.

laszip -i 00_tiles_raw\*.las ^
       -epsg 3059 ^
       -cores 2

Compression reduces the total size of the 25 tiles from over 4.1 GB to below 0.6 GB.

Next we use lasgrid to visualize the last return density which corresponds to the pulse density of the LiDAR survey. We map each 2 by 2 meter pixel where the last return density is 2 or less to blue and each 2 by 2 meter pixel it is 8 or more to red.

lasgrid -i 00_tiles_raw\*.laz ^
        -keep_last ^
        -step 2 ^
        -density_16bit ^
        -false -set_min_max 2 8 ^
        -odir 01_quality -odix _d_2_8 -opng ^
        -cores 3

This we follow by the mandatory lasoverlap check for flight line overlap and alignment where we map the number of overlapping swaths as well as the worst vertical difference between overlapping swaths to a color that allows for quick visual quality checking.

lasoverlap -i 00_tiles_raw\*.laz ^
           -step 2 ^
           -min_diff 0.1 -max_diff 0.2 ^
           -odir 01_quality -opng ^
           -cores 3

The results of the quality checks with lasgrid and lasoverlap are shown below.

Raster Derivative Generation

Now we use first las2dem to create a Digital Terrain Model (DTM) and a Digital Surface Model (DSM) in RasterLAZ format and then use blast2dem to create merged and hill-shaded versions of both. Because we will use on-the-fly buffering to avoid edge effects along tile boundaries we first spatially index the data using lasindex for more efficient access to the points from neighboring tiles.

lasindex -i 00_tiles_raw\*.laz ^
         -cores 3

las2dem -i 00_tiles_raw\*.laz ^
        -keep_class 2 9 ^
        -buffered 25 ^
        -step 1 ^
        -use_orig_bb ^
        -odir Latvia\02_dtm_1m -olaz ^
        -cores 3

blast2dem -i 02_dtm_1m\*.laz ^
          -merged ^
          -hillshade ^
          -step 1 ^
          -o dtm_1m.png

las2dem -i 00_tiles_raw\*.laz ^
        -drop_class 1 7 ^
        -buffered 10 ^
        -spike_free 1.5 ^
        -step 1 ^
        -use_orig_bb ^
        -odir 03_dsm_1m -olaz ^
        -cores 3

blast2dem -i 03_dsm_1m\*.laz ^
          -merged ^
          -hillshade ^
          -step 1 ^
          -o dsm_1m.png

Because the overlaid imagery does not look as nice in our new Google Earth installation, below are the DTM and DSM at versions down-sampled to 25% of their original size.

Many thanks to SunGIS from Latvia who tweeted us about the Open LiDAR after we chatted about it during the Foss4G 2019 gala dinner. Kudos to the Latvian Geospatial Information Agency (LGIA) for implementing a modern national geospatial policy that created opportunity for maximal return of investment by opening the expensive tax-payer funded LiDAR data for re-purposing and innovation without barriers. Kudos!

LASmoons: Gabriele Garnero

Gabriele Garnero (recipient of three LASmoons)
Interuniversity Department of Regional and Urban Studies and Planning
Politecnico e Università degli Studi, Torino
ITALY

Background:
Last spring, the LARTU research group produced a laser scanner survey of the Abbey of Sant’Andrea in Vercelli, on the occasion of the VIII centenary of the dedication (1219). The database produced with a topographic tool that integrates the potential of a total station with laser scanner and photogrammetric sensors (Trimble SX 10), has been used to produce representations that can be consulted in interactive mode, navigating within the point clouds and producing a consultation platform that can also be accessed by non-specialist users such as art historians or archaeologists.

lasmoons_gabriele_garnero_0

Goal:
The LAStools software will be used to improve both the point cloud produced by eliminating the remaining noises, and check other ways of publishing the data, so as to make it usable from outside, to the community of researchers.

Data:
+
laser scanner and photogrammetric acquisitions of the interior of the building (150 millions of points)
+ laser scanner and photogrammetric acquisitions of the outside of the building (210 millions of points)
+ drone-based shooting of outdoor areas processed with Pix4D (23 millions of points)

LAStools processing:
1) tile large point cloud into tiles with buffer [lastile]
2) mark set of points whose z coordinate is a certain percentile of that of their neighbors [lasthin]
3) remove isolated low points from the set of marked points [lasnoise]
4) classify marked points into ground and non-ground [lasground]
5) creates a LiDAR portal for 3D visualization of LAS files [laspublish]

In Sweden, all they Wanted for Christmas was National LiDAR as Open Data

Let’s heat up some sweet, warm and spicy Glögg in celebration! They must have been good boys and girls up there in Sweden. Because “Jultomten” or simply ”Tomten” – how Sweden’s Santa Clause is called – is assuring a “God Jul” for all the Swedish LiDAR lovers this Christmas season.

Only a few weeks ago this tweet of ours had (mistakenly) included Sweden in a list of European countries that had released their national LiDAR archives as open data for public reuse over the past six years.

Turns out we were correct after all. Sweden has just opened their LiDAR data for free and unencumbered download. To get the data simply create a user account and browse to the ftp site for download as shown in the image sequence below.

The released LiDAR data was collected with a density of 1 to 2 pulses per square meter and is distributed in LASzip compressed LAZ tiles of 2500 by 2500 meters. The returns are classified into four classes: ground (2), water (9), low noise (7) and high noise (18). All items that can not be classified as any of the first four classes coded as left unclassified (1). The LAZ files do not contain CRS information, but this can easily be added with horizontal coordinates in SWERED99 TM (EPSG code 3006) and elevations in RH2000 height (EPSG code 5613).

Below a look with lasview at a 5 km by 5 km area that composed of the four tiles ‘18P001_67100_5800_25.laz‘, ‘18P001_67100_5825_25.laz‘, ‘18P001_67125_5800_25.laz‘ and ‘18P001_67125_5825_25.laz‘ with several of the different color modes available.

 

Some more details: The data was acquired at flying altitude of around 3000 meter with a maximum scan angle of ± 20º and a minimum side overlap of 10% between the flightlines. The laser footprint on ground is below 75 centimeters with slight variation based on the flying altitude. The laser scanning survey was performed with LiDAR instruments that can provide at least three returns from the same pulse. All LiDAR returns are preserved throughout the entire production chain.

The LiDAR data comes with the incredibly Creative Commons – CC0 license, which means that you can use, disseminate, modify and build on the data – even for commercial purposes – without any restrictions. You are free to acknowledge the source when you distribute the data further, but it is not required.

The LiDAR data will eventually cover approximately 75% of Sweden and new point clouds will continuously be added as additional scanning is performed according to the schedule shown below. The survey will be returning to scan every spot again after about 7 years.

2018-2022 LiDAR acquisition plan for Sweden

Below a lasinfo report for tile ‘18P001_67125_5825_25.laz‘. One noticeable oddity is the distribution of intensities. The histogram across all intensities with bins of size 256 shows two clearly distinct sets of intensities each with their own peak and a void of values between 3000 and 10000.

lasinfo -i 18P001_67125_5825_25.laz -cd -histo intensity 256
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            1
  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000
  version major.minor:        1.2
  system identifier:          ''
  generating software:        'TerraScan'
  file creation day/year:     303/2018
  header size:                227
  offset to point data:       227
  number var. length records: 0
  point data format:          1
  point data record length:   28
  number of point records:    20670652
  number of points by return: 13947228 4610837 1712043 358397 42147
  scale factor x y z:         0.01 0.01 0.01
  offset x y z:               0 0 0
  min x y z:                  582500.00 6712500.00 64.56
  max x y z:                  584999.99 6714999.99 136.59
LASzip compression (version 3.2r2 c2 50000): POINT10 2 GPSTIME11 2
reporting minimum and maximum for all LAS point record entries ...
  X            58250000   58499999
  Y           671250000  671499999
  Z                6456      13659
  intensity          32      61406
  return_number       1          5
  number_of_returns   1          5
  edge_of_flight_line 0          1
  scan_direction_flag 0          1
  classification      1         18
  scan_angle_rank   -19         19
  user_data           0          1
  point_source_ID  1802       1804
  gps_time 222241082.251248 222676871.876191
number of first returns:        13947228
number of intermediate returns: 2110980
number of last returns:         13952166
number of single returns:       9339722
covered area in square units/kilounits: 5923232/5.92
point density: all returns 3.49 last only 2.36 (per square units)
      spacing: all returns 0.54 last only 0.65 (in units)
overview over number of returns of given pulse: 9339722 5797676 4058773 1263967 210514 0 0
histogram of classification of points:
        10888520  unclassified (1)
         9620725  ground (2)
           22695  noise (7)
          138147  water (9)
             565  Reserved for ASPRS Definition (18)
intensity histogram with bin size 256.000000
  bin [0,256) has 1753205
  bin [256,512) has 3009640
  bin [512,768) has 2240861
  bin [768,1024) has 1970696
  bin [1024,1280) has 1610647
  bin [1280,1536) has 1285858
  bin [1536,1792) has 974475
  bin [1792,2048) has 790480
  bin [2048,2304) has 996926
  bin [2304,2560) has 892755
  bin [2560,2816) has 164142
  bin [2816,3072) has 57367
  bin [3072,3328) has 18
         [void]
  bin [10752,11008) has 589317
  bin [11008,11264) has 3760
  bin [11264,11520) has 99653
  bin [11520,11776) has 778739
  bin [11776,12032) has 1393569
  bin [12032,12288) has 1356850
  bin [12288,12544) has 533202
  bin [12544,12800) has 140223
  bin [12800,13056) has 16195
  bin [13056,13312) has 2319
  bin [13312,13568) has 977
  bin [13568,13824) has 765
  bin [13824,14080) has 648
  bin [14080,14336) has 289
  bin [14336,14592) has 513
  bin [14592,14848) has 383
  bin [14848,15104) has 178
  bin [15104,15360) has 526
  bin [15360,15616) has 108
  bin [15616,15872) has 263
  bin [15872,16128) has 289
  bin [16128,16384) has 69
  bin [16384,16640) has 390
  bin [16640,16896) has 51
  bin [16896,17152) has 186
  bin [17152,17408) has 239
  bin [17408,17664) has 169
  bin [17664,17920) has 58
  bin [17920,18176) has 227
  bin [18176,18432) has 169
  bin [18432,18688) has 40
  bin [18688,18944) has 401
  bin [18944,19200) has 30
  bin [19200,19456) has 411
  bin [19456,19712) has 34
  bin [19712,19968) has 34
  bin [19968,20224) has 398
  bin [20224,20480) has 24
  bin [20480,20736) has 108
  bin [20736,20992) has 267
  bin [20992,21248) has 29
  bin [21248,21504) has 318
  bin [21504,21760) has 26
  bin [21760,22016) has 59
  bin [22016,22272) has 184
  bin [22272,22528) has 52
  bin [22528,22784) has 18
  bin [22784,23040) has 116
  bin [23040,23296) has 55
  bin [23296,23552) has 89
  bin [23552,23808) has 250
  bin [23808,24064) has 24
  bin [24064,24320) has 52
  bin [24320,24576) has 14
  bin [24576,24832) has 29
  bin [24832,25088) has 71
  bin [25088,25344) has 74
  bin [25344,25600) has 2
  bin [25600,25856) has 17
  bin [25856,26112) has 2
  bin [26368,26624) has 9
  bin [26624,26880) has 1
  bin [26880,27136) has 1
  bin [27136,27392) has 1
  bin [27392,27648) has 1
  bin [27648,27904) has 3
  bin [28416,28672) has 2
  bin [29184,29440) has 4
  bin [30720,30976) has 1
  bin [30976,31232) has 2
  bin [31232,31488) has 1
  bin [32512,32768) has 1
  bin [36864,37120) has 1
  bin [58368,58624) has 1
  bin [61184,61440) has 1
  average intensity 3625.2240208968733 for 20670652 element(s)