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 r.in.lidar and v.in.lidar modules. You can try them, e.g. usingr.in.lidar --help
Start GRASS GIS
Start GRASS - click on GRASS icon or typegrass83
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:
cd
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
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.out.file be_pointdensity
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
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.
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 input=elev_lid793016_be elevation=elev_lid793016_be_1m npmin=120 segmax=25 dmin=1
v.surf.rst 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
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.