GRASS Scripting: Contours and Reliefs

grass is a library of commands for manipulating, analyzing and displaying cartographic information. While the functionality itself has evolved now for many years, it was scattered over different organisations, programming languages and a zoo of different map formats.

grass' biggest achievement is to have brought all this functionality into a common framework and also to provide an interactive Tk application which allows amateurs like me to experiment with amazing features while remaining mostly ignorant about formats and other cartographic idiosyncracies.

The interactive application is really only a more or less thin layer on top of the set of grass commands, so whenever one triggers any action in the GUI, one or more grass commands are executed underneath. All of which is recorded in a log file. And that becomes incredibly useful once you start to automate things to run unattended in batch mode.

Like with many other open source projects (guilty myself), the documentation is brutally incomplete, self-referential, inconsistent and takes a lot of educated guesswork. So maybe the following is of help.

Feeding the Cranky Kid

Before doing anything mappish, it is necessary to prime your database. Grass (tears slowly fill my eyes) is using a file-based approach, although quite some work is invested to (somehow) integrate vector information from relational databases. To build this file database, best is to download the spearfish example and have it extracted somewhere.

Assuming that the directory is called spearfish60, then with

cat test.gr | grass -text spearfish60/PERMANENT/

you can invoke grass in batch mode. The file test.gr is piped into grass. It contains grass commands for it, but not quite, though, because grass itself is actually another UNIX shell. One, which extends the existing one with the PATH for the grass commands and a huuuuuuge list of other state variables. This will bite when you plan to use grass in a concurrent environment, such as a web server (tears running down my cheeks).

Inside the file database you then have to select a part in which you choose to work in, PERMANENT in our case. There is already stuff in there, but that does not bother us.

Setting the Stage

Given the above, the commands to be put into test.gr can be any shell command. We here are mainly interested in grass-related ones only, but it should be clear from the above that the syntax these follow is close to UNIX shell syntax (tears streaming down my face).

For instance, a simple

g.version

will welcome us with a version information. What I found useful when testing, is to coax grass to allow files to be overwritten automatically when you repeat the same commands over and over:

export GRASS_OVERWRITE=1

Otherwise it will refuse to do so.

Elevation Data Input

While the spearfish demo database contains quite a number of sample data, you might want to create some of your own. One kind of data set can carry elevation data. Such can be easily produced inside an ASCII file:

north: 4928000
south: 4914020
east: 609000
west: 590010
rows: 466
cols: 633
 * * * * 1112 1111 1111 1110 1111 1113 1113 ...
...

The header block first contains world coordinates of the geographical boundaries of the elevation data. I borrowed them from an existing map (I believe this is the area around San Francisco). The rows and cols attributes detail how many data rows will follow and how many colums each of them holds. Then follows the data itself whereby a * indicates to grass that the value is missing.

To load this data from the file rhabarba.elev into a raster map named Rhabarba within the grass database, you have to launch the command:

r.in.ascii input=rhabarba.elev output=Rhabarba

If you want to avoid that grass is driving you insane, it is a clever move to focus your current region onto this map:

g.region rast=Rhabarba

Otherwise some other region might be the global default and your display remains frustratingly empty.

Outputting

If you wanted to display your raster map, then the d.rast command is your choice:

d.rast map=Rhabarba -o

The -o switch instructs the program to look at non-NULL values only. That would show something like this:

1

but since we want to run this in batch-mode, it is better to display it in a virtual window, which can be later output in, say, PNG. So before running d.rast we insert:

d.mon start=png1

Then grass will write everything into a file map.png inconveniently placed into the current directory (whinge).

Better control I found to have with d.out.file. With it you can specify the format (PPM) and also the file location:

d.out.file format=ppm output=/tmp/Rhabarba.elev

Absurdedly, d.out.file appends another file extension, in this case .ppm, so that you end up with /tmp/Rhabarba.elev.ppm.

Before displaying any raster map, you can change a number of things about it. Among those is the color scheme to be used:

r.colors map=Rhabarba rules=elevation

or

r.colors map=Rhabarba color=aspect

The color tables can be one of the existing, but you can build also your own table as the format is quite straightforward to produce.

Adding a Relief

To add a shaded relief to an existing raster map is really easy:

r.shaded.relief map=Rhabarba shadedmap=RhabarbaRelief
                altitude=30 azimuth=270 zmult=1 scale=1

The newly generated RhabarbaRelief raster map contains the landscape seen with a sun position given by the azimuth and here 30 degrees altitude. I also found a nice trick to make the result optically pleasing:

d.his h_map=Rhabarba i_map=RhabarbaRelief

2

It takes the hue values from Rhabarba and the intensity values from RhabarbaRelief. After that you may want to save this as well:

d.out.file format=ppm output=/tmp/Rhabarba.rel

Adding Contours

Different to reliefs, contours will be created as vector map:

r.contour input=Rhabarba output=RhabarbaContour step=100
d.vect map=RhabarbaContour color=0:0:0 lcolor=0:0:0
    fcolor=170:170:170 display=shape
    type=point,line,boundary,centroid,area
    icon=basic/circle size=5 layer=1 lsize=8
    xref=left yref=center llayer=1
d.out.file format=ppm output=/tmp/Rhabarba.cont

The r.contour command creates this vector map, here with step size 100 (meters, elevation). After that, the vector map is rendered into the current display. Not sure what all the other parameters do.

Again, at the end we save this to a file.

All Together Now!

One of the niceties of the interactive application is that you can specify your layers for a display. Similar to Gimp or PhotoShop, you can specify what goes on top, how much transparency it should have and what is underneath.

3

Unfortunately, the underlying command - when executed - is referring to secretedly maintained temporary files. Worse, the command cloaks itself from end users:

g.pnmcomp isn't meant for end users. It's an internal tool for use by a Tcl/Tk GUI.

Feel free to make a deep sigh here.

It takes some guessing and source code reading to figure out how to actually use it. And no, the similarity with netpbm command of the same name is only superficial. Here is my coffee cup reading attempt:

The invocation

g.pnmcomp in=/tmp/Rhabarba.rel.ppm,/tmp/Rhabarba.cont.ppm
          mask=/tmp/Rhabarba.rel.pgm,/tmp/Rhabarba.cont.pgm
          opacity=1.0,0.3
          width=640 height=480
          background=255:255:255
          output=/tmp/Rhabarba.ppm

will:

  • create a background with the provided colors first (here white),
  • will create an image of geometry 640x480,
  • will then take the picture in /tmp/Rhabarba.cont.ppm,
  • overlay it with a gray mask in /tmp/Rhabarba.cont.pgm and attach the opacity 0.3 to it,
  • will then stack underneath it the picture in /tmp/Rhabarba.rel.ppm,
  • overlay that with a gray mask from /tmp/Rhabarba.rel.pgm with opacity 1.0,
  • and write the end result to /tmp/Rhabarba.ppm.

The opacity does not necessarily mean what the word says. Here an opacity of 1.0 means that everything underneath is invisible. Opacity 0.0 means that the layer itself is invisible and everything underneath is visible.

You may also wonder how these gray masks can be generated. I did not find a way with grass (there is probably one), but the netPBM library can do that easily:

pbmmake -white 640 480 | pbmtopgm 1 1 > /tmp/Rhabarba.cont.pgm
pbmmake -white 640 480 | pbmtopgm 1 1 > /tmp/Rhabarba.rel.pgm

All that produces

4

No more tears, but do not underestimate the computational costs.

AttachmentSize
Rhabarba.elev-small.jpg13.25 KB
Rhabarba.rel-small.jpg17.41 KB
grass-snap.jpg26.8 KB
Rhabarba.small.jpg19.83 KB
Posted In

Does OpenGL work with

Does OpenGL work with similar commands? I always wonder how one generates the grass and relief using that program.

appliance parts (not verified) | Mon, 06/23/2008 - 18:56