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.
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)
Pingback: Digitales Oberflächenmodell von Hamburg | Hannes ihm sein Blog