The GRASS is Greener
GRASS is the common name for Geographic Resources Analysis Support System. It has a huge GIS (Geographic Information Systems) toolkit and is under constant development. I’ve been meaning to take it out on a spin for some cluster analysis and surface classification. Better late than never.
A couple of months ago, I used RGoogleMaps to automate the download of some USHCN sites from Google Maps. We can use GRASS tools to classify the various pixel colors into an arbitrary set of classes using cluster analysis. I’ll use the same image, Canon City CO, that I used in the earlier post.
I am a two-day newbie to GRASS. The following steps were derived from instructions at the GRASS Image Classification wiki and the PerryGeo: Impervious surface deliniation with GRASS blog post.
Download and install GRASS. Download the above image into the GRASS GIS project directory. Open a command terminal, cd into the same project directory, and run the command:
r.in.gdal -e input=51294_18.png output=51294_18 location=51294_18
Start GRASS (command ‘grass’) and select the ‘51294_18’ project location and the ‘PERMANENT’ mapset. Select the ‘Enter GRASS’ button. Exit the GRASS gui. Your command line prompt should now read something like the following with the 51924_18 project displayed in the command prompt: GRASS 6.4.0RC5+39438 (51294_18):~/projects/grassdata
You can check that the image has been imported and the project initialized. The ‘##’ denotes output.
## raster files available in mapset :
## 51294_18.alpha 51294_18.blue 51294_18.green 51294_18.red
Group the raster color files from the image back together
i.group group=51294_18 subgroup=51294_18i input=51294_18.blue@PERMANENT,51294_18.green@PERMANENT,51294_18.red@PERMANENT
i.target -c group=51294_18
i.group -l -g group=51294_18 subgroup=51294_18
I has a slight hickup in the above, and when I reviewed the file 51294_18/PERMANENT/group/51294_18/REF, the r-g-b mappings were off. I corrected them by hand. Probably a syntax or ordering issue in the above command.
Now do the cluster analysis. We divide the image into two clusters (aka classes). This is an unsupervised classification.
i.cluster group=51294_18 subgroup=51294_18 classes=2 sigfile=51294_18_sig.txt reportfile=51294_18_rpt.txt
i.maxlik group=51294_18 subgroup=51294_18 sigfile=51294_18_sig.txt class=51294_18b_class reject=51294_18b_reject
# open a display window
# display the two classes
# print a copy of the display
We can smooth the image out by taking the most common pixel in a 3×3 group. This helps filter out shadows and other piecemeal noise.
r.neighbors input=51294_18b_class output=51294_18b_mode method=mode size=3
We can read the percent area for each classification.
r.stats -p 51294_18b_mode
And create mask for one classification or the other
# mask the “urban” definition
d.rgb b=51294_18.blue g=51294_18.green r=51294_18.red
d.rast 51294_18b_mode cat=2 -o
Note that not everything that is marked in red is a man-made structure of some sort. It just matches a r-g-b cluster which happens to be mostly man-made for classification 2 (red). But areas that are otherwise disturbed or denuded of vegetation also show up within that classification – which happens to be pretty convenient.
Gonna have to do something about the Google trademarks and taglines.