NCSU GIS/MEA582:
Geospatial Modeling and Analysis

Geomorphometry I: Terrain Modeling (data in ASCII files)

Resources:

Download the ASCI (x,y,z) lidar bare earth data lid_be_pts.txt
Download the ASCI (x,y,z) lidar multiple return data lid_mr_pts.txt

Start GRASS GIS

Start GRASS - click on GRASS icon or type
grass83

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

Bare earth data

Analyze bare earth and multiple return lidar data properties. First, download the ASCI (x,y,z) lidar data lid_be_pts.txt and then compute a raster map representing number of points per 2m and 6m grid cell to assess the point density. If you are unsure about the current working directory and where your text file with point is, run r.in.xyz from GUI to find the path to the external lid_be_pts.txt. Do not forget to zoom to computational region to check its extent and see the resulting data. You can use horizontal legends by resizing legend into wide short rectangle or by using at=6,10,2,30 in the command line. You can resize the Map Display to create white space above and below the raster map image.
g.region rural_1m res=2 -p
r.in.xyz input=lid_be_pts.txt output=lid_be_binn2m method=n
r.in.xyz input=lid_mr_pts.txt output=lid_mr_binn2m method=n
d.erase
d.rast lid_mr_binn2m
d.legend lid_mr_binn2m at=2,20,2,5
d.out.file mylidmrn2m
r.report lid_mr_binn2m unit=p
r.univar lid_mr_binn2m
d.rast lid_be_binn2m
d.legend lid_be_binn2m at=2,20,2,5
r.report lid_be_binn2m unit=p
r.univar lid_be_binn2m

Range of coordinates at lower resolution

Patch and overlay planimetry to see that there are no points in areas with buildings:

v.patch P079214,P079215,P079218,P079219 out=planimetry_rural
d.vect planimetry_rural
d.out.file mylidben2m

Decrease resolution and try the above steps with the r.in.xyz module again. Again, run it from GUI, define full path, or manage your current working directory.

g.region rural_1m res=6 -ap
r.in.xyz input=lid_be_pts.txt out=lid_be_binn6m meth=n
d.erase
d.rast lid_be_binn6m
d.legend lid_be_binn6m at=2,20,2,5
d.out.file mylidben6m
r.report lid_be_binn6m unit=p
r.univar lid_be_binn6m

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

r.in.xyz input=lid_be_pts.txt out=lid_be_binmean6m meth=mean
r.colors lid_be_binmean6m color=elevation
d.rast lid_be_binmean6m
d.legend lid_be_binmean6m at=2,40,2,5
r.in.xyz input=lid_mr_pts.txt out=lid_mr_binmean6m meth=mean
r.colors lid_mr_binmean6m co=elevation
d.rast lid_mr_binmean6m
d.legend lid_mr_binmean6m at=2,40,2,5
d.out.file mylidmrmean6m

Compute range:

r.in.xyz input=lid_be_pts.txt out=lid_be_binrange6m meth=range
r.in.xyz input=lid_mr_pts.txt out=lid_mr_binrange6m meth=range
d.erase
d.rast lid_be_binrange6m
d.legend lid_be_binrange6m at=2,40,2,5
d.rast lid_mr_binrange6m
d.legend lid_mr_binrange6m at=2,40,2,5

Identify the features that are associated with large range values. Display with vector data.

d.vect planimetry_rural
d.vect lakes type=boundary co=violet
d.vect streams co=blue
d.out.file mylidmrrange6m

Display only the high values of range (0.5-20m) overlayed over orthophoto. What landcover is associated with large range in multiple return data?

Which landscape features may be lost at 6m resolution?

g.region rural_1m -p

Do not forget to zoom/set the display to computational region to display only selected interval of values in GUI. Add raster > Required tab, select lid_be_binrange6m In Selection tab, set List of values to be displayed to 0.5-20.

d.erase
d.rast ortho_2001_t792_1m
d.rast lid_be_binrange6m val=0.5-20.
d.erase
d.rast ortho_2001_t792_1m
d.rast lid_mr_binrange6m val=0.5-20.

Interpolation

We now know how dense the data are and what is the range within cell.
If we need 1m resolution DEM for our application this analysis tells us that we need to interpolate.

Import the ascii lidar data as vector points.
Import only points in the rural area without building a table (-t flag in tab Points), number of column used as z is 3 (third column), and using z-coordinate for elevation (Tab Optional, Create 3D vector map).
We assume that the txt file is in your current working directory.

g.region rural_1m -p
v.in.ascii -ztr input=lid_be_pts.txt out=elev_lid_bepts z=3

Display bare ground and multiple return points over orthophoto:

d.erase
d.rast ortho_2001_t792_1m

Display our imported points:

d.vect elev_lid_bepts size=2 col=red

Display points available in the data set (elev_lid_bepts and elev_lid792_bepts should be identical):

d.vect elev_lidrural_mrpts size=4 col=0:100:0
d.vect elev_lid792_bepts size=2 col=yellow

Extract first return to get points for DSM.
Interpolate DEM and DSM (more about interpolation in the interpolation assignment).
We use default parameters except for number of points used for segmentation and interpolation - segmax and npmin and higher tension for multiple return data:

g.region rural_1m -p
v.extract elev_lidrural_mrpts out=elev_lidrur_first type=point where="Return=1"
v.surf.rst input=elev_lid792_bepts elevation=elev_lidrural_1m npmin=120 segmax=25
v.surf.rst input=elev_lidrur_first elevation=elev_lidrurfirst_1m npmin=120 segmax=25 tension=100 smooth=0.5
d.erase
d.rast elev_lidrural_1m
d.rast elev_lidrurfirst_1m

Remove all layers except for elev_lidrural_1m and elev_lidrurfirst_1m and switch to 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 z=90 for reference. Create crossections using cutting plane, shade the crossection using the color by bottom surface option.
Save image for report. If you don't remember this, see screen capture video for GRASS GIS 3D view.

Multiple returns

Find out where we have multiple returns:

g.region rural_1m -p
d.erase
d.rast ortho_2001_t792_1m

Condition for subset in GUI: Add vector > Selection > type return=1 into WHERE condition of SQL statement.
You need to add the points 4 times to create a map that will have
all sets of returns.

d.vect elev_lidrural_mrpts where=return=1 col=red size=2
d.vect elev_lidrural_mrpts where=return=2 col=green size=3
d.vect elev_lidrural_mrpts where=return=3 col=blue
d.vect elev_lidrural_mrpts where=return=4 col=yellow
d.out.file lidreturns

Can you guess what is in the area that does not have any 1st return points?