<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Geospatial data models
Resources
Step1:  In startup pannel set GIS Data Directory to path to datasets,
for example on MS Windows, C
Step2:  The -p flag for g.region is used to print the region
we just set.
Then we display the 30m resolution NED elevation raster.
Step3:  To resample it to 10m resolution, first set the computational region to resolution 10m,
then resample the raster using the nearest neighbor method.
Hint
Step4:  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
Step5:  The elevation map "elev_ned10m_nn" looks the same as the original one,
so now check the resampled elevation surface using the aspect map
Step6:  Display the resampled map by adding "aspect_ned10m_nn" to Layer Manager
or in command line using
Step7:  To save the map, click in Map Display to on the button
Save display to graphic file" or alternatively,
use the following command
Step8:  Now, reinterpolate DEMs using bilinear and bicubic interpolation.
Check the interpolated elevation surfaces using aspect maps.
Step9:  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
Step10:  Why is the aspect of elevation raster map computed by the nearest neighbor method different from the one computed by bilinear interpolation?
Resampling to lower resolution
Resample to lower resolution (30m -> 100m).
First, display the original elevation and land use maps
Step11:  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?
Step12:  Before the next computation, remove all map layers from the Layer Manager
because we don't need to see them anymore.
Step13:  Remove or switch off the land use, elevation and aspect maps.
Converting between vector data types
Convert census blocks polygons to points using their centroids
(useful for interpolating a population density trend surface)
Step14:  Display census boundaries using GUI
Step15:  Convert contour lines to points (useful for computing DEM from contours)
Step16:  Display the "elev_ned_contpts" points vector and zoom-in to very small area
to see the actual points.
Step17:  Convert from vector to raster
Convert vector data to raster for use in raster-based analysis.
First, adjust the computational region to resolution 200m
Step18:  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".
Step19:  Now convert lines in "streets" vector to raster.
Set the resolution to 30m and use speed limit attribute.
Step20:  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
Step21:  Save the displayed map.
In Map Display click button Save display to graphic file,
or use the following.
Step22:  Convert from raster to vector
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.
Step23:  Visually compare the result with streams digitized from airphotos.
Step24:  Save the displayed map (in Map Display click button Save display to graphic file).
Step25:  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
Step26:  Save the displayed map either using GUI or using the following in case
you are working in the command line.
<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)
# 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: Geospatial data models
Resources:
GRASS GIS overview and manual
Recommendations
and tutorial
how to use GUI from the first assignment
Start GRASS GIS
Start GRASS - click on GRASS icon or type
End of explanation
!g.region region=swwake_30m -p
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_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.
Resampling to higher resolution
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.
End of explanation
!d.rast elev_ned_30m
Image(filename="map.png")
Explanation: The -p flag for g.region is used to print the region
we just set.
Then we display the 30m resolution NED elevation raster.
End of explanation
!g.region res=10 -p
!r.resamp.interp elev_ned_30m out=elev_ned10m_nn method=nearest
Explanation: 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.
End of explanation
!d.rast elev_ned10m_nn
Image(filename="map.png")
Explanation: 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:
End of explanation
!r.slope.aspect elevation=elev_ned10m_nn aspect=aspect_ned10m_nn
Explanation: The elevation map "elev_ned10m_nn" looks the same as the original one,
so now check the resampled elevation surface using the aspect map:
End of explanation
!d.rast aspect_ned10m_nn
Image(filename="map.png")
Explanation: Display the resampled map by adding "aspect_ned10m_nn" to Layer Manager
or in command line using:
End of explanation
Image(filename="map.png")
Explanation: To save the map, click in Map Display to on the button
Save display to graphic file" or alternatively,
use the following command:
End of explanation
!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")
Explanation: Now, reinterpolate DEMs using bilinear and bicubic interpolation.
Check the interpolated elevation surfaces using aspect maps.
End of explanation
Image(filename="map.png")
Explanation: 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:
End of explanation
!d.rast elev_ned_30m
!d.rast landuse96_28m
Image(filename="map.png")
Explanation: Why is the aspect of elevation raster map computed by the nearest neighbor method different from the one computed by bilinear interpolation?
Resampling to lower resolution
Resample to lower resolution (30m -> 100m).
First, display the original elevation and land use maps:
End of explanation
!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")
Explanation: 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?
End of explanation
!d.erase
!r.resamp.stats landuse96_28m out=landuse96_100m method=mode
!d.rast landuse96_100m
Image(filename="map.png")
Explanation: Before the next computation, remove all map layers from the Layer Manager
because we don't need to see them anymore.
End of explanation
!v.to.points census_wake2000 type=centroid out=census_centr use=vertex
Explanation: Remove or switch off the land use, elevation and aspect maps.
Converting between vector data types
Convert census blocks polygons to points using their centroids
(useful for interpolating a population density trend surface):
End of explanation
!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")
Explanation: 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.
End of explanation
!v.to.points input=elev_ned10m_cont10m output=elev_ned_contpts type=line use=vertex
Explanation: Convert contour lines to points (useful for computing DEM from contours):
End of explanation
!d.vect elev_ned_contpts co=brown icon=basic/point size=3
Image(filename="map.png")
Explanation: Display the "elev_ned_contpts" points vector and zoom-in to very small area
to see the actual points.
End of explanation
!g.region swwake_30m res=200 -p
Explanation: Convert from vector to raster
Convert vector data to raster for use in raster-based analysis.
First, adjust the computational region to resolution 200m:
End of explanation
!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")
Explanation: 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".
End of explanation
!g.region res=30 -p
!v.to.rast streets_wake out=streets_speed_30m use=attr attrcol=SPEED type=line
Explanation: Now convert lines in "streets" vector to raster.
Set the resolution to 30m and use speed limit attribute.
End of explanation
!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")
Explanation: 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:
End of explanation
Image(filename="map.png")
Explanation: Save the displayed map.
In Map Display click button Save display to graphic file,
or use the following.
End of explanation
!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")
Explanation: Convert from raster to vector
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.
End of explanation
!d.vect streams_derived_t color=blue
!d.vect streams color=red
Image(filename="map.png")
Explanation: Visually compare the result with streams digitized from airphotos.
End of explanation
Image(filename="map.png")
Explanation: Save the displayed map (in Map Display click button Save display to graphic file).
End of explanation
!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")
Explanation: 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.
End of explanation
Image(filename="map.png")
# end the GRASS session
os.remove(rcfile)
Explanation: Save the displayed map either using GUI or using the following in case
you are working in the command line.
End of explanation
<END_TASK>