GIS/MEA 584:
Mapping and Analysis Using UAS

Using Python and GRASS GIS

• Import lidar points.
• Interpolate with parallelization.
• Process multiple maps at once and generate an HTML report.
• Use matplotlib to show terrain profiles.

• Lake Wheeler GRASS Location; be sure to use separate Mapset for this assignment
• selected LAS tiles from here (see below for the tiles we will need; you can use the provided Python code to download them)

• Your GRASS GIS 7 installation should also include libLAS library which is used by GRASS modules v.in.lidar and r.in.lidar (standalone GRASS GIS for MS Windows, OSGeo4W and Ubuntu packages contain libLAS).
• libLAS installation should include command line tool las2las.
• You will need to edit files in a text editor with Python support, e.g. Gedit or Notepad++, or a Python editor or IDE such as Spyder.

Running the code

On Linux run the script in the GRASS GIS command line (the system one with GRASS GIS written in ASCII) using python command and full path to the .py file, for example:

python /path/to/file.py

On Windows run the script through the GRASS GIS GUI. In Layer Manager main menu, go to File -> Launch script. This will ask you to find the file on your disk. It will ask you a question about paths for addon. If you are not sure, answer No. Then it will run it in the Console in the GUI.

To run the script again, use arrow up key and press enter. Note that the script might have created some maps or files which you might need to remove before running it again. See the end of this assignment for more information on executing Python code.

Downloading the data

import os

# use path which applies to you
# on Windows use something like 'D:\\path\\to\\files'
os.chdir('/path/to/las/files')

import urllib

base_url = 'http://fatra.cnr.ncsu.edu/lidar/'

urllib.urlretrieve(base_url + '0791_005.las', '0791_005.las')
urllib.urlretrieve(base_url + '0791_011.las', '0791_011.las')
urllib.urlretrieve(base_url + '0792_017.las', '0792_017.las')
urllib.urlretrieve(base_url + '0782_020.las', '0782_020.las')
urllib.urlretrieve(base_url + '0781_008.las', '0781_008.las')

Working with command line tools

The LAS files are in different coordinate system, so we need to change it.

First, we check if we are in the directory with the files:

print os.getcwd()

# on Windows use 'D:\\path\\to\\files'
os.chdir('/path/to/las/files')

Now change to coordinate system for all files in the current directory:

import os
from subprocess import call

for f in os.listdir('.'):
    if f.endswith('.las'):
        name, ext = os.path.splitext(f)
        call(['las2las', '--a_srs=EPSG:2264', '--t_srs=EPSG:3358', '-i',
              f, '-o', name + '_spm' + ext])

Import multiple LAS files to GRASS GIS and merging them

g.region n=219637 s=219254 w=636730 e=637193 -p

import os
import grass.script as gs

region = gs.region()

vectors = []
for lidar_file in os.listdir('.'):
    if lidar_file.endswith('_spm.las'):
        bbox = gs.read_command('r.in.lidar', input=lidar_file,
                               output='foo', flags='g').strip()
        bbox = gs.parse_key_val(bbox, vsep=' ', val_type=float)
        if (bbox['n'] < region['s'] or bbox['s'] > region['n']
           or bbox['e'] < region['w'] or bbox['w'] > region['e']):
            gs.info('Skipping tile %s' % lidar_file)
            continue
        name = 'tile_' + lidar_file.rsplit('.', 1)[0]
        vectors.append(name)
        gs.run_command('v.in.lidar', input=lidar_file, output=name,
                       flags='rt', class_filter=2)
gs.run_command('v.patch', input=vectors, output='merged_points',
               flags='b', overwrite=True)
gs.run_command('g.remove', type='vector', name=vectors, flags='f')

As you can see the code above goes over all files in a directory but it creates vector maps only for the files where the data are actually in the current region. The same approach can be used to import data in different format with different modules just for the area of interest (specified using g.region).

DSM interpolation with parallelization

g.region n=219637 s=219254 w=636730 e=637193 res=1 -p

We import GridModule and compute the size of the tiles if we want to run it on 4 cores at once:

import grass.script as gs
from grass.pygrass.modules.grid import GridModule

region = gs.region()
width = region['cols'] / 2 + 1
height = region['rows'] / 2 + 1

grd = GridModule('v.surf.rst', debug=True, width=width, height=height,
                 overlap=10, processes=4, npmin=100,
                 input='merged_points', elevation='mid_pines_small_dem_1m',
                 tension=30, smooth=1, overwrite=True)
grd.run()

Then run it in parallel (debug=False) and we can see the time needed reduced significantly. On Windows you need to enclose the entire code into main function to work properly:

import grass.script as gs
from grass.pygrass.modules.grid import GridModule

def main():
    region = gs.region()
    width = region['cols'] / 2 + 1
    height = region['rows'] / 2 + 1
    grd = GridModule('v.surf.rst', debug=False, width=width, height=height,
                    overlap=10, processes=4, npmin=100,
                    input='merged_points', elevation='mid_pines_small_dem_1m',
                    tension=30, smooth=1, overwrite=True)
    grd.run()

if __name__ == '__main__':
    main()

Note that parallelization is not always straightforward. Here we call v.surf.rst to interpolate raster from vector points in different parts of the area. These different calls run at once at different processor cores and they don't communicate with each other. Some modules, such as r.sim.water where the water needs to continuously flow through the whole landscape, cannot be called in parallel in this way. Also note that parallelization often comes with some cost. Here the different tiles must be patched together at the end. However, we actually get the speed improvement because running the interpolating at multiple cores at once is much greater improvement than the slow down cause by patching.

Working with large number of maps

g.region raster=2015_06_20_pix4d_11GCP_dsm

import grass.script as gs

# adjust the names according to what you used in the previous assignments
rasters = ['2015_06_20_pix4d_11GCP_dsm', '2015_06_20_DSM_Trimble_11GCP',
           '2015_06_20_DSM_agi_11GCP']
for raster_a in rasters:
    for raster_b in rasters:
        # the two for loops create combination of each map with each map
        # including the map itself and we don't want to create a difference
        # of map with itself, so we skip these combinations
        if raster_a == raster_b:
            continue
        # create the name for the difference raster
        difference = "diff_{a}_{b}".format(a=raster_a.replace('@', '_'),
                                           b=raster_b.replace('@', '_'))
        # compute difference between the rasters using r.mapcalc
        gs.mapcalc('{diff} = {a} - {b}'.format(
            a=raster_a, b=raster_b, diff=difference))
        # set color table to the resulting map
        gs.run_command('r.colors', map=difference, color='differences')

g.remove type=raster pattern="diff_*" -f

Generating a report

First, we get a list of raster maps we are interested in. We use g.list wrapper with pattern option. Then we compute univariate statistics and we also compute shaded relief for each raster map. Finally, we render the raster map together with shaded relief using d.* modules to a PNG file. Note that we are using d.shade to avoid creating another raster which would combine the original raster and its shaded relief.

import grass.script as gs

# get list of rasters we are interested in using a search pattern
rasters = gs.list_strings('raster', pattern='2015_06_*')
# initialize a dictionary to hold statistics for each raster
stats = {}
# iterate through the list of rasters
for raster in rasters:
    # save univariate statistics for raster to dictionary
    stats[raster] = gs.parse_command('r.univar', map=raster, flags='g')
    # compute shaded relief, use name of the original raster including mapset
    shaded_relief = raster.replace('@', '_') + '_shade'
    gs.run_command('r.relief', input=raster, output=shaded_relief,
                   overwrite=True)
    gs.run_command('r.colors', map=raster, color='elevation')
    # render the raster (geographical extent follows current region)
    gs.run_command('d.mon', start='cairo', output=raster + '.png',
                   width=400, height=400, overwrite=True)
    gs.run_command('d.shade', shade=shaded_relief, color=raster)
    gs.run_command('d.mon', stop='cairo')

# a template for an HTML report for one raster
template = """<h2>Raster map {name}</h2>
<h3>Statistics</h3>
<table>
<tr><td>Mean</td><td>{mean}</td>
<tr><td>Variance</td><td>{var}</td>
</table>
<h3>Image</h3>
<img src="{name}.png">
"""

# write to a file using a template
with open('report.html', 'w') as output:
        for raster in sorted(rasters):
            stat = stats[raster]
            output.write(template.format(
                name=raster, mean=stat['mean'], var=stat['variance']))

g.remove type=raster pattern="*_shade" -f

Terrain profiles

Compute a profile for multiple rasters and plot them in matplotlib:

import grass.script as gs
# on Windows, uncomment the two following lines
# import matplotlib  # required by Windows
# matplotlib.use('wx')  # required by Windows
import matplotlib.pyplot as plt

# adjust the names according to what you used in the previous assignments
rasters = ['2015_10_06_DSM_agi_8GCPs', 'mid_pines_lidar2013_dem']
coordinates = [637160.446919, 219373.236976,
               637105.693795, 219392.416036,
               637072.439779, 219400.613538]
profiles = {}
distance = []
for raster in rasters:
    profile = gs.read_command('r.profile', input=raster,
                              coordinates=coordinates, quiet=True).strip()
    # parse output
    if not distance:
        for line in profile.splitlines():
            distance.append(float(line.split()[0]))
    profile_elev = []
    for line in profile.splitlines():
        profile_elev.append(float(line.split()[-1]))
    profiles[raster] = profile_elev

for raster in profiles:
    plt.plot(distance, profiles[raster], label=raster)
    plt.legend(loc=0)
plt.show()

import numpy as np

stats = {}
for raster in profiles:
    profile = np.array(profiles[raster])
    maxim = np.max(profile)
    minim = np.min(profile)
    mean = np.mean(profile)
    stddev = np.std(profile)
    median = np.median(profile)
    roughness = np.std(np.diff(profile))
    stats[raster] = (minim, maxim, mean, stddev, median, roughness)

with open('profile_stats.csv', 'w') as f:
    f.write(','.join(['name', 'minim', 'maxim', 'mean',
                      'stddev', 'median', 'roughness']))
    f.write('\n')
    for raster in profiles:
        f.write(raster + ',')
        f.write(','.join([str(i) for i in stats[raster]]))
        f.write('\n')

Tasks and deliverables

• Try the provided code snippets. Modify them as you wish.
• Run the script for creating profiles several times with different coordinates.
• Extent the generated report by computing viewshed from a selected point.
• Deliver the generated report as PDF.
  • To get the PDF, print the web page as PDF in a web browser (it should be also possible to copy the HTML through the clipboard to a word processor like LibreOffice Writer).
  • Alternatively, generate report in LaTeX, or another markup, and deliver generated PDF.
  • In any case, include also the source text file (but you don't have to include the individual images).
• Alternatively, you can run your own code or use your own data.

Downloading all files linked on a website

import urllib
import re

# note the slash at the end
base_url = 'http://fatra.cnr.ncsu.edu/lidar/'

tile_file = 'tile_list.html'
# download the file (in one step)
urllib.urlretrieve(base_url, filename=tile_file)

# open the file and read it line by line
# we expect to get HTML with a list of filenames in the given directory
with open(tile_file) as html:
    for line in html:
        # search for the links (a tags with href attribute)
        # we expect one link per line in the HTML source code
        match = re.search(r'<a href="(.*?)">', line)
        if not match:
            continue
        # get the content of the href attribute captured before
        filename = match.group(1)
        # download only files we are interested in
        if not filename.endswith('.las'):
            continue
        # we expect the link to be relative to our base URL
        urllib.urlretrieve(base_url + filename, filename=filename)

Script with parameters

import sys

if len(sys.argv) != 3:
    sys.exit("Please, provide two parameters")
a = float(sys.argv[1])
b = float(sys.argv[2])
c = a + b
print "{a} + {b} = {c}".format(a=a, b=b, c=c)

GRASS GIS provides its own way of handling parameters which is specialized for geospatial data and standardized across GRASS modules. Interface defined like this can be used to create command line help, documentation or GUI. However, even when used just in a script, first advantage being the usage of key-value syntax in the command line as opposed to relying on the order of parameters. In comparison to Python packages for command line parsing, GRASS GIS provides predefined types related to various values and data types used in GRASS GIS. Read more about it in the documentation.

A proper structure of a Python script

#!/usr/bin/env python

import sys

def main():
    a = 3
    b = 5
    print a + b

if __name__ == "__main__":
    main()

More ways of executing Python code with GRASS GIS

On Linux (and Mac), all scripts should contain a shebang line which tell the system which interpreter to use. The shebang is a first line in a file and starts with #!. For your Python scripts, use:

#!/usr/bin/env python

chmod u+x /path/to/file.py

/path/to/file.py

If you want to run the code which uses GRASS GIS packages from your Python editor you can start the editor from the system command line which opens with GRASS GIS. Then the editor is part of the GRASS GIS environment and GRASS Python packages will be accessible and GRASS modules will recognize the current Location and Mapset.

If you want to start GRASS GIS without GUI to just run your Python scripts, your can run the grass command in terminal on Linux or using the GRASS GIS Command Line application on Windows.

If you want to use GRASS Python packages and GRASS modules without starting the GRASS GIS application at all, you can do so by setting up the path to GRASS Python packages and doing an initialization procedure to setup runtime environment for modules and also connect to a given GRASS Database, Location and Mapset. Please refer to the documentation for details.

Alternatively, you can set your system environment to be the proper environment GRASS modules will work in without a need to actually start GRASS GIS. To run something you would need to just set GRASS Database, Location and Mapset. There is a reason why GRASS GIS is not doing that, although there is a lot of applications which are invading your system in this way or force you to do it, especially on Windows. If you wish GRASS GIS would be the same, use the documentation or the wiki as a guide. This solution should be fine in some isolated environment such as virtual machine. Less invasive version of this is to set just Python packages path (PYTHONPATH) and GRASS installation directory (GISBASE) so that you can import GRASS Python packages and leave the rest on the GRASS setup functions. Note this can still result in problems with different versions of Python on your system (not generally an issue on Linux).

Note that last two ways are more intended for application developers using Python and GRASS GIS. GIS professionals generally more benefit from running the scripts from within GRASS GIS as they can review the results right-away and such scripts can be easily turned into real GRASS modules with a consistent user interface.