We often process LiDAR in tiles for two reasons: first, to keep the number of points per file low and use main memory efficient, and second, to speed up the computation with parallel tile processing and keep all cores of a modern CPU busy. However, it is very (!!!) important to take the necessary precautions to avoid “edge artifacts” when processing LiDAR in tiles. We have to include points from neighboring tiles during certain LAStools processing steps to avoid edge artifacts. Why? Here is an illustration from our PHIL LiDAR tour earlier this year:
What you see is the temporary TIN of ground points created (internally) by las2dem or blast2dem that is then rastered at the user-specified step size onto a grid. Without a buffer (right side) there will not always be a triangle to cover every pixel. Especially in the corners of the tile you will often find empty pixels. Furthermore the poorly shaped “sliver triangles” along the boundary of the TIN do not interpolate the ground elevations properly. In contrast, with buffer (left side) the TIN generously covers the entire area that is to be rastered with nicely shaped triangles.
Here the christmas cookies analogy: You need to roll out the dough larger than the cookies you will cut to make sure your cookies will have nice edges. Think of the TIN as the dough and the square tile as your cookie cutter. You need to use a sufficiently large buffer when you roll out your TIN to assure an edge without crumbles when you cut out the tile … (-: … otherwise you are pretty much guaranteed to get results that – upon closer inspection – have these kind of artifacts:
Without buffers processing artifacts also happen when classifying points with lasground or lasclassify, when calculating height above ground or height-normalizing LiDAR tiles with lasheight, when removing noise with lasnoise, when creationg contours with las2iso or blast2iso, or any other operation where an incomplete neighborhood of points can affect the results. Hence, we need to surround each tile with a temporary buffer of points. Currently there are two ways of working with buffers with LAStools:
- creating buffered tiles during the initial tiling step with the ‘-buffer 25’ option of lastile, maintaining buffered tiles throughout processing and finally using the ‘-use_tile_bb’ option of lasgrid, las2dem, blast2dem, or lascanopy to raster the tiles without the temporary buffer.
- creating buffered tiles from non-overlapping (= unbuffered) tiles with “on-the-fly” buffering using the ‘-buffered 25’ option of most LAStools such as lasground, lasheight, or las2dem. For some workflows it is useful to also add ‘-remain_buffered’ if buffers are needed again in the next step. Finally, we use the ‘-use_orig_bb’ option of lasgrid, las2dem, blast2dem, or lascanopy to raster the tiles without the temporary buffer.
In the following three (tiny) examples using the venerable ‘fusa.laz’ sample that is distributed with LAStools to illustrate the two types of buffering as well as to show what happens when no buffers are used. In each example we will first cut the small ‘fusa.laz’ sample into nine smaller tiles and then process these separately on 4 cores in parallel.
1. Initial buffer creation with lastile
This is what most of my tutorials teach. It assumes you are the one creating the tiling in the first place. If you do it with lastile and add a buffer right from the start things are pretty easy.
lastile -i ..\data\fusa.laz ^ -set_classification 0 -set_user_data 0 ^ -tile_size 100 -buffer 20 ^ -odir 1_raw -o futi.laz
We cut the input into 100 meter by 100 meter tiles but add a 20 meter buffer around each tile. That means that each tile on disk will contain the points for an area of up to 140 meter by 140 meter. The GUI for LAStools shows the overlap and if you scrutinize the bounding box values that the cursor points to you notice the extra 20 meters in each direction.
lasground -i 1_raw\futi*.laz ^ -city ^ -odir 1_ground -olaz ^ -cores 4 lasheight -i 1_ground\futi*.laz ^ -drop_above 50 ^ -odir 1_height -olaz ^ -cores 4 lasclassify -i 1_height\futi*.laz ^ -odir 1_classify -olaz ^ -cores 4
las2dem -i 1_classify\futi*.laz ^ -keep_class 2 6 ^ -step 0.25 -use_tile_bb ^ -odir 1_dbm -obil ^ -cores 4
We created a digital terrain model with buildings (DBM) by keeping the points with classification 2 (ground) and 6 (building). After loading the resulting 9 tiles into QGIS and generating a virtual raster we see a nice seamless DBM without any edge artifacts.
If you need to deliver the LiDAR files you should remove the buffers with lastile and option ‘-remove_buffer’.
lastile -i 1_classify\futi*.laz ^ -remove_buffer ^ -odir 1_final -olaz ^ -cores 4
2. On-the-fly buffering
Now assume you are given LiDAR tiles without buffers. We generate them here with lastile.
lastile -i ..\data\fusa.laz ^ -set_classification 0 -set_user_data 0 ^ -tile_size 100 ^ -odir 2_raw -o futi.laz
The only difference is that we do not request the 20 meter buffer and the result is a typical tiling as you may receive it from a vendor or download it from a LiDAR portal. The GUI for LAStools shows that there is no overlap and if you scrutinize the bounding box values that the cursor points to, you see that the tiles is exactly 100 meters bty 100 meters.
Now we have to think about buffers a lot. When using on-the-fly buffering we should first spatially index the tiles with lasindex for faster access to the points from neighbouring tiles.
lasindex -i 1_raw\futi*.laz -cores 4
Below in red are the modifications for on-the-fly buffering to the standard workflow of lasground, lasheight, and lasclassify. The first lasground run uses ‘-buffered 20’ to add buffers to each tile and ‘-remain_buffered’ to write those buffers to disk. This way they do not have to created again by lasheight and lasclassify.
lasground -i 2_raw\futi*.laz ^ -buffered 20 -remain_buffered ^ -city ^ -odir 2_ground -olaz ^ -cores 4 lasheight -i 2_ground\futi*.laz ^ -remain_buffered ^ -drop_above 50 ^ -odir 2_height -olaz ^ -cores 4 lasclassify -i 2_height\futi*.laz ^ -remain_buffered ^ -odir 2_classify -olaz ^ -cores 4
At the end we have to remember that the tiles still have on-the-fly buffers and them cut off with option ‘-use_orig_bb’ of las2dem.
las2dem -i 2_classify\futi*.laz ^ -keep_class 2 6 ^ -step 0.25 -use_orig_bb ^ -odir 2_dbm -obil ^ -cores 4
Again, we created a digital terrain model with buildings (DBM) by keeping the points with classification 2 (ground) and 6 (building). The resulting hillshade computed from a virtual raster that combines the 9 BIL rastera into one looks perfectly smooth in QGIS.
If you need to deliver the LiDAR files you should probably remove the buffers first … but that is not yet implemented. (-:
lastile -i 2_classify\futi*.laz ^ -remove_buffer ^ -odir 2_final -olaz ^ -cores 4
3. Bad: No buffering
Here what you are *not* supposed to do. Assuming you get unbuffered tiles.
lastile -i ..\data\fusa.laz ^ -set_classification 0 -set_user_data 0 ^ -tile_size 100 ^ -odir 3_raw -o futi.laz
Bad. You do not take care about buffering when processing the tiles.
lasground -i 3_raw\futi*.laz ^ -city ^ -odir 3_ground -olaz ^ -cores 4 lasheight -i 3_ground\futi*.laz ^ -drop_above 50 ^ -odir 3_height -olaz ^ -cores 4 lasclassify -i 3_height\futi*.laz ^ -odir 3_classify -olaz ^ -cores 4
Bad. You do not take care about buffering when generating the DBM.
las2dem -i 3_classify\futi*.laz ^ -keep_class 2 6 ^ -step 0.25 ^ -odir 3_dbm -obil ^ -cores 4
Bad. You get crappy results with edge artifacts clearly visible in the hillshade.
Bad. If you zoom in on a corner where 4 tiles meet you find missing pixels and incorrect elevation values. Bad. Bad. Bad. So please folks. Try this on your own data. Notice the horrible edge artifacts. Then always use buffers … (-: