Analyzing acacia tree health in the Arava with GRASS GIS

Using the GRASS GIS image analysis modules and color aerial photos, it is possible to create sets of landscape ecology statistics regarding the spread of acacia trees. A  comparison of these statistics for different years will indicate how stable and sustainable the acacia growth is along a particular wadi.

The workflow:

  1. Import a tiff format aerial ortho-photo into GRASS
  2. Using the ortho-photo image as a background, digitize a few training areas
  3. Run the SMAP classifier to create a new landcover map
  4. Isolate the “trees” class into a new map, and run a landscape ecology analysis on that map

Detailed instructions:

Note: commands to be typed in at the terminal are in blue mono typeface

Step 1

Import the original tiff file

GRASS 6.3.0 (ITM):~ > ortho/hatzeva_4_2007.tif out=hatzeva -o

The import creates three raster layers for each of the three color bands, red, green, blue.

GRASS 6.3.0 (ITM):~ > g.mlist rast

We set the region based on one of the rasters.


Now let’s display a composite of the rasters (first start a monitor):

d.mon x0
using default visual which is TrueColor
ncolors: 16777216
Graphics driver [x0] started

Now we choose the analysis area, and reset the region to cover just that area. Run d.where and by clicking in the monitor, choose the north, south, east, west extents to cover an area of one square kilometer.

g.region -p w=223700 e=224700 s=517200 n=518200
projection: 99 (Transverse Mercator)
zone:       0
datum:      towgs84=-48,55,52
ellipsoid:  grs80
north:      518200
south:      517200
west:       223700
east:       224700
nsres:      0.5
ewres:      0.5
rows:       2000
cols:       2000
cells:      4000000

We’ll create a permanent composite raster, clipped to the analysis region, for use in digitizing:

 r.composite out=hatzeva

The new raster will be clipped to the small region. Just to check, let’s display it:

d.rast hatzeva

Step 2

Next we digitize some sample training areas for trees, wadi, and savanna. Run the digitizer starting a new vector layer and with the composite raster as the background.

v.digit -n train bgcmd="d.rast hatzeva"

Some important tips when working with the digitizer.

  1. First, in the v.digit control panel window click the “Settings” button and go to the “Table” tab. Here enter the following columns for the training vector layer: “label” as type character, “label_id” as type integer, and grassrgb as type character.

    v.digit setting windows
    v.digit setting windows
  2. Because of the way the GRASS vector model is designed, each area feature is composed of both a boundary line and a centroid. Both must be digitized. First click the “Zoom in” button to zoom to a specific tree (note the use of the left and center mouse buttons for zooming). Then click the “Digitize new boundary” button.
  3. The area attributes are attached to the area centroid, and not the boundary. Notice in the v.digit control panel window there is a drop down labeled “Mode”. When digitizing boundaries, this will be set to “No category”, and when digitizing the centroids, reset it to “Next not used”. Now  begin marking the edges of a tree with the mouse. Click the right mouse button to finish the feature. screenshot-vdigit
  4. The boundary must be closed. The line will become green if the boundary is closed properly. If the end vertices stay red, then the boundary is not closed. Use the zoom in button to zoom very close to the end points, then use the “Move vertex” button to shift one exactly on top of the other. Now the line should become green.

    A closed boundary
    A closed boundary
  5. Now choose the “Digitize new centroid” button, change the mode to “Next not used” and click to mark a centroid inside the tree boundary. Now click the “Display attributes” button and click on the new centroid. In the attributes window enter “tree” as the label, “1” as the label_id and “050:160:085” as the grassrgb value.screenshot-form
  6. Continue in a similar fashion to digitize a few other trees, wadi areas and savanna areas. For wadis, enter the following values as centroid attributes: label=”wadi”, label_id=”2″, grassrgb=”255:255:210″. And for savanna areas enter the centroid attribute values: label=”savanna”, label_id=3, grassrgb=”250:150:040″.

Now we can view the digitized training areas, overlayed on the aerial photo:

g.region rast=hatzeva
d.rast hatzeva
d.vect -a train rgb_col=grassrgb

Finish this step by converting the vector to a raster map (by the same name): in=train out=train use=attr type=area col=label_id rgbcolumn=grassrgb labelcolumn=label

Step 3

Before running the classifier, a few preparatory steps are required. First define the group of raster images that will be used in the classification: group=hatzeva subgroup=hatzeva,,
group [hatzeva] - does not yet exist. Creating...
Adding files to group [hatzeva]
Adding raster map []
Adding raster map []
Adding raster map []
Writing group REF

Next create a signature file using our training map.

i.gensigset trainingmap=train group=hatzeva subgroup=hatzeva signaturefile=hatzeva.sig

The process will run for a few moments. After it finishes we start the actual classification:

i.smap group=hatzeva subgroup=hatzeva signaturefile=hatzeva.sig out=hatzeva_class

This runs for a while, produing many rows of output as it progresses, and finally creates a classified raster called hatzeva_class with three landcover classifications. We set reasonable colors for the classes:

r.colors hatzeva_class col=rules
Enter rules, "end" when done, "help" if you need it.
Data range is 1 to 3
> 1 50:160:85
> 2 255:255:110
> 3 250:150:40
> end
Color table for <hatzeva_class> set to rules

and display the landcover map to check results:

g.region -p rast=hatzeva_class
d.rast hatzeva_class

Step 4

At this point we can already see some basic statistical information about the classification:

GRASS 6.3.0 (ITM):~ > r.stats -a -p hatzeva_class
1 16787.500000   1.68%
2 606434.750000  60.64%
3 376777.750000  37.68%

Here we see the total area and percent area for each category where: 1=trees, 2=wadis, 3=savanna.

But we want to extract specific statistics attuned to landscape ecology analysis. The GRASS suite r.le.* offers several functions along this line. We begin by extracting a map of only the tree patches.

r.mapcalc trees="if( hatzeva_class==1, 1, null() )

This function creates a new raster with pixel value ‘1’ where the hatzeva classification value = 1, and NULL everywhere else.

Now we run the r.le.setup module to create a configuration file:

r.le.setup map=trees

This command opens the setup window where we choose the required options. The module allows defining a variety of analysis windows. We will simply use the (default) whole map as the analysis region.  First press ‘2’ “Setup a sampling frame”. At the first question “Will the sampling frame use the whole map” press ‘y’. Press ‘y’ a second time to refresh, and then press ‘7’ to exit the setup.

Now we run the r.le.patch module. We will choose, as an example, to extract the standard deviation of patch size.

GRASS 6.3.0 (ITM):~ > r.le.patch map=trees sam=w siz=s1,s2

    MAP:      trees
    SAMPLE:      whole map        TRACING:  8 neighbor
          mean patch size
          st. dev. patch size


Tracing patch    1327

The resulting calculation in placed into a directory “r.le.out” as a file named “s1-2.out” . The contents of the file look like:

GRASS 6.3.0 (ITM):~ > cat r.le.out/s1-2.out
Scale  Unit  MN. PATCH SIZE   S.D. PATCH SIZE -- in pixels
    1     1           50.603           91.082

We see that the mean patch (tree) size is 50 sq.m. and the standard deviation is 91.  If a full report of statistics for each indiviual patch is required, we add to the r.le.patch command the “out=…” option to create detailed statistics for each tree:

GRASS 6.3.0 (ITM):~ > r.le.patch map=trees sam=w siz=s1,s2 out=trees.txt

Now, in the r.le.out directory we find an additional text file “trees.txt” with a row for each patch including the patch size, perimeter, perimeter/area ratio and another landscape ecology measurement, the core perimeter/area ratio. This file can be imported into a spreadsheet program for further calculations and graphing.

ale  it    num   att    row  col     size     size     size      per    P/A   CP/A    RCC   number index
  1   1      1     1.0    6  202      112        0        0       50  0.446  1.332    inf        0 0.000
  1   1      2     1.0    8 1203       73        0        0      101  1.384  3.334    inf        0 0.000
  1   1      3     1.0    5 1476       95        0        0       48  0.505  1.389    inf        0 0.000
  1   1      4     1.0    4  608       12        0        0       18  1.500  1.465    inf        0 0.000

Leave a Reply

Your email address will not be published. Required fields are marked *