GIS/MEA582: Geospatial Modeling and Analysis

This is an unmaintained course material, please see current material at:

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:

Start GRASS GIS

Start GRASS - click on GRASS icon or type

grass70

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

Change working directory:
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI command console and hit Enter:

cd

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 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 it.):

lasinfo tile_0793_016_spm.las

Class 2 represents bare earth points.
Now we import bare earth points and first return points separately:

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.
Compare point densities for bare earth, first return.

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

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.rast lid_be_binmean1m
d.legend lid_be_binmean1m at=2,40,2,5
d.rast lid_1r_binmean1m
d.legend lid_1r_binmean1m at=2,40,2,5
d.out.file mylidmrmean1m

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.rast lid_be_binmean2m
d.legend lid_be_binmean2m at=2,40,2,5
d.rast lid_1r_binmean2m
d.legend lid_1r_binmean2m at=2,40,2,5
d.out.file mylidmrmean2m

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.in.lidar input=tile_0793_016_spm.las output=lid_1r_binrange2m return_filter=first method=range
d.rast lid_be_binrange2m
d.legend lid_be_binrange2m at=2,40,2,5
d.rast lid_1r_binrange2m
d.legend lid_1r_binrange2m at=2,40,2,5

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.
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
v.surf.rst elev_lid793016_be elevation=elev_lid793016_be_1m npmin=120 segmax=25 dmin=1
v.surf.rst 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 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.