Resources:
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 [ ]:
# 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_data_models).
Click Start GRASS.
If you prefer to work in GUI, you should be able to find out yourself the GUI equivalents for the tasks below. Some hints for GUI are included, but from now on, most of the instructions will be provided as commands for command line. Hint for running most of the commands in GUI - type or paste the name of the module into the command console in the Console tab and then hit Enter to open the GUI dialog. Read the manual page for each command you are using for the first time to learn what it is doing and what the parameters mean.
Resample the given raster map to higher and lower resolution (30m->10m, 30m->100m) and compare resampling by nearest neighbor with bilinear and bicubic method. First, set the computation region extent to our study area and set resolution to 30 meters. The computational region (region for short) is set using g.region module. Here for convenience we use named region which defines both the extent and the resolution. This named region is included in the data (location) we are using but it is possible to create new named regions and use them to bookmark different study areas.
In [ ]:
!g.region region=swwake_30m -p
The -p flag for g.region is used to print the region
we just set.
Then we display the 30m resolution NED elevation raster.
In [ ]:
!d.rast elev_ned_30m
Image(filename="map.png")
To resample it to 10m resolution, first set the computational region to resolution 10m, then resample the raster using the nearest neighbor method. Hint: To open the r.resamp.interp in GUI, type or paste the module name into the Console tab, then Enter to open the GUI dialog, don't forget to set the method to nearest under Optional tab.
In [ ]:
!g.region res=10 -p
!r.resamp.interp elev_ned_30m out=elev_ned10m_nn method=nearest
Display the resampled map by adding "elev_ned10m_nn" to Layer Manager in case you don't have it in the Layer Manager already. Alternatively, use in command line the following:
In [ ]:
!d.rast elev_ned10m_nn
Image(filename="map.png")
The elevation map "elev_ned10m_nn" looks the same as the original one, so now check the resampled elevation surface using the aspect map:
In [ ]:
!r.slope.aspect elevation=elev_ned10m_nn aspect=aspect_ned10m_nn
Display the resampled map by adding "aspect_ned10m_nn" to Layer Manager or in command line using:
In [ ]:
!d.rast aspect_ned10m_nn
Image(filename="map.png")
To save the map, click in Map Display to on the button Save display to graphic file" or alternatively, use the following command:
In [ ]:
Image(filename="map.png")
Now, reinterpolate DEMs using bilinear and bicubic interpolation. Check the interpolated elevation surfaces using aspect maps.
In [ ]:
!r.resamp.interp elev_ned_30m out=elev_ned10m_bil meth=bilinear
!r.resamp.interp elev_ned_30m out=elev_ned10m_bic meth=bicubic
!r.slope.aspect elevation=elev_ned10m_bil aspect=aspect_ned10m_bil
!r.slope.aspect elevation=elev_ned10m_bic aspect=aspect_ned10m_bic
!d.rast aspect_ned10m_bil
!d.rast aspect_ned10m_bic
Image(filename="map.png")
Save the displayed maps and in your report, compare the results with the previously computed nearest neighbor result. In Map Display click button Save display to graphic file, or use the following in the command line:
In [ ]:
Image(filename="map.png")
In [ ]:
!d.rast elev_ned_30m
!d.rast landuse96_28m
Image(filename="map.png")
Then change the region resolution and resample elevation (which is a continuous field) and land use (which has discrete categories). Explain selection of aggregation method. Can we use average also for landuse? What does mode mean?
In [ ]:
!g.region res=100 -p
!r.resamp.stats elev_ned_30m out=elev_new100m_avg method=average
!d.rast elev_new100m_avg
Image(filename="map.png")
Before the next computation, remove all map layers from the Layer Manager because we don't need to see them anymore.
In [ ]:
!d.erase
!r.resamp.stats landuse96_28m out=landuse96_100m method=mode
!d.rast landuse96_100m
Image(filename="map.png")
In [ ]:
!v.to.points census_wake2000 type=centroid out=census_centr use=vertex
Display census boundaries using GUI: Add vector "census_wake2000" Selection > Feature type > boundary (switch off the other types). Save the displayed map in Map Display click button Save display to graphic file. Alternatively, use the following commands to control display.
Note that in both command line and GUI you must either enter the full path to the file you are saving the image in, or you must know the current working directory.
In [ ]:
!d.vect census_centr icon=basic/circle fill_color=green size=10
!d.vect census_wake2000 color=red fill_color=none
!d.legend.vect
Image(filename="map.png")
Convert contour lines to points (useful for computing DEM from contours):
In [ ]:
!v.to.points input=elev_ned10m_cont10m output=elev_ned_contpts type=line use=vertex
Display the "elev_ned_contpts" points vector and zoom-in to very small area to see the actual points.
In [ ]:
!d.vect elev_ned_contpts co=brown icon=basic/point size=3
Image(filename="map.png")
In [ ]:
!g.region swwake_30m res=200 -p
Then remove all layers from the Layer Manager.
Convert vector points "schools" to raster. As value for raster use attribute column "CORECAPACI" for core capacity. To add legend in GUI use Add map elements > Show/hide legend and select "schools_cap_200m".
In [ ]:
!d.vect schools_wake
!v.info -c schools_wake
!v.to.rast schools_wake out=schools_cap_200m use=attr attrcol=CORECAPACI type=point
!d.rast schools_cap_200m
!d.vect streets_wake co=grey
!d.legend schools_cap_200m at=70,30,2,6
Image(filename="map.png")
Now convert lines in "streets" vector to raster. Set the resolution to 30m and use speed limit attribute.
In [ ]:
!g.region res=30 -p
!v.to.rast streets_wake out=streets_speed_30m use=attr attrcol=SPEED type=line
If you haven't done this already, add remove all other map layers from Layer Manager and add "streets_speed_30m" raster layer. Add legend for "streets_speed_30m" raster using GUI in Map Display: Add legend > Set Options > Advanced > List of discrete cat numbers and type in speed limits 25,35,45,55,65; move legend with mouse as needed.
Alternatively, use the following commands:
In [ ]:
!d.erase
!d.rast streets_speed_30m
!d.legend streets_speed_30m at=5,30,2,5 use=25,35,45,55,65
Image(filename="map.png")
Save the displayed map. In Map Display click button Save display to graphic file, or use the following.
In [ ]:
Image(filename="map.png")
Convert raster lines to vector lines.
First, set the region and remove map layers from Layer Manager. Then do the conversion.
Explain why we are using r.thin module. You may want to remove all previously used layers from the Layer Manager before you start these new computations.
In [ ]:
!d.erase
!g.region raster=streams_derived -p
!d.rast streams_derived
!r.thin streams_derived output=streams_derived_t
!r.to.vect streams_derived_t output=streams_derived_t type=line
Image(filename="map.png")
Visually compare the result with streams digitized from airphotos.
In [ ]:
!d.vect streams_derived_t color=blue
!d.vect streams color=red
Image(filename="map.png")
Save the displayed map (in Map Display click button Save display to graphic file).
In [ ]:
Image(filename="map.png")
Convert raster areas representing basins to vector polygons.
Use raster value as category number (flag -v) and display vector polygons filled with random colors. In GUI: Add vector > Colors > Switch on Random colors. You may want to remove all previously used layers from the Layer Manager before you start these new computations.
In [ ]:
!g.region raster=basin_50K -p
!d.erase
!d.rast basin_50K
!r.to.vect -sv basin_50K output=basin_50Kval type=area
!d.vect -c basin_50Kval
!d.vect streams color=blue
Image(filename="map.png")
Save the displayed map either using GUI or using the following in case you are working in the command line.
In [ ]:
Image(filename="map.png")
In [ ]:
# end the GRASS session
os.remove(rcfile)