GIS/MEA582: Geospatial Modeling and Analysis

Geomorphometry I: Terrain Modeling

Resources: GRASS GIS overview and manual, libLAS tools for lidar data conversions

Download the Raleigh 2013 lidar data as LAS file and orthophoto (download individual files to avoid issues with unzipping):

Software support

The GRASS GIS installation on your computer may not support r.in.lidar and v.in.lidar modules. You can try them, e.g. using

r.in.lidar --help

On Windows platform, please use the standalone version of GRASS GIS we recommend in Course Logistics. If the modules don't work for you, please use instead of this assignment the ASCII file based version of the assignment.

Start GRASS GIS

Create a new mapset (called e.g. HW_terrain_modeling) in nc_spm_08_grass7 project and change working directory:
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI Console and hit Enter:

cd

Download all files one-by-one (see above) to the selected directory. Now you can use the commands from the assignment requiring the file without the need to specify the full path to the file.

Analyze bare earth and multiple return lidar data properties by binning

Import the points using v.in.lidar. We can specify which class and which return (first, middle, last) we want to import.
We can see the classification either in metadata distributed with the lidar data or it can be displayed with lasinfo tool (in case lasinfo command is not available, skip this step.):

lasinfo tile_0793_016_spm.las

Class 2 represents bare earth points.
Now we import bare earth points and first return points separately (add -o flag to the command if you are getting error that the projection does not match the project projection):

v.in.lidar -t input=tile_0793_016_spm.las output=elev_lid793016_be class_filter=2
v.in.lidar -t input=tile_0793_016_spm.las output=elev_lid793016_1r return_filter=first

Note: if you have any problems with the files, review the instructions above or use v.in.lidar GUI dialog to browse to get the path to the files.

Set the region from the imported point file with resolution of 1 meter. Compute raster maps (with r.in.lidar) representing number of points per 1 m grid cell (add -o flag to the r.in.lidar command if you are getting error that the projection does not match the project projection): Compare point densities for bare earth, first return. How do the point densities for bare earth and first return compare?

g.region vect=elev_lid793016_1r res=1 -p
r.in.lidar input=tile_0793_016_spm.las output=lid_be_binn1m class_filter=2 method=n
r.in.lidar input=tile_0793_016_spm.las output=lid_1r_binn1m return_filter=first method=n
r.colors lid_be_binn1m,lid_1r_binn1m color=bcyr -e
d.erase
d.rast lid_be_binn1m
d.legend lid_be_binn1m at=2,50,2,9
r.report lid_be_binn1m unit=p,c
r.univar lid_be_binn1m
d.out.file be_pointdensity
d.erase
d.rast lid_1r_binn1m
d.legend -d -s lid_1r_binn1m at=2,50,8,12
r.report lid_1r_binn1m unit=p
r.univar lid_1r_binn1m
d.out.file first_pointdensity

Compute a raster map representing mean elevation for each 1m cell.
It will have holes.

r.in.lidar input=tile_0793_016_spm.las output=lid_be_binmean1m class_filter=2 method=mean
r.in.lidar input=tile_0793_016_spm.las output=lid_1r_binmean1m return_filter=first method=mean
r.colors lid_be_binmean1m color=elevation
r.colors lid_1r_binmean1m color=elevation
r.info map=lid_1r_binmean1m
d.erase
d.rast lid_be_binmean1m
d.legend lid_be_binmean1m at=2,40,2,5
d.out.file lidbemean1m
d.erase
d.rast lid_1r_binmean1m
d.legend lid_1r_binmean1m at=2,40,2,5
d.out.file lid1rmean1m

Compute a raster map representing mean elevation for each 2m cell.
Result is almost good enough for 1st return, but there are still many holes for bare earth.

g.region res=2 -ap
r.in.lidar input=tile_0793_016_spm.las output=lid_be_binmean2m class_filter=2 method=mean
r.in.lidar input=tile_0793_016_spm.las output=lid_1r_binmean2m return_filter=first method=mean
r.colors lid_be_binmean2m color=elevation
r.colors lid_1r_binmean2m color=elevation
d.erase
d.rast lid_be_binmean2m
d.legend lid_be_binmean2m at=2,40,2,5
d.out.file lidbemean2m
d.erase
d.rast lid_1r_binmean2m
d.legend lid_1r_binmean2m at=2,40,2,5
d.out.file lid1rmean2m

Why are there holes in the binned DEMs? Why is the first return map mostly green (where is the orange-brown elevation?)

Compute range of elevation values in each grid cell:

r.in.lidar input=tile_0793_016_spm.las output=lid_be_binrange2m class_filter=2 method=range
r.colors lid_be_binrange2m color=bcyr -e
r.in.lidar input=tile_0793_016_spm.las output=lid_1r_binrange2m return_filter=first method=range
r.colors lid_1r_binrange2m color=bcyr -e
d.erase
d.rast lid_be_binrange2m
d.legend lid_be_binrange2m at=2,40,2,5
d.out.file lid_be_binrange2m
d.erase
d.rast lid_1r_binrange2m
d.legend lid_1r_binrange2m at=2,40,2,5
d.out.file lid_1r_binrange2m

Import and display orthophoto, switch off all layers except for orthophoto.

r.in.gdal ortho_0793_016_ncspm.tif output=ortho_2013_0793
d.rast ortho_2013_0793

Identify the features that are associated with large range values.
Display only the high values of range (objects taller than 10m).
What landcover is associated with large range in multiple return data?

d.rast lid_1r_binrange2m values=10.-200.
d.out.file mylid_1rrange2m

We now know how dense the data are at 1m and 2m resolution and what is the range of elevation values within a grid cell.
If we need a 1m resolution DTM or DSM for our application this analysis tells us that the fast binning is not enough and we need to interpolate it from the point cloud.

Interpolate DTM and DSM

To interpolate DTM and DSM using the RST spline function we use default parameters except for number of points used for interpolation (npmin), minimum distance between points and higher tension and smoothing for multiple return data.
Be patient, it can take a few minutes to run depending on your computer power.
We derive shaded relief for the interpolated DTM and DSM and display the result as colored shaded topography to highlight the terrain features.

g.region res=1 -ap
v.surf.rst input=elev_lid793016_be elevation=elev_lid793016_be_1m npmin=120 dmin=1
v.surf.rst input=elev_lid793016_1r elevation=elev_lid793016_1r_1m npmin=120 tension=100 smooth=0.5 dmin=1
r.relief input=elev_lid793016_be_1m output=elev_lid793016_be_1msh
r.relief input=elev_lid793016_1r_1m output=elev_lid793016_1r_1msh
d.erase
d.shade shade=elev_lid793016_be_1msh color=elev_lid793016_be_1m
d.legend elev_lid793016_be_1m
d.out.file elevation_be_1m
d.erase
d.shade shade=elev_lid793016_1r_1msh color=elev_lid793016_1r_1m
d.legend elev_lid793016_1r_1m
d.out.file elevation_1r_1m

How does the binned first return DSM compare to interpolated first return DSM and what is the reason for the differences in range and coverage?

Compute a map of above ground height of vegetation and buildings from your interpolated bare ground and first return DEMs using map algebra. Include the output map with appropriate color ramp in your report. What is the mean height computed from this map?

Optional: Visualize the interpolated DEMs in 3D with cutting planes

Remove legend and switch off all map layers except for the last 2 interpolated ones.
Use 3D view with cutting planes to compare the bare earth and terrain surface.
Make sure fine resolution is set to 1 for all surfaces.
Assign each surface constant color, add constant plane at 75m elevation for reference.
Shade the crossection using the color by bottom surface option.
If you don't remember this, see screen capture video for 3d view.

Save image for your report.

Optional: Visualize the point cloud using plas.io

Visualizing the raw points in point clouds can be time consuming both in 2D and 3D, however several tools for doing so in an efficient way exist. Try to use plas.io, an online point cloud visualization tool to visualize the LAS file from this assignment. Explore the provided visualization tools in plas.io and create images showing different properties of the point cloud or the visualization techniques and how they would be useful for further analysis.