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 - click on GRASS icon or type
In [ ]:
# Obtain sample data and set new Grass mapset
import urllib
from zipfile import ZipFile
import os.path
zip_path = "/home/jovyan/work/tmp/nc_spm_08_grass7.zip"
mapset_path = "/home/jovyan/grassdata"
if not os.path.exists(zip_path):
urllib.urlretrieve("https://grass.osgeo.org/sampledata/north_carolina/nc_spm_08_grass7.zip", zip_path)
if not os.path.exists(mapset_path):
with ZipFile(zip_path, 'r') as zip_ref:
zip_ref.extractall(mapset_path)
In [ ]:
# a proper directory is already set, download files
import urllib
urllib.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/lid_be_pts.txt", "lid_be_pts.txt")
urllib.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/lid_mr_pts.txt", "lid_mr_pts.txt")
urllib.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/lid_be_pts.txt", "lid_be_pts.txt")
In [ ]:
# using Python to initialize GRASS GIS
import os
import sys
import subprocess
from IPython.display import Image
# create GRASS GIS runtime environment
gisbase = subprocess.check_output(["grass", "--config", "path"]).strip()
os.environ['GISBASE'] = gisbase
sys.path.append(os.path.join(gisbase, "etc", "python"))
# do GRASS GIS imports
import grass.script as gs
import grass.script.setup as gsetup
# set GRASS GIS session data
rcfile = gsetup.init(gisbase, "/home/jovyan/grassdata", "nc_spm_08_grass7", "user1")
In [ ]:
# using Python to initialize GRASS GIS
# default font displays
os.environ['GRASS_FONT'] = 'sans'
# overwrite existing maps
os.environ['GRASS_OVERWRITE'] = '1'
gs.set_raise_on_error(True)
gs.set_capture_stderr(True)
In [ ]:
# using Python to initialize GRASS GIS
# set display modules to render into a file (named map.png by default)
os.environ['GRASS_RENDER_IMMEDIATE'] = 'cairo'
os.environ['GRASS_RENDER_FILE_READ'] = 'TRUE'
os.environ['GRASS_LEGEND_FILE'] = 'legend.txt'
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
Console and hit Enter:
In [ ]:
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. 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.
In [ ]:
!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
Image(filename="map.png")
!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
Image(filename="map.png")
In [ ]:
!v.patch P079214,P079215,P079218,P079219 out=planimetry_rural
!d.vect planimetry_rural
Image(filename="map.png")
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.
In [ ]:
!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
Image(filename="map.png")
!r.report lid_be_binn6m unit=p
!r.univar lid_be_binn6m
Image(filename="map.png")
Compute a raster map representing mean elevation for each 6m cell. It will still have a few holes.
In [ ]:
!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
Image(filename="map.png")
Compute range:
In [ ]:
!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
Image(filename="map.png")
Identify the features that are associated with large range values. Display with vector data.
In [ ]:
!d.vect planimetry_rural
!d.vect lakes type=boundary co=violet
!d.vect streams co=blue
Image(filename="map.png")
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?
In [ ]:
!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.
In [ ]:
!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.
Image(filename="map.png")
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.
In [ ]:
!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:
In [ ]:
!d.erase
!d.rast ortho_2001_t792_1m
Image(filename="map.png")
Display our imported points:
In [ ]:
!d.vect elev_lid_bepts size=2 col=red
Image(filename="map.png")
Display points available in the data set (elev_lid_bepts and elev_lid792_bepts should be identical):
In [ ]:
!d.vect elev_lidrural_mrpts size=4 col=0:100:0
!d.vect elev_lid792_bepts size=2 col=yellow
Image(filename="map.png")
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:
In [ ]:
!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
Image(filename="map.png")
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.
Find out where we have multiple returns:
In [ ]:
!g.region rural_1m -p
!d.erase
!d.rast ortho_2001_t792_1m
Image(filename="map.png")
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.
In [ ]:
!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
Image(filename="map.png")
Can you guess what is in the area that does not have any 1st return points?
In [ ]:
# end the GRASS session
os.remove(rcfile)