I love GRASS… GIS

But this wasn’t always the case.

GRASS GIS was, for a long time, something I dismissed as ‘too complex’ for my everyday geospatial operations. I formulated any number of excuses to work around the software and could not be convinced it had practical use in my daily work. It was ‘too hard to set-up’, ‘never worked well with QGIS’, and ‘made my scripting processes a nightmare’.

I am here to officially say I’ve been very wrong for a very long time. GRASS GIS is pretty amazing and is a wonderfully easy tool to script.

Recently I scripted a model for extracting river centrelines from high resolution elevation data. During the process I thought how useful it would have been, way back, to have a simple example for setting up an environment and scripting a BASH process using GRASS. So, I built one for myself.

The following example provides the steps for a simple catchment extraction on a piece of LINZ 8m elevation data. The example is only for demonstration of running a basic BASH/GRASS set up with a hydrology command and NOT a demonstration for how to do hydrology using GRASS. The 8m elevation data is not the best data for hydrological extraction; however, these data do provide a nicely sized dataset to use as practice and will give results. Also please note, catchment extraction on a square raster is will not give accurate results at the edges of the raster.

In this example we will:

1. Download a small piece of elevation data from the LINZ Data Service
2. Build a GRASS environment to process these data
3. Build a BASH script to process the catchments
4. Import the elevation into the GRASS environment
5. Perform some basic GRASS operations (fill and watershed)
6. Export raster format for viewing
7. Export the vector catchments to shapefile

The following assumes you are working in a Linux environment and have a basic knowledge of BASH scripting. I tested this process using Ubuntu 18.04 and GRASS 7.4. I have 16GB RAM on my machine.

If you do not all ready have it, you can install GRASS as follows:

$ sudo apt-get install grass-core

If interested, but not necessary, you can install functionality for building GRASS plug-ins:

$ sudo apt-get install grass-dev

Create yourself a directory to work in:

$ mkdir grass_test

Download a piece of elevation data from the LINZ Data Service, place in your directory:

https://data.linz.govt.nz/layer/51768-nz-8m-digital-elevation-model-2012/data/

From the Tiles Table, I downloaded the JM tile in EPSG:2193.

Within your directory you will need to build a ‘PERMANENT’ folder in order for GRASS to do its magic. This folder will be set to operate in the projection of your data. Be sure your data is in the same projection as the environment you built. The data is downloaded in EPSG:2193, New Zealand Transverse Mercator (NZTM), so we set the GRASS environment to work in this projection:

$ grass -c epsg:2193 -e grass_test/GRASS_ENV

‘-c’ will create your directory using epsg:2193 and ‘-e’ will exit once this operation is complete.

Running this command will build a PERMANENT folder in:

grass_test/GRASS_ENV

providing all the bits GRASS needs to operate. Your folder structure will look like:

grass_test/GRASS_ENV/PERMANENT

Have a look inside the folder and see what GRASS built for itself.

The key from now on is to always run your BASH script through this environment. Running your script in this environment will allow you to perform BASH and GRASS commands at the same time.

First, let’s look at the the basic command to run a BASH script using a GRASS environment.

$ grass grass_test/GRASS_ENV/PERMANENT --exec sh grass_test/catchment.sh

The above says, launch GRASS using this environment

grass_test/GRASS_ENV/PERMANENT

and execute, –exec, this script

/grass_test/catchment.sh

Let’s build the BASH script. In this we will:

1. Import the elevation data. GRASS likes to work in its own data formats
2. Set the region for where the operation will be performed. GRASS needs to know where the operation is going to be performed.
3. Perform the operations, r.fill.dir and r.watershed
4. Export a raster to .tif format
5. Export the vector outputs to .shp format

#!/bin/bash

# set base path
outDir=grass_test

# Set raster as variable
raster=${outDir}/JM.tif

# Set a base name for the data. This is used to demonstrate that normal
# BASH commands can be used in this process, along side GRASS
rasterName=$( basename $raster | sed 's/.tif//g' )

# Import raster data
r.in.gdal input=$raster output=$rasterName --overwrite

# Set region. IMPORTANT so GRASS knows where the data is located.
# This region is set for the duration of the following commands
g.region rast=$rasterName

# Fill sinks
fillDEM=${rasterName}_filldem
directionDEM=${rasterName}_directiondem
areasDEM=${rasterName}_areasDEM
r.fill.dir input=$rasterName output=$fillDEM direction=$directionDEM areas=$areasDEM --overwrite

# Export a raster for viewing
areaOut=${outDir}/${rasterName}_areas.tif
r.out.gdal input=$areasDEM output=$areaOut

# Run watershed operation on fill sink raster
threshold=100000
accumulation=${rasterName}_accumulation
drainage=${rasterName}_drainage
stream=${rasterName}_stream
basin=${rasterName}_basin
r.watershed elevation=$fillDEM threshold=$threshold accumulation=$accumulation drainage=$drainage stream=$stream basin=$basin --overwrite

# Convert Basin (watershed) to vector format
basinVect=${rasterName}_basinVect
r.to.vect input=$basin output=$basinVect type=area column=bnum --overwrite

# Export catchment to vector format
basinVectOut=${outDir}/${rasterName}_basinVectOut.shp
v.out.ogr input=$basinVect output=$basinVectOut type=area format=ESRI_Shapefile --overwrite

There you have it. A simple BASH script, set up for running some GRASS commands, running on your LINUX machine.

Putting it all together

Set up you environment, copy the above script, paste it into a text editor, save it as:

grass_test/catchment.sh

Now run

$ grass grass_test/GRASS_ENV/PERMANENT —exec sh grass_test/catchment.sh

Your output catchments shapefile should be similar to the image below

 

Click here for the next step in this process; how to clip your raster using the watershed layer using GDAL and loops.

Let me know if this was helpful or if you would like to see any changes.