<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Geomorphometry I
Step1:  In startup pannel set GIS Data Directory to path to datasets,
for example on MS Windows, C
Step2:  Range of coordinates at lower resolution
Patch and overlay planimetry to see that there are
no points in areas with buildings
Step3:  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.
Step4:  Compute a raster map representing mean elevation for each 6m cell.
It will still have a few holes.
Step5:  Compute range
Step6:  Identify the features that are associated with large range values.
Display with vector data.
Step7:  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?
Step8:  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.
Step9:  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.
Step10:  Display bare ground and multiple return points over orthophoto
Step11:  Display our imported points
Step12:  Display points available in the data set
(elev_lid_bepts and elev_lid792_bepts should be identical)
Step13:  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
Step14:  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
Step15:  Condition for subset in GUI
Step16:  Can you guess what is in the area that does not have any 1st return points?
<ASSISTANT_TASK:>
Python Code:
# 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)
# 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")
# 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")
# 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)
# 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'
Explanation: 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
End of explanation
!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")
Explanation: 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:
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.
End of explanation
!v.patch P079214,P079215,P079218,P079219 out=planimetry_rural
!d.vect planimetry_rural
Image(filename="map.png")
Explanation: Range of coordinates at lower resolution
Patch and overlay planimetry to see that there are
no points in areas with buildings:
End of explanation
!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")
Explanation: 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.
End of explanation
!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")
Explanation: Compute a raster map representing mean elevation for each 6m cell.
It will still have a few holes.
End of explanation
!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")
Explanation: Compute range:
End of explanation
!d.vect planimetry_rural
!d.vect lakes type=boundary co=violet
!d.vect streams co=blue
Image(filename="map.png")
Explanation: Identify the features that are associated with large range values.
Display with vector data.
End of explanation
!g.region rural_1m -p
Explanation: 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?
End of explanation
!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")
Explanation: 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.
End of explanation
!g.region rural_1m -p
!v.in.ascii -ztr input=lid_be_pts.txt out=elev_lid_bepts z=3
Explanation: 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.
End of explanation
!d.erase
!d.rast ortho_2001_t792_1m
Image(filename="map.png")
Explanation: Display bare ground and multiple return points over orthophoto:
End of explanation
!d.vect elev_lid_bepts size=2 col=red
Image(filename="map.png")
Explanation: Display our imported points:
End of explanation
!d.vect elev_lidrural_mrpts size=4 col=0:100:0
!d.vect elev_lid792_bepts size=2 col=yellow
Image(filename="map.png")
Explanation: Display points available in the data set
(elev_lid_bepts and elev_lid792_bepts should be identical):
End of explanation
!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")
Explanation: 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:
End of explanation
!g.region rural_1m -p
!d.erase
!d.rast ortho_2001_t792_1m
Image(filename="map.png")
Explanation: 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:
End of explanation
!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")
Explanation: 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.
End of explanation
# end the GRASS session
os.remove(rcfile)
Explanation: Can you guess what is in the area that does not have any 1st return points?
End of explanation
<END_TASK>