Pre-Processing Mobile Rail LiDAR with LAStools

The majority of LAStools users are processing airborne LiDAR. That should not surprise as airborne is by far the most common form of LiDAR in terms of square kilometers covered. The availability of LiDAR as “open data” is also pretty much restricted to airborne surveys, which are often tax-payer funded and then distributed freely to achieve maximum return of investment.

But folks are increasingly using our software to do some of the “heavy lifting” for mobile LiDAR, either mounted on a truck for scanning cities or on a train for capturing railroad infrastructure. The LiDAR collected for the cities of Budapest and Singapore, for example, was pre-processed by multi-core scripted LAStools when the scanning trucks returned with their daily trajectories worth of point clouds captured by a RIEGL VMX-450 mobile mapping system.

One customer who was recently scanning railroad infrastructure wanted to do automatic ground classification as a first step prior to further segmentation of the data. We were asked for advice because on such data the standard settings of lasground left too many patches of ground unclassified. Also the uniform tiling lastile generates by default is not a good way to break such data into manageable pieces given the drastically varying point densities in mobile scanning.

We obtained a 217 MB file in LAZ format with 40 million points corresponding to a 2.7 km stretch of railway track. We first run a quick lasindex (with the options for ‘mobile’) on the file that creates a spatial indexing LAX file with maximally 10 meter resolution. This not only allows faster area-of-interest queries but also gives us a more detailed preview than just the bounding box of where the LiDAR points actually are in the GUI of LAStools.

mobile_rail_lidar_01

Presence of LAX files results in actual extend of LiDAR being shown in GUI.

lasindex -i segment.laz -tile_size 10 -maximum -100

We then run lastile four times to create an adaptive tiling in which no tile has more than 6 million points. The first call creates the initial 1000 by 1000 meter tiles. The following three calls refine all those tiles that still have more than 6 million points first into 500 by 500 meter, then 250 by 250 meter, and finally 125 by 125 meter tiles in parallel on 4 cores. Note the ‘-refine_tiling’ option is used in the first call to lastile and the ‘-refine_tiles’ option in all subsequent calls.

lastile -i segment.laz ^
        -tile_size 1000 ^
        -buffer 10 -flag_as_withheld ^
        -refine_tiling 6000000 ^
        -odir tiles_raw -o rail.laz
lastile -i tiles_raw\*_1000.laz ^
        -flag_as_withheld ^
        -refine_tiles 6000000 ^
        -olaz ^
        -cores 4
lastile -i tiles_raw\*_500.laz ^
        -flag_as_withheld ^
        -refine_tiles 6000000 ^
        -olaz ^
        -cores 4
lastile -i tiles_raw\*_250.laz ^
        -flag_as_withheld ^
        -refine_tiles 6000000 ^
        -olaz ^
        -cores 4

The resulting tiles all have fewer than 6 million points but still have the initial 10 meter buffer that was specified by the first call to lastile. Two tiles were sufficiently small after the 1st call, three tiles after the 2nd call, eleven tiles after 3rd call, and three tiles after the 4th.

contents of tile shown in blue in adaptive tiling below

points of adaptive tile (high-lighted in blue below) colored by intensity

Adaptive tiling created with four calls to lastile.

Adaptive tiling created with four calls to lastile. Scale factors of 0.00025 (see mouse cursor) implies that point coordinates are stored with quarter millimeter resolution. Lowering them to 0.001 would result in better compression and lower I/O.

Noise in the data – especially low noise – can lead lasground into choosing the wrong points during ground classification by latching on to those low noise points. We first classify the noise points into a different class (7) using lasnoise so we can later ignore them. These particular settings were found by experimenting on a few tiles with different values (see the README file) until visual inspection showed that most low points had been classified as noise.

lasnoise -i tiles_raw\*.laz ^
         -step_xy 0.5 -step_z 0.1 ^
         -odir tiles_denoised -olaz ^
         -cores 4
noise points shown in violett

noise points shown in violett

The points classified as noise will not be considered as ground points during the next step. For this it matters little that lamp posts, wires, or vegetation are wrongly marked as noise now. We can always undo their noise classification once the ground points were classified. Important is that those pointed to by the mouse cursor, which are below the desired ground, are excluded from consideration during the ground classification step. Here those low points are not actually noise but returns generated wherever the laser was able to “peek” through an opening to a lower surface.

lasground -i tiles_denoised\*.laz ^
          -ignore_class 7 ^
          -step 1 -sub 3 -bulge 0.1 -spike 0.1 -offset 0.02 ^
          -odir tiles_ground -olaz ^
          -cores 4

For classification with lasground there are a number of options to play with  (see the README file) but the most important is the correct step size. It is terrain along the railway track bed that is supposed to get represented well. The usual step of 5 to 40 meter for lasground aim at the removal of vegetation and man-made structures from airborne LiDAR. They are not the right choice here. A step of 1 and the parameters shown above gives us the ground shown below.

Classification of terrain along railway track using lasground with '-step 1'

Classification of terrain along railway track bed using lasground with ‘-step 1’

The new ‘-flag_as_withheld’ option in lastile that flags each point in the buffer with the withheld flag is useful in case we want to remove all buffer points on-the-fly, for example, in order to create a DTM hillshade of 25 cm resolution for a visual quality check of the entire 2.7 km track using blast2dem from the BLAST extension of LAStools.

blast2dem -i tiles_ground\*.laz -merged ^
          -drop_withheld -keep_class 2 ^
          -hillshade -step 0.25 ^
          -o dtm_hillshaded.png
Small 600 x 600 pixel detail of hill-shaded 5663 x 9619 pixel DTM raster generated by blast2dem

Small 600 x 600 pixel detail of hill-shaded 5663 x 9619 pixel DTM raster generated by blast2dem.

Worldwide LiDAR of Rainforest Biomass for REDD+

We at rapidlasso have long been big fans of the biomass and biodiversity work done by Greg Asner’s group and their Carnegie Airborne Observatory. They were, in fact, the earliest “power-users” of the BLAST extension of LAStools and helped finding all the bugs when rasterizing billions of LiDAR points collected during large-scale surveys in the Amazon into Digital Elevation Models (DEMs) with blast2dem. Below a video fly-through of some of the LiDAR they collected.

A few days ago, Greg Asner together with his colleagues Joseph Mascaro, Stuart Davies, Alex Dehgan, and Sassan Saatchi published a thought-provoking article called “These are the days of lasers in the jungle” which is essentially a “call for action” to map the world’s tropical forests with a fleet of airplanes outfitted with advanced LiDAR to rapidly and to accurately assess global forest carbon stocks.

Why would anyone want to do this? In order to properly quantify actual emissions savings for REDD+ (Reduced Emissions from Deforestation and Degradation). REDD+ is a tropical forest carbon accounting program of the United Nations Framework Convention on Climate Change that aims to compensate tropical countries for reducing carbon emissions from deforestation and forest degradation that account for roughly 10 percent of global greenhouse gas emissions. The key to implement such a program is the ability to accurately and affordably estimate the actual emissions savings and a worldwide LiDAR inventory of tropical forests will accomplish just that argues the paper. This “call for action” has since been picked up by Mongabay – a site that examines emerging trends in climate, technology, and finance on conservation and development.

Interesting is the price tag that they estimate, which is a fraction of the cost of a typical Earth observation satellite mission, They claim: “Our ambitious plan can be accomplished for far less than what we have already spent on avoided deforestation. Aircraft leasing, data collection and processing costs for 30 days of flying can reasonably be limited to USD 500,000 Using this monthly sampling unit, collecting at an average of 100,000 hectares per day, a fleet of ten aircraft could do the job in four years at USD 250 million, or just 5% of pledged REDD + funding.”

The authors state “The time has come for a brute-force effort to directly assess the carbon stock for all of the world’s tropical forests by 2020” because “airborne LiDAR is uniquely suited for this role because it can be collected, standardized, reported and verified in a simple manner by both a landholder and any third party”. Should such a campaign turn out to be a viable option to implement REDD+ we hope full waveform LiDAR – not just discrete returns – will be collected.

Waveform digitizers have become popular because they can capture the reflection of the emitted laser pulse with much more detail than a discrete return system. The intensity of the signal returning to the plane is digitized up to one billion times per second, giving a vertical resolution of one digitized amplitude each 15 centimeters. As this can capture the interactions between each laser pulse and the vegetation with much greater detail, it would seem that having full waveform instead of discrete returns will prove especially useful for biomass studies. A future blog post will talk about our own experiences of scanning tropical rainforest to produce full waveform LiDAR in PulseWaves format.

Another approach with a similar objective is being taken by NASA with their future GEDI space LiDAR. The Global Ecosystem Dynamics Investigation (GEDI) instrument will be the first to systematically probe the depths of the forests from space to reveal their 3D structure, as depicted in the artist’s concept below, and provide crucial information about the impact that trees have on the amount of carbon in the atmosphere. “One of the most poorly quantified components of the carbon cycle is the net balance between forest disturbance and regrowth”, said Ralph Dubayah, the GEDI principal investigator at the University of Maryland. “GEDI will help scientists fill in this missing piece by revealing the vertical structure of the forest, which is information we really can’t get with sufficient accuracy any other way”. The instrument will be built at NASA’s Goddard Space Flight Center in Greenbelt, Maryland. Not sure about the price tag.

nasa_gedi

Image Credit: NASA’s Goddard Space Flight Center

Until either of these projects become reality, “enjoy” this video by Carnegie Airborne Observatory of a 3D flyover that shows rapidly expanding palm oil plantations in the Peruvian Amazon rainforest that are contributing to deforestation.

locating German bunkers concealed by canopy

After accidentally finding Russian tanks in Polish forests I was curious to see if there was something else hiding under the forest canopy. Remember, I randomly picked a 500 by 500 meter LiDAR tile as example data to introduce a group of forestry students to LiDAR processing with LAStools during the ForseenPOMERANIA camp. After extracting ground points with “lasground.exe“, strange bumps appeared in the bare-earth hillshades generated with “las2dem.exe” for terrain that was supposed to be completely flat … they turned out to be Russian WW-II positions.

I met Achim when returning to teach the next two groups of students. His hobby is a mix between geo-caching and conflict-archaeology: locating old German bunkers based on approximate coordinates available in historic records and tourist maps and then mapping them precisely with GPS. Achim had a list of longitude/latitude positions as KML files where he was planning to search for known bunkers. I used “lasboundary.exe” to create polygonal outlines in KML format for all areas where we had LiDAR from the forestry project. With Google Earth it was easy to find overlaps between his target areas and our LiDAR coverage.

I extracted the ground points and created bare-earth DTMs of the relevant area with a LAStools batch processing pipeline of “lastile.exe“, “lasground.exe” and “las2dem.exe” and used “blast2dem.exe” to create a seamless hillshading with proper KML referencing (here is a tutorial for such a pipeline). What I found was pretty amazing.

At first glance it looks like a maze of little creeks that are running alongside the ridges of the hillsides but we know that water flows downhill and not “along-hill”. What we see is a network of WW-II trenches that are connecting the bunkers Achim is looking for. A closer look also reveals the likely location of those bunkers.

I placed a pin on each of them and exported their longitude and latitude coordinates for upload into Achim’s GPS device. The next day we set out to verify our LiDAR findings on the ground.

It was a rainy day. Walking through this maze of green and overgrown trenches from one moss-covered bunker ruin to the next felt oddly quiet and peaceful. Achim explained that these bunkers were originally built to defend the border with Poland – long before the Second World War broke out. Only when Russian soldiers were advancing on Germany after the collapse of the Eastern front, the young boys of the Hitler Youth were commanded to dig this network of trenches in order to fortify the bunker and stop enemy lines from gaining ground. After the war the Russians blew up all bunkers that were facing East so that their troops would not ever have to face them again.

Tutorial: derivative production

This is the final part of a three-part tutorial on how to use LAStools to implement a pipeline that (1) quality checks a newly received set of raw LiDAR flight strips, (2) tiles and prepares the LiDAR for subsequent exploitation, and (3) generates raster and vector derivatives from the LiDAR points such as DTMs, DSMs, contour lines and building footprints with multi-core batch processing and Tin streaming with BLAST.

It is maybe helpful to start with the tutorial on quality checking. But it is mandatory to first complete the tutorial on LiDAR preparation because that is where you generate the cleaned and classified LiDAR tiles in folders ’tiles_classified’ and ’tiles_final’ that we are using as input in this tutorial.

Generating matching DTM tiles without edge artifacts

Create a new folder called ‘.\lastools\bin\tiles_dtm’ for storing the bare-earth digital terrain models (DTMs). Then run ‘las2dem.exe’ in GUI mode and load the 4 tiles from the folder ’tiles_classified’. You can either do this by double-clicking ‘las2dem.exe’ and using the ‘browse …’ menu or by entering the command below:

C:\lastools\bin>las2dem -i tiles_classified\fitch*.laz -gui

The las2dem tool creates rasters from LiDAR by triangulating the points and then sampling the resulting TIN at the center of each raster pixel. Remember, the tiles in folder ’tiles_classified’ have a buffer of 50 meter that was initially created with ‘lastile’ and that has been carried through ‘lasground’, ‘lasheight’, and ‘lasclassify’. The presence of these buffers allows us now to avoid artifacts along the tile boundaries because our TINs will be 50 meters wider in all directions than the raster that we are sampling. Rasterizing only the extend of the original tile without the surrounding buffer is accomplished with the ’tile_bb_only’ command line parameter.

Set the GUI options as shown below to create one DTM per tile. The DTMs will seamlessly fit together. Check the button ‘use tile bounding box’ for the reason explained above and check the button ‘extra pass’ which makes two passes over the data. In the first pass it only counts the number of points that pass through the filter so that in the second it can optimize the memory used to triangulate them. Choose ‘bil’ as the output format. It is much much much more efficient than the ‘asc’ format many people like to use. Select 4 cores if you have that many and set the output directory to ’tiles_dtm’. Use the ‘filter…’ menu to choose two filters: the first keeps only the ground points with filter ‘keep_classification 2’ and the second keeps only one point every 0.5 by 0.5 meter area with the filter ‘thin_with_grid 0.5’. The order is really important (excercise: explain why!). Because we create DTMs with the default ‘step’ size of 1 meter there is no sense in constructing a temporary TIN that is more detailed than one ground point per 0.5 by 0.5 meter area. Constructing a large TIN requires a lot of main memory. We should only use as many points as necessary to sample the TIN at our target resolution ‘step’. Adding a ‘thin_with_grid’ filter with half the ‘step’ size seems always a good choice for ‘las2dem’.

tutorial3 las2dem GUI generate DTM

Press ‘RUN’ and the ‘RUN’ window will pop up. You may have to close the ‘output’ menu before you can actually see the ‘RUN’ button. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Alternatively, if you prefer to work in the command-line, you can use this equivalent command here.

C:\lastools\bin>las2dem -i tiles_classified\fitch*.laz ^
                        -keep_class 2 -thin_with_grid 0.5 ^
                        -extra_pass ^
                        -use_tile_bb ^
                        -odir tiles_dtm -obil ^
                        -cores 4

Generating matching DSM tiles without edge artifacts

Create a new folder called ‘.\lastools\bin\tiles_dsm’ for storing the first-return digital surface models (DSMs). You could re-use ‘las2dem’ if you have not closed the application yet. Otherwise re-run it in GUI mode and load the 4 tiles from the folder ’tiles_classified’ again. Do this by double-clicking ‘las2dem.exe’ and using the ‘browse …’ menu or by entering the command below:

C:\lastools\bin>las2dem -i tiles_classified\fitch*.laz -gui

Set the GUI options as shown below. The only changes are that the output directory name is now ’tiles_dsm’ and that the filtering is now keeping first returns instead of ground points. If you re-used the GUI from the last step you need to delete the old filters and enter the new ones in the correct order (exercise: why is the order important?). You can delete old filter entries by double clicking them. Select 4 cores only if you have a computer with that many cores. This will assign tiles to different cores and process them in parallel.

tutorial3 las2dem GUI dsm

Press ‘RUN’ and the ‘RUN’ window will pop up. You may have to close the ‘output’ menu before you can actually see the ‘RUN’ button. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Alternatively, if you prefer to work in the command-line, you can use this equivalent command here:

C:\lastools\bin>las2dem -i tiles_classified\fitch*.laz ^
                        -first_only -thin_with_grid 0.5 ^
                        -extra_pass ^
                        -use_tile_bb ^
                        -odir tiles_dsm -obil ^
                        -cores 4

Visualizing and processing DEM rasters as if they were points

Here a little trick. All LAStools can read the BIL format “as if” it was a point format. The other two raster formats that are supported this way are ESRI’s ASC and FUSION’s DTM. But if possible don’t use ASC because it stores the values as ASCII text which is terribly inefficient and much much slower to read than BIL or DTM. Run lasview in GUI mode and load the four BIL rasters from the folder ’tiles_dsm’. You need to check the button labelled ‘.bil’ in the ‘browse …’ menu to be able to select the BIL rasters as input. Alternatively you can use this command-line:

C:\lastools\bin>lasview -i tiles_dsm\fitch*.bil -gui

Choose the option ‘selected file only’ and pick one of the four tiles. Play with the ‘color_by’ or ‘render_only’ options. But these can also be changed once ‘lasview’ is displaying the point cloud.

tutorial3 lasview gui bil

Visualize the ‘fitch_273920_4714320.bil’ from the ’tiles_dsm’ folder. After all points have loaded press ‘t’ on the keyboard to triangulate them. Press ‘-‘ to make the points disappear and press ‘h’ multiple times to iterate over the possible ways to shade the triangulation.

tutorial3 lasview bil dsm

Using lasgrid to merge DEM tiles into one large DEM raster

Run ‘lasgrid.exe’ in GUI mode and load the 4 raster DSM tiles in BIL format from the ’tiles_dsm’ folder by double-clicking ‘lasgrid.exe’ and using the ‘browse …’ menu. Check the button labelled ‘.bil’ to be able to select the BIL rasters as input. Or use the following in the DOS command-line:

C:\lastools\bin>lasgrid -i tiles_dsm\fitch*.bil -gui

Configure the GUI with the settings shown below. Check the ‘merge files into one’ button, specify the output file as ‘dsm.tif’, and specify the northern UTM zone 19 and NAVD88 as the projection information because the BIL files we generated earlier do (currently) not store it.

tutorial3 lasgrid gui dsm merge

Press ‘RUN’. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Or use this in the DOS command-line:

C:\lastools\bin>lasgrid -i tiles_dsm\fitch*.bil -merged ^
                        -utm 19N -vertical_navd88 ^
                        -o dsm.tif

Load the merged raster ‘dsm.tif’ into a GIS package such as QGIS or ArcGIS and repeat the exercise by creating – in the exact same way – a merged ‘dtm.tif’ from the 4 raster DTM tiles in the ’tiles_dtm’ folder.

Using the BLAST extension for seamless generation of rasters

Run ‘blast2dem.exe’ in GUI mode and load the 4 raster DTM tiles in BIL format by double-clicking ‘blast2dem.exe’ and using the ‘browse …’ menu.  You need to check the button labelled ‘.bil’ to be able to select the BIL rasters as input. Or enter this in the DOS command-line window:

C:\lastools\bin>blast2dem -i tiles_dtm\fitch*.bil -gui

Configure the GUI with the settings shown below. Check the ‘merge files into one’ button, specify the output file as ‘dtm_hillshade_raster.png’, select ‘hillside shading’, and choose ‘png’ as the output format. Use the ‘projection …’ menu to specify the northern UTM zone 19 and NAVD88 as the projection information because the BIL files we generated earlier do (currently) not store it.

tutorial3 blast2dem GUI dtm merge hillshade

Press ‘RUN’. If you forgot or need to fix a setting you can directly modify the shown command-line before you press START. Or use this in the DOS command-line:

C:\lastools\bin>blast2dem -i tiles_dtm\fitch*.bil -merged ^
                          -hillshade ^
                          -utm 19N -vertical_navd88 ^
                          -o dtm_hillshade_raster.png

Whenever ‘blast2dem’, ‘las2dem’, ‘lasgrid’, or ‘lasoverlap’ generate rasters that are standard color or gray-scale PNG, TIF, or JPG that can be loaded as an image overlay into Google Earth, a KML is automatically generated – assuming there is projection information either available in the input or specified in the command-line. If you double-click the KML file that we just created you will see the hillshading overlaid and nicely aligned with Fitchburg airport.

tutorial3 blast2dem dtm merge hillshade

Repeat the exercise by creating – in the exact same way – a merged and hillshaded ‘dsm_hillshade_raster.png’ with KML file from the 4 raster DSM tiles in the ’tiles_dsm’ folder.

Using the BLAST extension for seamless generation of contours

Run ‘blast2iso.exe’ in GUI mode and load the 4 BIL raster tiles from the ’tiles_dtm’ folder by double-clicking and using the ‘browse …’ menu or by entering the command below:

C:\lastools\bin>blast2iso -i tiles_dtm\fitch*.bil -gui

Set the GUI as shown below to compute 1 meter contours from on-the-fly merged bare-earth DTM rasters in BIL format. Check the ‘merge files into one’ button, specify the output file name as ‘dtm_contours_raster_1m.shp’, set the isovalue extraction to every 1 unit (here: meter), and check the buttons for ‘simplify short segments’, ‘simplify small bumps’ and, ‘clean short lines’. Read the corresponding README.txt file if you want more details on what these options exactly do. Add the northern UTM zone 19 as the projection information because the BIL files we generated earlier do (currently) not store it. This is less important when generating SHP output but would be required to produce KML instead.

tutorial3 blast2iso GUI bil merged contours 1m

Press ‘RUN’. If you forgot or need to fix a setting you can directly modify the command-line before you press START. Or use this command-line:

C:\lastools\bin>blast2iso -i tiles_dtm\fitch*.bil -merged ^
                          -iso_every 1 ^
                          -simplify_length 0.5 -simplify_area 1 -clean 5 ^
                          -utm 19N ^
                          -o dtm_contours_raster_1m.shp

Visualize the result with the convenient ShapeViewer.exe program and check if the 1 meter contours that ‘blast2dem’ has produced are what you expected:

tutorial3 blast2iso bil merged 1m contours

Repeat the exercise by creating – in the exact same way – a merged and hillshaded ‘dsm_contours_raster_3m.shp’ with from the 4 raster DSM tiles in the ’tiles_dsm’ folder but with a spacing of 3 meters instead of 1 meter between the contour lines.

Extracting building footprints from on-the-fly merged tiles

Run ‘lasboundary.exe’ in GUI mode and load the 4 classified tiles with:

C:\lastools\bin>lasboundary -i tiles_final\fitch*.laz -gui

Set the GUI as shown below to generate compute polygonal outlines around the buildings points classified in the previous tutorial. Check the ‘merge files into one’ button to create one seamless output file. Add a ‘keep_classification 6’ filter that filters out all points except those classified as building. Check the ‘disjoint’ button – otherwise a funny-looking single polygonal hull will be created connecting all buildings. Set the ‘concavity’ to 1.5 meters to specify the smallest feature that lasboundary can grow into. This value needs to be 2 to 3 times the average point spacing. If you set it too small, lasboundary will “eat” its way into the buildings. If you set it too high nearby buildings will be joined together and the outlines will become coarser. Play with a few different settings such as 1.0, 2.0, 5.0, and 10.0 to get an intuitive idea how this parameter affects the result.

tutorial3 lasboundary GUI merged building footprints

Press ‘RUN’. If you forgot or need to fix a setting you can directly modify the command-line before you press START. Or use the equivalent command-line below directly in the DOS window.

C:\lastools\bin>lasboundary -i tiles_final\fitch*.laz -merged ^
                            -keep_class 6 ^
                            -concavity 1.5 -disjoint ^
                            -o buildings.shp

Visualize the result with the convenient ShapeViewer.exe program and check if the building footprints that lasboundary has produced are what you expected:

tutorial3 lasboundary merged building footprints

By generating KML output instead and overlaying the result on Google Earth imagery you can easily check which and how many buildings were missed and go back to adjust the parameters for ‘lasclassify’ in the previous tutorial to maybe achieve better results.

C:\lastools\bin>lasboundary -i tiles_final\fitch*.laz -merged ^
                            -keep_class 6 ^
                            -concavity 1.5 -disjoint ^
                            -o buildings.kml

tutorial3 lasboundary merged buildings kml

LAStools’ BLAST extension can process billions of LiDAR points

Often I get asked about the difference between las2dem.exe and blastdem.exe – the latter being part of the BLAST extension of LAStools. In the academic video from 2006 below, I outline the core technology details behind BLAST, which can seamlessly compute gigantic Delaunay Triangulations from billions of LiDAR points – using very little main memory. The blast2dem.exe module does not store the gigantic TINs that BLAST produces, but immediately rasters the finalized triangles streaming out of BLAST into equally massive DEMs. This avoids to ever having to store such huge triangulations – that are only needed temporarily – to disk. Currently blast2dem.exe can “only” process up to 2 billion points per file – a limitation that shall be lifted soon.