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
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:
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.
cd
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?