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.
- Import a tiff format aerial ortho-photo into GRASS
- Using the ortho-photo image as a background, digitize a few training areas
- Run the SMAP classifier to create a new landcover map
- Isolate the “trees” class into a new map, and run a landscape ecology analysis on that map
Note: commands to be typed in at the terminal are in blue mono typeface
Import the original tiff file
GRASS 6.3.0 (ITM):~ >r.in.gdal 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 d.rgb red=hatzeva.red green=hatzeva.green blue=hatzeva.blue
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)
We’ll create a permanent composite raster, clipped to the analysis region, for use in digitizing:
r.composite red=hatzeva.red green=hatzeva.green blue=hatzeva.blue out=hatzeva
The new raster will be clipped to the small region. Just to check, let’s display it:
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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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:
d.vect -a train rgb_col=grassrgb
Finish this step by converting the vector to a raster map (by the same name):
v.to.rast in=train out=train use=attr type=area col=label_id rgbcolumn=grassrgb labelcolumn=label
Before running the classifier, a few preparatory steps are required. First define the group of raster images that will be used in the classification:
i.group group=hatzeva subgroup=hatzeva in=hatzeva.red,hatzeva.green,hatzeva.blue group [hatzeva] - does not yet exist. Creating... Adding files to group [hatzeva] Adding raster map [hatzeva.red] Adding raster map [hatzeva.green] Adding raster map [hatzeva.blue] Writing group REF Done.
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
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:
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 PARAMETER CHOICES: MAP: trees SAMPLE: whole map TRACING: 8 neighbor SIZE MEASURES: mean patch size st. dev. patch size R.LE.PATCH IS WORKING....; Tracing patch 1327 R.LE.PATCH IS DONE; OUTPUT FILES IN SUBDIRECTORY "r.le.out"
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