GIS/MEA582: Geospatial Modeling and Analysis

Spatial interpolation and approximation I: methods

Resources:

GRASS GIS overview and manual
Recommendations and tutorial how to use GUI from the first assignment

Start GRASS GIS

Create a new mapset (called e.g. HW_interpolation_1) 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:

cd

Compute Voronoi polygons

Display the polygons with centroids.
Find the column name where z is stored and convert the polygons to raster.
Compute aspect to evaluate the surface geometry.

g.region rural_1m -p
d.erase
v.voronoi elev_lid792_randpts output=elev_vor
d.vect elev_vor size=1 -c type=area,centroid
v.info -c elev_vor
v.to.rast elev_vor output=elev_vor_1m attrcolumn=value use=attr
r.colors elev_vor_1m color=elevation
r.slope.aspect elevation=elev_vor_1m aspect=asp_vor_1m

Display the resulting elevation map and aspect map.

d.erase
d.rast elev_vor_1m
d.rast asp_vor_1m
d.vect elev_lid792_randpts size=1 color=red
d.out.file aspect_voronoi

Optionally, you can view it in 3D perspective (switch off all layers except for elev_vor_1m and switch to 3D view).

Interpolation using IDW

Set region and resolution, find a column name where z is stored.
Interpolate DEM using IDW, check the result using aspect.

g.region rural_1m -p
v.info -c elev_lid792_randpts
v.surf.idw elev_lid792_randpts output=elev_idw_1m column=value
r.colors elev_idw_1m color=elevation
r.slope.aspect elevation=elev_idw_1m aspect=asp_idw_1m
d.erase
d.rast elev_idw_1m
d.rast asp_idw_1m
d.out.file aspect_idw

Design experiment that elucidates the impact of IDW parameters on the surface, focus on the impact of:

exponent e.g., power=0.5, 1, 5 (2 is the default)
number of neighboring points e.g., npoint=1, 5, 20, 60 (12 is the default)

Include selected images (e.g. hillshade or aspect) and relevant stats (e.g., mean, min, max from r.univar, histogram) that highlight the differences in the resulting surfaces into your report.

Check the surface interpolated with default parameters using 3D view.
Do not forget to switch off everything except for the interpolated elevations and set fine resolution to 1.
You can use constant color for the surface to highlight its structure.
Save an image for your report.

Compute DEM from contours

Compute DEM from contours using linear interpolation between isolines:

g.region rural_1m -p
v.to.rast elev_lid792_cont1m output=el_lid792_cont1m attrcolumn=level use=attr
r.surf.contour el_lid792_cont1m output=el_rcont
r.colors el_rcont color=elevation

Check the result using a 2D aspect map or view el_rcont in 3D.
In 3D set view from SE and light from NW to reveal subtle geometry.

r.slope.aspect elevation=el_rcont aspect=asp_rcont
d.rast el_rcont
d.rast asp_rcont
d.vect elev_lid792_cont1m col=white
d.out.file asp_rcont

Optional: create TIN model

Convert z-value stored as attribute "value" to z-coordinate.
Compute TIN:

v.to.3d elev_lid792_randpts output=elev_lid792_randpts3d column=value
v.delaunay elev_lid792_randpts3d output=elev_rand_tin
v.colors elev_rand_tin use=z color=elevation
r.mapcalc "level90 = 90"

Visualize the TIN as 3D vector data:
Keep only "level90" and "elev_rand_tin" switched on (remove or uncheck everything else).
Switch the view from 2D to 3D. Go to Data > Vector and unckeck Show vector points. In Vector lines, change color from black to gray and set Display from on surface to as 3D.

Optional: Use Python to create the data for IDW comparison

A part of the first interpolation assignment (4A) for GRASS GIS is computing IDW with different parameters. The task is to compare the different results in any way you think is appropriate and comment on the results. To comment on the results, you need to create them. You can simply execute the module with different parameters and then compute all the subsequent analyses which will help you determine what are the properties of the surfaces generated by IDW with different settings. However, you can also use Python to compute this data if you know Python at least a little bit.

The simplest way to execute the Python code which uses GRASS GIS packages is to use Simple Python editor integrated in GRASS GIS accessible from the toolbar or the Python tab in the Layer Manager. Another option is to write the Python code in your favorite plain text editor like Notepad++ (note that Python editors are plain text editors). Then run the script in GRASS GIS using the main menu File -> Launch script.

The way to call GRASS modules in Python is very similar to what we use in the assignments and it is described in the documentation. The page also contains several examples how to write simple and more advanced scripts. Here are some examples relevant to the assignment.

Interpolation using different number of points (modify the list of numbers of points):

import grass.script as gs

for npoints in [1, 20]:
    name = 'elev_idw_1m_npoints_{}'.format(npoints)
    stats = gs.parse_command('v.surf.idw', input='elev_lid792_randpts',
                             output=name, column='value', npoints=npoints)

Computing statistics but showing only some for different number of points (you can combine the code with the code above):

import grass.script as gs

for npoints in [1, 20]:
    name = 'elev_idw_1m_npoints_{}'.format(npoints)
    print("\n\n")
    print(name)
    print(len(name) * "=")
    stats = gs.parse_command('r.univar', map=name, flags='eg')

    print(stats['min'], stats['max'])

Now write your own code to compute IDW interpolation for different values of power parameter (e.g., 0.5, 1, 2, 5) by modifying the previous code. Once you have the elev_idw_1m_power_05, .... computed you can set the color table and compute shaded relief for your outputs:

import grass.script as gs

for power in [0.5, 1, 2, 5]:
    gs.run_command('r.colors',
        map='elev_idw_1m_power_{}'.format(power),
        color='elevation')
    gs.run_command('r.relief',
        input='elev_idw_1m_power_{}'.format(power),
        output='elev_idw_1m_power_{}_relief'.format(power))
    gs.run_command('r.shade',
        color='elev_idw_1m_power_{}'.format(power),
        shade='elev_idw_1m_power_{}_relief'.format(power),
        output='elev_idw_1m_power_{}_shaded'.format(power))

Creating a PNG image with histogram for changing power:

import grass.script as gs

for power in [0.5, 1, 2, 5]:
    gs.run_command('d.mon', start='cairo',
        output='elev_idw_1m_power_{}_histogram.png'.format(power))
    gs.run_command('d.histogram',
        map='elev_idw_1m_power_{}'.format(power))
    gs.run_command('d.mon', stop='cairo')

Here are two commands often used when using the scripts. First is setting the computational. We can do that in a script, but it better and more general to do it before executing the script:

g.region region=rural_1m

When we want to run the script again, we first need to remove the created raster maps:

g.remove type=raster pattern="elev_idw_1m_npoints_*"

In case you don't know anything about Python scripting but you still want to try something this might be a good start together with some (free) courses at Codecademy. To learn more about using Python in GRASS GIS, see the introduction to the grass.script package.