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:

Software support

The GRASS GIS installation on your computer may not support and modules. You can try them, e.g. using --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 - click on GRASS icon or type

In startup pannel set GRASS GIS Database Directory to path to datasets, for example on MS Windows, C:\Users\myname\grassdata.
For GRASS Location select nc_spm_08_grass7 (North Carolina, State Plane, meters) and for GRASS Mapset create a new mapset (called e.g. HW_terrain_modeling).
Click Start GRASS session.

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:

Download all files (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 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 it.):
lasinfo tile_0793_016_spm.las
Class 2 represents bare earth points.
Now we import bare earth points and first return points separately: -t input=tile_0793_016_spm.las output=elev_lid793016_be class_filter=2 -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 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 representing number of points per 1 m grid cell.
Compare point densities for bare earth, first return.

g.region vect=elev_lid793016_1r res=1 -p input=tile_0793_016_spm.las output=lid_be_binn1m class_filter=2 method=n 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.rast lid_be_binn1m
d.legend lid_be_binn1m at=2,50,2,9 lid_be_binn1m unit=p,c
r.univar lid_be_binn1m
d.out.file be_pointdensity
d.rast lid_1r_binn1m
d.legend -d -s lid_1r_binn1m at=2,50,8,12 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. input=tile_0793_016_spm.las output=lid_be_binmean1m class_filter=2 method=mean 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 map=lid_1r_binmean1m
d.rast lid_be_binmean1m
d.legend lid_be_binmean1m at=2,40,2,5
d.out.file lidbemean1m
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 input=tile_0793_016_spm.las output=lid_be_binmean2m class_filter=2 method=mean 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.rast lid_be_binmean2m
d.legend lid_be_binmean2m at=2,40,2,5
d.out.file lidbemean2m
d.rast lid_1r_binmean2m
d.legend lid_1r_binmean2m at=2,40,2,5
d.out.file lid1rmean2m

Compute range of elevation values in each grid cell: input=tile_0793_016_spm.las output=lid_be_binrange2m class_filter=2 method=range
r.colors lid_be_binrange2m color=bcyr -e input=tile_0793_016_spm.las output=lid_1r_binrange2m return_filter=first method=range
r.colors lid_1r_binrange2m color=bcyr -e
d.rast lid_be_binrange2m
d.legend lid_be_binrange2m at=2,40,2,5
d.out.file lid_be_binrange2m
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. 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.
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 and what is the range within cell.
If we need a 1m resolution DEM or DSM for our application this analysis tells us that we need to interpolate it from the point cloud. What steps would you begin with when processing point cloud data you are not familiar with?

Interpolate DEM and DSM

To interpolate DEM and DSM we use default parameters except for number of points used for segmentation and interpolation
(segmax and npmin and higher tension for multiple return data).
You can set dmin=1 to make the interpolation run faster (see the manual).
Be patient, it can take a few minutes to run depending on the computer power.
g.region res=1 -ap input=elev_lid793016_be elevation=elev_lid793016_be_1m npmin=120 segmax=25 dmin=1 input=elev_lid793016_1r elevation=elev_lid793016_1r_1m npmin=120 segmax=25 tension=100 smooth=0.5 dmin=1
r.colors elev_lid793016_be_1m color=elevation
r.colors elev_lid793016_1r_1m color=elevation
d.rast elev_lid793016_be_1m
d.rast elev_lid793016_1r_1m
Hide 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

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, an online point cloud visualization tool to visualize the LAS file from this assignment. Explore the provided visualization tools in and create images showing different properties of the point cloud or the visualization techniques and how they would be useful for further analysis.