LASmoons: Geoffrey Ower

Geoffrey Ower (recipient of three LASmoons)
School of Biological Sciences
Illinois State University, Normal, USA

Background:
The spatial distribution and abundance of mosquitoes is important because biting mosquitoes can acquire and transmit pathogens such as viruses that infect humans or domestic animals. There have been 1,263 confirmed human infections with mosquito-borne West Nile virus in Illinois since 2003. In 2004 there was an epidemic of mosquito-borne La Crosse encephalitis, which resulted in 9 human cases in Illinois (USGS, HHS & CDC 2015). Understanding the geographical factors that influence mosquito species’ distributions could help to predict at-risk areas, which could in turn help to prioritize public health efforts to control mosquitoes and the diseases that they transmit more effectively.

lasmoons_Marta_Grech_0

The spatial distribution of mosquitoes depends in a large part upon the availability of water sources because mosquito larval and pupal stages are aquatic. Some mosquito species are specialized to use water-filled natural (e.g., treeholes, rockholes) or artificial (e.g., buckets, discarded tires) containers. Mosquitoes also require access to blood meals, access to carbohydrates (e.g., nectar), and resting sites. Spatial factors thought to be important in determining the spatial distribution and abundance of container-dwelling mosquitoes include land cover and land use (Diuk-Wasser et al. 2006), temperature, precipitation (Ruiz et al. 2010), elevation (Sun et al. 2009), human population density (Higa et al. 2010), and socioeconomic status (Dowling et al. 2013).

Goal:

The objective of this project is to determine what spatial factors predict the distribution and abundance of mosquito species in Bloomington-Normal, Illinois. Species distribution maps will be produced for each species of mosquito that colonized oviposition traps (water-filled plastic cups lined with paper on which mosquito eggs are laid) placed on sampling transects during three sampling periods in August and September 2015. Poisson regression models will be used to produce maps predicting the occurrence of each mosquito species for the full 509 square kilometre study area.

Data:
+
509 square kilometres of LiDAR data including Bloomington-Normal, Illinois, U.S.A. and surrounding areas with an average point density of 3.12 points/square metre classified into LAS Specification v1.2 codes: 1 (unclassified), 2 (ground), 7 (noise/low points), 9 (water), 10 (ignored ground: breakline proximity).

LAStools processing:
1)
check the quality of the LiDAR data [lasoverlap, lascontrol, lasinfo, lasgrid]
2)
merge and retile the original data [lastile]
3) classify point clouds into ground and non-ground [lasground]
4) create digital terrain (DTM) and digital surface models (DSM) [las2dem, blast2dem]
5) classify building and vegeration points [lasclassify]
6) extract building footprints [lasboundary]
7)
.produce height normalized tiles [lasheight]
8) generate a Canopy Height Model (CHM) with the workflow described here using the pit-free algorithm of Khosravipour et al. (2014) [lasthin, las2dem, lasgrid]

References:
Diuk-Wasser, M. A., Brown, H. E., Andreadis, T. G., Fish, D. 2006. Modeling the spatial distribution of mosquito vectors for West Nile virus in Connecticut, USA. Vector-Borne and Zoonotic Diseases 6: 283-295.
Dowling, Z., Ladeau, S. L., Armbruster, P., Biehler, D., Leisnham, P. T. 2013. Socioeconomic status affects mosquito (Diptera: Culicidae) larval habitat type availability and infestation level. Journal of Medical Entomology 50: 764-772.
Higa, Y., Yen, N. T., Kawada, H., Son, T. H., Hoa, N. T., Takagi, M. 2010. Geographic distribution of Aedes aegypti and Aedes albopictus collected from used tires in Vietnam. Journal of the American Mosquito Control Association 26: 1-9.
Khosravipour, A., Skidmore, A. K., Isenburg, M., Wang, T. J., Hussin, Y. A. 2014. Generating pit-free Canopy Height Models from Airborne LiDAR. Photogrammetric Engineering and Remote Sensing 80: 863-872.
Ruiz, M. O., Chaves, L. F., Hamer, G. L., Sun, T., Brown, W. M., Walker, E. D., Haramis, L., Goldberg, T. L. Kitron, U. D. 2010. Local impact of temperature and precipitation on West Nile virus infection in Culex species mosquitoes in northeast Illinois, USA. Parasites & vectors 3: 19.
Sun, X., Fu, S., Gong, Z., Ge, J., Meng, W., Feng, Y., Wang, J., Zhai, Y., Wang, H. H., Nasci, R. S., Tang, Q., Liang, G. 2009. Distribution of arboviruses and mosquitoes in Northwestern Yunnan Province, China. Vector-Borne and Zoonotic Diseases 9: 623-630.
USGS, HHS & CDC. 2015. Disease maps. http://diseasemaps.usgs.gov/mapviewer

LASmoons: Kiti Suomalainen

Kiti Suomalainen (recipient of three LASmoons)
Energy Centre
University of Auckland, New Zealand

Background:
Auckland enjoys 2050 hours of sunshine annually, comparable to Melbourne (2100) and Istanbul (2026), yet it is lagging behind in solar power installations. However, Auckland Council is committed to a sustainable pathway in mobility and energy consumption, aiming at solar photovoltaic (PV) installations powering an equivalent of over 176 500 homes by 2040 (approx. 37% of all homes), among other sustainability targets. The council has recently (2013/2014) conducted a LiDAR survey of the Auckland region (image 0), which they have provided to me for free for this project.

lasmoons_kiti_suomalainen_0

Goal:
This project aims to use this data, to quantitatively and accurately assess the solar potential of Auckland region rooftops, and provide the results free for the public. Image 1 gives a crude example of what the results may look like. Using LAStools I expect to efficiently get a more detailed DSM and corresponding elevation rasters. The goal is for any resident or user to be able to zoom in to any property within the extent of the collected LiDAR data and get an idea of the solar potential on that particular rooftop. Results will be given in average winter day, average summer day and average annual solar radiation per rooftop (or optimal x sq metres of rooftop – e.g. optimally placed 4 typical sized panels’ area).

lasmoons_kiti_suomalainen_1

Data:
+
appoximately 2250 aquare kilometres of LiDAR data collected in 2013.
+ average point density: 1.5 points per sq metre.

LAStools processing:
1)
create DTM tiles with 0.5 step, ground points only (classification 2), using the more efficient ‘.bil’ format [las2dem]
2) create DSM tiles with 0.5 step, first returns only, using the ‘.bil’ format [las2dem]
3) merge DTM and DSM tiles into single elevation raster [las2grid]
4).extract building footprints from the classified .las tiles (classification 6) for visualisation of final results, and reality checks [lasboundary]
5) use the raster files to calculate daily (winter day, summer day) and annual solar radiation on each rooftop (Solar Roof Tools by esri)

Reference:
Auckland Council, Low Carbon Auckland – Auckland’s Energy Resilience and Low Carbon Action Plan, July 2014.
NZ Aerial Mapping and Aerial Surveys Limited, LiDAR Flyover 2013/14 Project Final Report for Auckland Council, June 2015.

Removing low noise from Semi-Global Matching (SGM) Points

At PhoWo and INTERGEO 2015 rapidlasso was spending quality time with VisionMap who make the A3 Edge camera that provides fine resolution images from high altitudes and can quickly cover large areas. Under the hood of their LightSpeed software is the SURE dense-matching algorithm from nframes that turns images into photogrammetric point clouds. We were asked whether LAStools is able to create bare-earth DTM rasters from such points. If you have read our first, second, or third blog post on the topic you know that our asnwer was a resounding “YES!”. But we ran into an issue due to the large amount of low noise. Maybe the narrow angle between images at a high flying altitude affects the semi-global matching (SGM) algorithm. Either way, in the following we show how we use lascanopy and lasheight to mark low points as noise in a preprocessing step.

We obtained a USB stick containing a 2.42 GB file called “valparaiso_DSM_SURE_100.las” containing about 100 million points spaced 10 cm apart generated by SURE and stored with an (unnecessary high) resolution of millimeters (aka “resolution fluff”) as the third digit of all coordinates was always zero:

las2txt -i F:\valparaiso_DSM_SURE_100.las -stdout | more
255991.440 6339659.230 89.270
255991.540 6339659.240 89.270
255991.640 6339659.240 88.660
255991.740 6339659.230 88.730
[...]

We first compressed the bulky 2.42 GB LAS file into a compact 0.23 GB LAZ to our local hard drive – a file that is 10 times smaller and that will be 10 times faster to copy:

laszip -i F:\valparaiso_DSM_SURE_100.las ^
       -rescale 0.01 0.01 0.01 ^
       -o valparaiso_DSM_SURE_100.laz ^

Then we tiled the 100 million points into 250 meter by 250 meter tiles with 25 meter buffer using lastile. We use the new option ‘-flag_as_withheld’ to mark all buffer points with the withheld flag so they can be easily removed on-the-fly via the ‘-drop_withheld’ command-line filter (also see the README file).

lastile -i valparaiso_DSM_SURE_100.laz ^
        -tile_size 250 -buffer 25 ^
        -flag_as_withheld ^
        -odir valparaiso_tiles_raw -o val.laz
250 meter by 250 meter tiling with 25 meter buffer

250 meter by 250 meter tiles with 25 meter buffer

Before processing millions to billions of points we experiment with different options to find what works best on a smaller area, namely the tile “val_256750_6338500.laz” pointed to above. Using the workflow from this blog posts did not give perfect results due to the large amount of low noise. Although many low points were marked as noise (violett) by lasnoise, too many ended up classified as ground (brown) by lasground as seen here:
excessive low noise affects ground classification

excessive low noise affects ground classification

We use lascanopy – a tool very popular with forestry folks – to compute four BIL rasters where each 5m by 5m grid cell contains the 5th, 10th, 15th, and 20th percentile of the elevation values from all points falling into a cell (also see the README file):
lascanopy -i val_256750_6338500.laz ^
          -height_cutoff -1000 -step 5 ^
          -p 5 10 15 20 ^
          -obil
The four resulting rasters can be visually inspected and compared with lasview:
lasview -i val_256750_6338500_*.bil -files_are_flightlines
comparing 5th and 10th elevation percentiles

comparing the 5th and the 10th elevation percentiles

By pressing the hot keys <0>, <1>, <2> and <3> to switch between the different percentiles and <t> to triangulate them into a surface, we can see that for this example the 10th percentile works well while the 5th percentile is still affected by the low noise. Hence we use the 10th percentile elevation surface and classify all points below it as noise with lasheight (also see the README file).
lasheight -i val_256750_6338500.laz ^
          -ground_points val_256750_6338500_p10.bil ^
          -classify_below -0.5 7 ^
          -odix _denoised -olaz
We visually confirm that all low points where classified as noise (violett). Note the obvious “edge artifact” along the front boundary of the tile. This is why we always recommend to use a buffer in tile-based processing.
points below 10th percentile surface marked as noise

points below 10th percentile surface marked as noise

At the end of the blog post we give the entire command sequence that first computes a 10th percentile raster with 5m resolution for the entire file with lascanopy and then marks all points of each tile below the10th percentile surface as noise with lasheight. When we classify all points into ground and non-ground points with lasground we ignore all points classified as noise. Here are the results:

DTM extracted from SGM points despite low noise

DTM extracted from dense-matching points despite low noise

corresponding DSM with all buildings and vegetaion included

corresponding DSM with all buildings and vegetaion included

Above you see the generated DTM and the corresponding DSM. So yes, LAStools can create DTMs from points that are result of dense-matching photogrammetry … even when there is a lot of low noise. There are many other ways to mix and match the modules of LAStools for more refined workflows. Sometimes declaring all points below the 10th percentile surface as noise may be too agressive. In a future blog post we will look how to combine lascanopy and lasnoise for a more adaptive approach.

:: compute 10th percentile for entire area
lascanopy -i valparaiso_DSM_SURE_100.laz ^
          -height_cutoff -1000 -step 5 ^
          -p 10 ^
          -obil

:: tile input into 250 meter tiles with buffer
lastile -i valparaiso_DSM_SURE_100.laz ^
        -tile_size 250 -buffer 25 ^
        -flag_as_withheld ^
        -odir valparaiso_tiles_raw -o val.laz

:: mark points below as noise
lasheight -i valparaiso_tiles_raw/*.laz ^
          -ground_points valparaiso_DSM_SURE_100_p10.bil ^
          -classify_below -0.5 7 ^
          -odir valparaiso_tiles_denoised -olaz ^
          -cores 4

:: ground classify while ignoring noise points
 lasground -i valparaiso_tiles_denoised\*.laz ^
          -ignore_class 7 ^
          -town -bulge 0.5 ^
          -odir valparaiso_tiles_ground -olaz ^
          -cores 4 

:: create 50 cm DTM rasters in BIL format
las2dem -i valparaiso_tiles_ground\*.laz ^
        -keep_class 2 ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir valparaiso_tiles_dtm -obil ^
        -cores 4 

:: average 50 cm DTM values into single 1m DTM 
lasgrid -i valparaiso_tiles_dtm\*.bil -merged ^
        -step 1.0 -average ^
        -o valparaiso_dtm.bil

:: create hillshade adding in UTM 19 southern
blast2dem -i valparaiso_dtm.bil ^
          -hillshade -utm 19M ^
          -o valparaiso_dtm_hill.png

:: create DSM hillshade with same three steps
las2dem -i valparaiso_tiles_raw\*.laz ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir valparaiso_tiles_dsm -obil ^
        -cores 4
lasgrid -i valparaiso_tiles_dsm\*.bil -merged ^
        -step 1.0 -average ^
        -o valparaiso_dsm.bil
blast2dem -i valparaiso_dsm.bil ^
          -hillshade -utm 19M ^
          -o valparaiso_dsm_hill.png

Creating DTMs from dense-matched points of UAV imagery from SenseFly’s eBee

Tim Sutton and his team at Kartoza work on flood modelling and risk assessment using Inasafe. They have been trying to generate a DTM from point cloud data derived via dense-matching from UAV imagery collected by an eBee of SenseFly in the “unplanned developments” or “slums” North West of Dar es Salaam, the capital city of Tanzania. Tim’s team was stuck after “other software” produced this result:

results for ground points classification with other software

poor ground classification of “Tandale” with “other software”

Tim reached out to us at rapidlasso asking whether LAStools could handle this better. After all, we had published two blog articles – namely this one and that one – showing how to generate DTMs from the point clouds generated by the dense-matching photogrammetry software of Pix4D. Below the workflow we devised and the results we produced for Tim and his team.

We obtained 3 different data sets of areas called “Tandale”, “Borahatward”, and “Bugurunni”. We added one new option to our lasground software called ‘-bulge 1.0’ (see README) to improve the removal of smaller buildings and got this result for “Tandale”.

ground classification with LAStools

DTM of “Tandale” from ground points classified with LAStools

Before you point out the “facetted” look of this DTM keep in mind that “Tandale” is a densely populated poor area. A first hand account of the rough life in this area can be found here. Most dense-matching points are on corrugated roofs that become voids that need to be interpolated across in the DTM. Take a look at the corresponding DSM where all objects are still present.

original data

DSM of “Tandale” from all dense-matching points

Below we give a detailed description at the example of the “Bugurunni” data set of the workflow that was used to generate DTMs for the three data sets. At the end of this article you will see some more results.

We first use lassort to quantize, sort, and compress on 4 cores the seven spatially incoherent LAS files of the “Bugurunni” data set (totalling 4.5 GB with excessive resolution of millimeters) into LASzip-compressed files with a more reasonable resolution of centimeters and points ordered along a space-filling curve. We also add the missing projection information with ‘-utm 37M’. The resulting 7 LAZ files occupy only 0.7 GB meaning we get a compression of 9 : 1. The option ‘-odir’ specifies the output directory.

lassort -i bugurunni_densified_*.las ^
        -rescale 0.01 0.01 0.01 ^
        -utm 37M ^
        -odir bugurunni_strips -olaz ^
        -cores 4

Next we tile the sorted strips into 500 meter by 500 meter tiles with 50 meter buffer using lastile. We use the new option ‘-flag_as_withheld’ to mark all buffer points with the withheld flag so they can easily be removed on-the-fly with the ‘-drop_withheld’ command-line filter (see the README file for more on this).

lastile -i bugurunni_strips\*.laz ^
        -files_are_flightlines ^
        -tile_size 500 -buffer 50 ^
        -flag_as_withheld ^
        -o bugurunni_raw\bugu.laz
Using lasnoise on 4 cores we classify isolated points that might hinder ground-classification as noise (class 7). The parameters ‘-isolated 15’ means that all points surrounded by less than 15 other points in their 3 by 3 by 3 = 27 cells neighborhood in a 3D grid are considered isolated. The size of each grid cell is specified with ‘-step_xy 2 -step_z 1’  as 2 meter by 2 meter by 1 meter. These parameters were found experimentally (see the README file for more on this).
lasnoise -i bugurunni_raw\*.laz ^
         -step_xy 2 -step_z 1 ^
         -isolated 15 ^
         -odir bugurunni_noise -olaz ^
         -cores 4
Then we run lasground on 4 cores to classify the ground points with options ‘-metro’ and ‘-bulge 1.0’. The option ‘-metro’ is a convenient short-hand for ‘-step 50’ that will remove all objects on the terrain (e.g. large buildings) that have an extend of 50 meters or less. The option ‘-bulge 1.0’ instructs lasground to be conservative and only add points that are 1 meter or less above a smoothed version of the initial ground estimate (see the README file for more on this)..
lasground -i bugurunni_noise\*.laz ^
          -ignore_class 7 ^
          -metro -bulge 1.0 ^
          -odir bugurunni_ground -olaz ^
          -cores 4
Now we use las2dem to raster a DTM from only those points that were classified as ground. The option ‘-step 0.5’ sets the output grid resolution to 0.5 meters, ‘-kill 200’ interpolates across voids of up to 200 meters, and ‘-use_tile_bb’ rasters only the original 500 meter by 500 meter tile interior but not the 50 meter buffer. This assures that the resulting raster tiling aligns without artifacts across tile boundaries. The option ‘-obil’ chooses BIL as the output raster format.
las2dem -i bugurunni_ground\*.laz ^
        -keep_class 2 ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir bugurunni_dtm -obil ^
        -cores 4
As a simply form of anti-aliasing we average each four pixels of 0.5 meter resolution into one pixel of 1.0 meter resolution with lasgrid as all LAStools can read BIL files via on-the-fly conversion to points.
lasgrid -i bugurunni_dtm\*.bil -merged ^
        -step 1.0 -average ^
        -o bugurunni_dtm.bil

Finally we create a hillshade of the DTM adding back the projection that was “lost” in the BIL file generation so that blast2dem – the extremely scalable BLAST version of las2dem – can automatically produce a KML file for display in Google Earth.

blast2dem -i bugurunni_dtm.bil ^
          -hillshade -utm 37M ^
          -o bugurunni_dtm_hill.png

For comparison we also create a DSM with the same three steps but using all points.

las2dem -i bugurunni_raw\*.laz ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir bugurunni_dsm -obil ^
        -cores 4
lasgrid -i bugurunni_dsm\*.bil -merged ^
        -step 1.0 -average ^
        -o bugurunni_dsm.bil
blast2dem -i bugurunni_dsm.bil ^
          -hillshade -utm 37M ^
          -o bugurunni_dsm_hill.png
DTM of "Bugurunni" from ground points classified with LAStools

DTM of “Bugurunni” from ground points classified with LAStools

Above you see the generated DTM and below the corresponding DSM. So yes, LAStools can create DTMs from points that are result of dense-matching photogrammetry … under one assumption: there is not too much vegetation.

DSM of "Bugurunni" from all dense-matching points

DSM of “Bugurunni” from all dense-matching points

Below also the results for the “Borahatward” data. In a future blog post we will see how to deal with the excessive low noise sometimes present in dense-matching points.

DTM of "Bo"

DTM of “Borahatward” from ground points classified with LAStools

DSM of "Borahatward" from all dense-matching points

DSM of “Borahatward” from all dense-matching points

Use Buffers when Processing LiDAR in Tiles !!!

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:

Buffers are important to avoid edge artifacts along tile boundaries during DTM creation.

Buffers are important to avoid edge artifacts along tile boundaries.

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.

The

Christmas cookie analogy: buffers are like generously rolling out the dough

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:

  1. 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.
  2. 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.

tiles_buffered_with_lastile

Now we can forget about the buffers and run the standard workflow consiting of lasground, lasheight, and lasclassify to distinguish ground, vegetation, and building points in the LiDAR tiles.

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

At the end – when we generate raster products – we have to remember that the tiles were buffered by lastile and cut off the buffers when we raster the TIN with option ‘-use_tile_bb’ of las2dem.

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.

The DEM of the 9 tiles computed with buffers created by lastile has no edge artifacts acoss tile boundaries.

The DBM of the 9 tiles computed with buffers created by lastile has no edge artifacts acoss tile boundaries.

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 by 100 meters.

tiles_without_buffer

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.

The hillshaded DEM of the 9 tiles computed with on-the-fly buffering has no edge artifacts acoss tile boundaries.

The hillshaded DBM of 9 tiles computed with on-the-fly buffering has no edge artifacts acoss tile boundaries.

If you need to deliver the LiDAR files you should probably remove the buffers first … (-:

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.

The hillshaded DBM of 9 tiles computed WITHOUT using buffers has severe edge artifacts acoss tile boundaries.

The hillshaded DBM of 9 tiles computed WITHOUT using buffers has severe edge artifacts acoss tile boundaries.

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 … (-:

PS: Usually no buffers are needed for running lasgrid, lasoverlap, or lascanopy as they perform simple binning operations that do not make use of neighbour information.

LASmoons: Moreblessings Shoko

Moreblessings Shoko (recipient of three LASmoons)
School of Architecture, Planning and Geomatics
University of Cape Town, SOUTH AFRICA

Background:
Over one billion of the world’s population live in slums, lacking access to water, tenure, electricity, sanitation and basic services. Informal settlements are a growing challenge for urban governance especially in developing countries. These settlements not only thwart plans for coherent urban expansion, but also create a population that is vulnerable to risks in health, environmental, socially, politically and economically matters. This project seeks to provide improved spatial documentation for urban managemen to understand the spatial dynamics that lead to slum growth.

lasmoons_moreblessings_shokoGoal:
The aim of the research is to assess the efficiency of LiDAR in classification of an urban environment for the detection of informal settlements within South Africa. By creating more robust research tools that use three dimensional data such as LiDAR we create a platform for improving urban management especially regards the problem of slums which are a global concern. This work is part of an ongoing PhD research study in Cape Town, South Africa.

Data:
+ airborne LiDAR data covering entire extent of Cape Town, South Africa
+ Average point density: 4 points per square metre

LAStools processing:
1)
check data consistency [lasvalidate, lasinfo, lasprecision, lasduplicate]
2) organize raw LiDAR data into sufficiently small tiles with buffer [lastile]
3)
 distinguish ground from non-ground points [lasground]
4) remove outliers (points higher than 50 meters above ground) [lasheight]
5) classify points into building and vegetation returns [lasclassify]
6) triangulate and rasterize points into elevation, intensity and slope raster DTMs and DSMs [las2dem]

Reference:
image source: Google Earth

Discriminating Vegetation from Buildings

I came across an interesting blog article by Jarlath O’Neil-Dunne from the University of Vermont on how LiDAR return information can be used as a simple way to discriminate vegetated areas from buildings. He first computes a normalized first-return DSM and a normalized last-return DSM that he subtracts from another to highlight the vegetation. He writes “This is because the height difference of the first and last returns for buildings is often identical, whereas for trees it is typically much greater.”

Side note: I am not entirely happy with the terminology of a “Normalized Digital Terrain Model (nDTM)”. Jarlath writes: “A similar approach is used to create a Normalized Digital Terrain Model (nDTM).  A DTM is generated from the last returns. The DEM is then subtracted from the DTM to create the nDTM.” I like to reserve the term “Digital Terrain Model (DTM)” for bare-earth terrain computed from returns classified as ground.

Below I radically simplify Jarlath workflow by eliminating the two normalization steps. This not only saves the creation of 3 temporary rasters but also removes the requirement to have ground-classified LiDAR:

  1. Create a first-return frDSM
  2. Create a last-return lrDSM
  3. Subtract the lrDSM from the frDSM to get a return-difference rdDEM

This rdDEM has non-zero heights in all areas where the LiDAR produced more than one return. This happens most often and most pronounced in vegetated areas. Here is how to implement this with las2dem of LAStools:

las2dem -i ..\data\fusa.laz -first_only -o frDSM.bil
las2dem -i ..\data\fusa.laz -last_only -o lrDSM.bil
lasdiff -i frDSM.bil -i lrDSM.bil -o rdDEM.laz
lasview -i rdDEM.laz
rdDEM

The return-difference rdDEM shows the height difference between first and last returns.

Does this work well for you? The results on the “fusa.laz” data set are not entirely convincing … maybe because the vegetation was too dense (leaf-on?) so that the LiDAR penetration is not as pronounced. You can switch back and forth between the first-return and the last-return DSM by loading both *.bil files into lasview with the ‘-files_are_flightlines’ option and then press hotkeys ‘0’ and ‘1’ to toggle between the points and ‘t’ to triangulate the selected DSM.

lasview -i frDSM.bil lrDSM.bil -files_are_flightlines
first-return DSM

first-return DSM

last-return DSM

last-return DSM

We should point out that for Jarlath the return difference raster rdDEM is just one part of the pipeline that is followed by an object-based approach in which they integrate the spectral information from aerial imagery and then use iterative expert systems to further improve the tree canopy classification.

Nevertheless, we believe that our way of classifying vegetation and buildings via a pipeline of lasground, lasheight, and lasclassify gives a better and more robust initial guess than multi-return height differences towards what is vegetation and what are buildings. Below you see this is implemented using the new LASlayers concept:

lasground -i ..\data\fusa.laz -city -extra_fine -olay
lasheight  -i ..\data\fusa.laz -ilay -olay
lasclassify -i ..\data\fusa.laz -ilay -olay 
lasview -i ..\data\fusa.laz -ilay
Automated building and vegetation classification with lasclassify.

Automated building and vegetation classification with lasclassify.

Using lasgrid there are many ways that can easily turn the classified point cloud into a raster so that it can be used for subsequent exploitation together with other image data using a raster processing software. An example is shown below.

lasgrid -i ..\data\fusa.laz -ilay -keep_class 5 ^
        -step 0.5 -subcircle 0.1 -occupancy -fill 1 -false ^
        -use_bb -o vegetation.tif
lasgrid -i ..\data\fusa.laz -ilay -keep_class 6 ^
        -step 0.5 -subcircle 0.1 -occupancy -fill 1 -gray ^
        -use_bb -o buildings.tif
gdalwarp vegetation.tif buildings.tif classified.tif

classified

Alternatively we can use lasboundary to create a shapefile describing either the vegetation or the buildings.

lasboundary -i ..\data\fusa.laz -ilay -keep_class 5 ^
            -disjoint -concavity 1.5 -o vegetation.shp
lasboundary -i ..\data\fusa.laz -ilay -keep_class 6 ^
            -disjoint -concavity 1.5 -o buildings.shp
SHP file generated with lasboundary with polygons describing the vegetation.

SHP file generated with lasboundary with polygons describing the vegetation.

SHP file generated with lasboundary with polygons describing the buildings.

SHP file generated with lasboundary with polygons describing the buildings.