Geographic Resources Analysis Support System (GRASS): More Than a Mapping Tool

More-complex analysis that combines various types of data and examines their spatial relationships still requires considerable expertise and extensive GIS capabilities. Among many open source GIS tools, Geographic Resources Analysis Support System (GRASS) is one of the few systems that provides a comprehensive set of tools for working with georeferenced data.

Originally designed at the U.S. Army Construction Engineering Research Laboratory as land management support software for military installations, GRASS has evolved into a powerful, multipurpose GIS for geospatial data management, analysis, modeling, and visualization. The system was released under GNU GPL in 1999 and is currently developed by the international GRASS Development Team, coordinated at the Centro per la Ricerca Scientifica e Tecnologica (ITC-irst), in Italy. The binaries and the source code can be downloaded from the GRASS Web site or any of its 25 worldwide mirror sites.

GRASS stores the geometric and attribute data in its own database, which is adequate for average GIS users. More sophisticated applications benefit greatly from the ability to connect GRASS projects to an external centralized database such as an Oracle database.

Install GRASS with Oracle Support

Compiling GRASS from its sources requires installation of several supporting libraries. Most are standard for various Linux distributions, so usually only the Cartographic Projections Library (PROJ) and Geospatial Data Abstraction Library (GDAL)/OGR Simple Features Library need to be downloaded and installed separately. Prebuilt binaries can be used for PROJ, but GDAL/OGR need to be compiled from sources to include Oracle support. GDAL/OGR is a translator library that handles most raster and vector formats used for geospatial data. GRASS uses GDAL/OGR as an import/export back end for maximum interoperability. Several geospatial applications, such as Google Earth, also use this library for their data handling.

To access and write data sets stored in Oracle Database from within GRASS, you need to build OGR using the Oracle Call Interface (OCI), which, by default, is not compiled in OGR. When building GDAL/OGR, you need to define the path to the full Oracle client library, which is available with a standard Oracle installation. Make sure to download the full client library, because also some header files are needed for compiling.

Downloading and building GDAL/OGR with Oracle support is quite straightforward (you may want to select this or a later version):

wget ftp://ftp.gdal.org/gdal/gdal-1.3.1.tar.gz 
tar xvfz gdal-1.3.1.tar.gz cd gdal-1.3.1 
./configure 

--with-oci=/path/to/oci-library 
make 
make install

If the build was successful, you can use the ogrinfo command to list supported vector formats, as in the shortened example below:

ogrinfo --formats 
Loaded OGR Format Drivers:
-> "ESRI Shapefile" (read/write) 
[...]
 -> "SQLite" (read/write) 
 -> "ODBC" (read/write) 
 -> "OCI" (read/write) 
 -> "OGDI" (readonly) 
 -> "PostgreSQL" (read/write) 
[...]

Now the OCI bindings are compiled into OGR and you can proceed to build GRASS by using the just installed GDAL/OGR library. To follow the latest developments, we suggest to download the weekly CVS source code snapshot of GRASS 6.1, which includes support for using data sets stored in Oracle Database. Build it as follows.

wget http://grass.itc.it/grass61/source/snapshot/grass-6.1.cvs_src_snapshot_[date].tar.gz 
tar xvfz grass-6.1.cvs_src_snapshot_[date].tar.gz 
cd grass-6.1.cvs 
./configure \ 
[...] --with-gdal=/usr/local/bin/gdal-config \ 
--with-proj
make 
make install

Everything is compiled now, but you still need to define the path where your Oracle database is installed, by defining the related environment variable ORACLE_HOME:

export ORACLE_HOME=/path/to/oracle-installation

If you plan to work with Oracle database on a daily basis, you can define the environment variable ORACLE_HOME within your environment permanently:

echo "export ORACLE_HOME=/path/to/oracle-installation" >> ~/.bashrc

You can now test the Oracle bindings, by listing all spatially enabled layers of the database instance, using

ogrinfo -summary -ro "OCI:user/pw@db_instance"

At this point, the GRASS GIS and Oracle database are ready to be used together for geospatial data processing and analysis.

Getting Started with an Existing Data Set

The easiest way to learn GRASS is to start with an existing, ready-to-use data set. Several are available at the GRASS Web site. This article uses the data set for an area in South Dakota, USA, called Spearfish. First create a subdirectory in which you will store all your GRASS data (you can name it grassdata), and then download and extract the Spearfish data set into it:

cd $HOME
mkdir grassdata
cd grassdata
wget http://grass.itc.it/sampledata/spearfish_grass60data.tar.gz
tar xvfz spearfish_grass60data.tar.gz

Structure of GRASS database and commands. After extracting the Spearfish data set, you can see the structure of the GRASS database (see Figure 1). GRASS data is stored in a directory referred to as a database (also called GISDBASE); in your case, it will be the grassdata directory you have just created above. Within this database, the projects are organized by locations (subdirectories of the database): The downloaded data set contains the location spearfish60. It is important to know that each location is defined by its coordinate system, map projection, and geographical boundaries.

Figure 1: Structure of a GRASS 6 database

Each location can have several “mapsets” that are used to subdivide the project into different topics, subregions, or workspaces for individual team members. Each mapset includes subdirectories for raster, vector, and attribute data, plus a working spatial extent definition file, named WIND. These subdirectories and files are hidden from the user while working in GRASS. When defining a new location, GRASS automatically creates a special mapset called PERMANENT, which is used to store the default spatial extent, the coordinate system definitions and which may be used to also store geodata for the team.

The mapset concept provides a convenient environment for teamwork, by allowing any number of users to work in a single location simultaneously by using their mapsets. Users can read maps from other mapsets only if permissions are granted, but they can never write to other mapsets. Maps to be shared by all team members should be stored in the PERMANENT mapset. This simple scheme permits easy management of large GIS projects, even with remotely mounted network file systems such as NFS.

GRASS has more than 350 modules, so it is useful to get familiar with its naming convention—it is very simple, as shown in Table 1. The prefix indicates the functional group (type of operation the command performs), and the word after the dot describes what the command does or what type of data it works with.

Table 1: Command structure of GRASS module names

Start GRASS and display data. Now you are ready to launch GRASS 6.1; depending on the local installation, you can do it from either the menu or a terminal window, by typing

grass61

A graphical user interface should open, as shown in Figure 2. Enter the path to the database into the first field; in our case, it is /home/user1/grassdata. Then select spearfish60 as your location and user1 as your mapset.

Figure 2: GRASS 6 startup screen with selection of database, location, and mapset

Then click Enter GRASS at the bottom left. After some basic information about the GRASS version and how to access help:

Welcome to GRASS 6.1.cvs (2006)
GRASS homepage:                                 http://grass.itc.it/
This version running through:                           Bash Shell (/bin/bash)
Help is available with the command:                     g.manual -i
See the license terms with:                     g.version -c
If required, restart the graphical user interface with:d.m &
When ready to quit enter:                               exit
GRASS 6.1.cvs (spearfish60):~ >

you get the GRASS shell prompt and the GUI display manager (d.m) opens (see Figure 3):

Figure 3: GRASS GIS manager (d.m)

If you know the path to your database, location, and mapset, you can start GRASS directly from the command line:

grass61 /home/user1/grassdata/spearfish60/user1

Although GRASS can be used with several types of GUIs, such as the default tcltkgrass with d.m display manager (Figure 3) or a separately installed QGIS (Figure 4, www.qgis.org), command-line is used in this article for more-compact text.

Figure 4: QGIS interface with GRASS data visualized

Getting Help

To get an overview of all available commands and to access their manual pages in a HTML browser, type

g.manual -i

If you know the name of the command for which you need help, a short description of the command usage can be displayed if you use the -help parameter. For example:

d.rast -help

To find out what raster and vector data is accessible in your location; what coordinate system, cartographic projection, and units are used; and what your spatial extent (coordinates for your working region) and the current raster resolution type are, type

g.list rast 
g.list vect 
g.proj -w 
g.region -p

If you don’t know what a projection is and what all the coordinates mean, don’t worry about it now. You won’t need it for the examples that are used in this article; however, you will have to learn at least some basics if you want (or need) to do any serious work with geospatial data.

Display Data

The GRASS d.m or QGIS provide convenient GUIs for displaying complex sets of raster and vector data, but for our introductory examples, using a simple command-line display is sufficient:

# set the current spatial extent to the extent of the elevation raster map 
g.region rast=elevation.dem -p


# start a graphic-monitor and display raster digital elevation model (DEM), 
# and roads (vector line data) 
d.mon x0 
d.rast elevation.dem d.vect roads

You can also explore the data projected in 3-D space (Figure 5), by using the GRASS visualization module, NVIZ. To learn more about the use of the visualization GUI, select NVIZ help from the help menu, in the upper right corner of the NVIZ window that opens after you type

nviz elevation=elevation.dem vector=roads

Figure 5: Visualization in NVIZ viewer

Analyze Relationships

Geospatial analysis is one of the main strengths of GRASS GIS, and you can start exploring its capabilities by computing the basic terrain parameters, such as slope and aspect, and extracting stream networks and basin boundaries from a digital elevation model (DEM). These parameters are used in many applications, especially in land management, and they provide helpful information for personal use—for example, when planning a hiking or bicycle trip, you can find how many streams you will cross and how steep different segments of your trip will be. Our example uses terrain analysis combined with land ownership data to find landowners who control the land along the streams. Such information is often needed for conservation, water quality protection, or water resources research.

# compute and display slope and aspect (terrain orientation) maps 
r.slope.aspect elevation=elevation.dem slope=myslope aspect=myaspect 
d.rast myslope 
d.rast myaspect

# extract and display stream network and basins 
r.watershed elevation=elevation.dem basin=mybasins \
stream=mystreams threshold=1000 
d.rast mystreams 
d.rast mybasins

You can create a more sophisticated representation of the extracted streams, by thinning their raster representation to a single cell width; converting them to vector lines; and deriving some of their attributes, such as stream segment lengths:

# convert the raster stream map to vector line representation 
# and create an attribute table to store stream length attributes 
r.thin in=mystreams out=mystreams_thinned 
r.to.vect -s in=mystreams_thinned out=streams_vector 
v.db.addtable map=streams_vector table=streams_vector col="length int" 
d.vect streams_vector width=2

# create a new vector map with new categories added (vector IDs) 
# in this step an ID and table row is added to each vector line: 
v.category in=streams_vector out=streams_vector2 op=add 
v.to.db map=streams_vector2 op=cat column=cat 

# verify that all IDs appear: 
v.db.select streams_vector2

# populate stream segment length to DB and display the length values on the map 
v.to.db map=streams_vector2 option=length column=length units=me 
v.db.select streams_vector2 
d.vect streams_vector2 
d.vect streams_vector2 display=attr attrcol=length

# optionally you can view the new basin and stream maps draped over terrain 
# in 3D space 

nviz elevation.dem color=mybasins vect=streams_vector2

# clean up 
g.copy vect=streams_vector2,streams_final 
g.mremove vect="streams_vector*"

# store resulting stream-vectors in Oracle 
v.out.ogr in=streams_final format=OCI dsn="OCI:user/pw@db_instance" \
olayer=streams_oracle

# verification 
ogrinfo -summary -ro "OCI:user/pw@db_instance" streams_oracle

In the next step, you can combine the stream segments with the landowner data to find the owners with land adjacent to the stream segments and accumulate the total stream length per landowner.

# calculate intersection between landowners and stream-segments resulting in
# a new data set  streams_fields
v.overlay ain=streams_final atype=line bin=fields out=streams_fields operator=and

# add another column for segment-length per landowner
v.db.addcol streams_fields col="length int"

# add length to the newly created segments per fields
v.to.db map=streams_fields option=length column=length units=me

# add another column for accumulated segment-length per landowner
v.db.addcol streams_fields col="acclength int"

# export your intersected data set into Oracle
v.out.ogr in=streams_fields format=OCI dsn="OCI:user/pw@db_instance" \
  olayer=streams_fields_oracle

The recently added OGR-binding in GRASS still needs some fixing to do the export properly, so you need to manually adjust the feature-ID-column using the following SQL-Statement:

echo "UPDATE streams_fields_oracle SET \"OGR_FID\" = \"cat\";" \
| sqlplus user/pw@db_instance

Now you can accumulate the stream segment lengths for each landowner. Since you have exported the data set into Oracle, do it directly in Oracle using the sqlplus commandline tool.

# grab the minimum/maximum of category numbers from fields for updating the 
# segments
v.univar fields col=cat

# update the accumulated length for each landowner
for i in `seq 1 63`; do \
  echo "UPDATE streams_fields_oracle SET \"acclength\" = \
  (SELECT sum(\"length\") FROM streams_fields_oracle WHERE \"b_cat\"= $i) \
  WHERE \"cat\" = (SELECT \"cat\" FROM streams_fields_oracle WHERE \"b_cat\"= $i \
  AND rownum = 1);" >> /tmp/acclength.sql ; \
done

# update SQL in oracle
cat /tmp/acclength.sql | sqlplus user/pw@db_instance

# set the accumulated length to all related stream segments per landowner
for i in `seq 1 63`; do \
  echo "UPDATE streams_fields_oracle SET \"acclength\" = (SELECT \
  \"acclength\" FROM streams_fields_oracle WHERE \"b_cat\"=$i AND rownum = 1) \
  WHERE \"acclength\" = 0 AND \"b_cat\" = $i;" >> /tmp/acclength_all.sql ;\
done

# update SQL for all stream segments
cat /tmp/acclength_all.sql | sqlplus user/pw@db_instance

After updating the attribute columns in Oracle Database, you can work with the data directly in GRASS, using the v.external command, which “links” external data sets into GRASS. You can now display the results, using the following commands:

# re-link your modified data from Oracle back into GRASS 
v.external dsn=OCI:user/pw@db_instance out=streams_fields_external \
layer=streams_fields_oracle


# display the streams with over 2000m of accumulated length for a landowner
# in red with less than 2000m length in green along with the property
# boundaries (field) and basins
d.rast mybasins -o
d.vect fields type=boundary,centroid att=label lsize=10 xref=left \
  yref=bottom display=shape,attr type=point,line,boundary,centroid size=1 \
  layer=1 color=0:0:0 lcolor=0:0:0 fcolor=170:170:170 llayer=1 icon=basic/x 
d.vect streams_fields_external width=3 where="\"acclength\" >= 2000 AND \
  \"b_label\" != 'Black Hills Natl. Forest' " col=red
d.vect streams_fields_external width=3 where="\"acclength\" < 2000 OR \
  \"b_label\" = 'Black Hills Natl. Forest' " col=green

# display table
v.db.select streams_fields_external col=\"acclength\",\"b_label\" \
  where="\"acclength\" >= 2000" |sort -nu

# cleanup
rm -f /tmp/acclength*

# exit GRASS
exit

The results show that in the Spearfish area, the highest accumulated stream length is located within the Black Hills National Forest, providing most of the streams with adequate protection. However, there are several landowners who can have significant impact on water quality (Figure 6). You can try to learn more about their possible impact by combining the land ownership map fields with the land cover map landcover.30m and the extracted streams, to identify those landowners who should be encouraged to put the agricultural land along the stream into the conservation program.

Figure 6: Subsection of the resulting map showing the relation between the basins, streams, and land ownership

If you have your own spatial data stored in the Oracle database, you should now be able to read and import it into GRASS and output the results of geospatial analysis back into Oracle. Your data will most likely be in a different coordinate system than our spearfish60 location (UTM), so you will have to learn how to create a new location in GRASS. A quick description is in the document “GRASS 6 in a Nutshell”.

A Parting Note

This brief introduction to GRASS GIS has explored the system from the user’s side and demonstrated a simple spatial analysis. GRASS, as open source software, provides full access to its code; it is therefore possible to modify the modules or add new ones. The code is written in C, and an introduction to programming in GRASS would require a separate article. If you are interested in development, you can find the GRASS Programmer’s Manual on the GRASS Web site (PDF, HTML); it is updated weekly, and the code is documented inline in source code files (“doxygen”-based). The GRASS GIS development mailing list is a good place to look for help.

Helena Mitasova is a researcher at the Department of Marine, Earth and Atmospheric Sciences at North Carolina State University, and an active member of GRASS Development Team and Open Source Geospatial Foundation.

Markus Neteler received his MSc degree in Physical Geography and Landscape Ecology from the University of Hanover in Germany in 1999. He is researcher at ITC-irst and CEA, Trento, Italy, and is co-author of two books on GRASS.

Stephan Holl received his MSc deegree in Physical Geography and Landscape Ecology from University of Hannover in Germany in 2003. He works as a project manager in several projects dealing with Open Source GIS software at GDF Hannover, Germany.