GRASS GIS

Scripting

Under a LINUX operating system it is possible to create a list, also known as a script, of commands for GRASS GIS to run, this can make long processes a lot quicker. In this section we will go through a basic one command script, and develop to arrays, loops and gui's.

Scripting in LINUX is often conducted in a shell environment, to create your first script open your favorite text editor (emacs, nano, vi etc...). Firstly you need to tell it that it is a BASH script by entering the line #!/bin/bash. All your code will be placed on the following rows.

The first script

You are working on a project and you want to define a region which you can recalculate another time. The normal command is g.region n=0 s=0 w=100 e=100 res=1 -p this will set the region and display it. In a script:

#!/bin/bash
g.region n=100 s=0 w=0 e=100 res=1 -p
exit

You will see that the first lines has #!/bin/bash this is the directory of the BASH shell, to locate yours on LINUX type: which bash

g.region is the GRASS command

and exit exits the script

save this as region1.sh

Next you need to allow the file to be read and executed, navigate to the file and type:

chmod u+x region1.sh

To invoke the script from within GRASS navigate to the file and type ./region1.sh

Looping - Counting to 10

A loop is simply the same as saying: do this then keep doing it until a certain point.

In a script it can be very useful, for instance the simplest example is to count to 10:

for (i=1; i<=10; i++)
do
echo"i";
done

Essentially start at 1, if the current number is less than or equal to 10 then add 1 and do the following, in this case add 1.

It is also possible to count in 10s to 100:

for (i=1; i<=100; i+10)
do
echo"i";
done

Arrays

An array is just a list of values, for instance:

A
B
C
D
E
F

as an array:

\[A\]\
\[B\]\
\[C\]\
\[D\]\
\[E\]\
\[F\]

These are useful for scripting as they can be used in looping, again taking the region as an example, you want to divide your map into 10m squares with a 1m resolution do some analysis then move north to the next 10 metre square. Normally you would enter each one of these commands then move to the next one, like so:

g.region n=10 s=0 e=10 w=0 res=1
g.region n=20 s=10 e=10 w=0 res=1
g.region n=30 s=20 e=10 w=0 res=1
g.region n=40 s=30 e=10 w=0 res=1

With an array and a loop:

[1]=10
[2]=20
[3]=30
[4]=40

for (i=1; i<=4; i++)
do
j=$i-10
g.region n=$i s=$j e=10 w=0 res=1

done

exit

This can get more complex if you start moving the region east and west as well. For the previous example you could have used a loop set at intervals of 10, but Arrays are usefull when using strings (text).

Below is a script I am using to automate the creating of distribution maps

#!/bin/bash
####################################
#Automation of distribution mapping
#uses v.out.ascii r.in.xyz
####################################
#Created 17/03/2010
#By Gary R Nobles
#To create a raster distribution of 
#animal bones form a Neolithic Settlement
####################################

##################
#Array species class
#
s[1]=AMP
s[2]=AVE
s[3]=MAM
s[4]=MOL
s[5]=REP
s[6]=PIS
###########
#End Array#
###########

############
#Parameters#
############
mapin=Bones
col=Count
mapout=Bones_
meth=sum
####################################
#Set the Region
echo "Setting the region"
g.region n=40 s=0 e=40 w=0 res=1
g.region -p
#
#Begin Loop
for (i=1; i<=6; i++)
do

#v.out.ascii pipes to r.in.xyz output raster
v.out.ascii input=$mapin columns=$col where="Class='${s[i]}'" | r.in.xyz input=- output=$mapout${s[i]} method=$meth z=4
echo "${s[i]}"

#stats
r.stats -1 input=$mapout${s[i]}>/home/$mapout${s[i]}.txt 

#need to set null=0
echo "turning null to Zero"
r.mapcalc "$mapout${s[i]}=if( isnull($mapout${s[i]}),0,$mapout${s[i]})"

#export to ArcGIS
r.out.arc input=$mapout${s[i]} output=/home/$mapout${s[i]}.asc dp=8


done
echo "---------------COMPLETE---------------"
exit

This script combines a few GRASS modules to get vector data from a GRASS compatible database (in this case MySQL) turn it into a raster without creating additional files (v.out.ascii output goes directly into the r.in.xyz module). Basic statistics are stored in a file (by using the > symbol). The null vales are turned to 0 and exported for use in ArcGIS.

This shows the power of scripting especially when considering this dataset consisted of 30,000 bone fragments, the time take to process the script was a few seconds at most and it produced 6 distribution maps and basic statistics.