Home > GRASS, ushcn > The GRASS is Greener

The GRASS is Greener

2010 October 26

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.

USHCN Canon City

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.

g.list rast
## 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
d.mon x0
# display the two classes
d.rast.leg 51294_18b_class
# print a copy of the display
d.out.file output=51294_18b_class

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
d.erase -f
d.rast.leg 51294_18b_mode
d.out.file output=51294_18b_mode

51924_18 smoothed

We can read the percent area for each classification.

r.stats -p 51294_18b_mode
##1 77.32%
##2 22.68%

And create mask for one classification or the other

# mask the “urban” definition
d.erase -f
d.rgb b=51294_18.blue g=51294_18.green r=51294_18.red
d.rast 51294_18b_mode cat=2 -o
d.out.file output=51294_18b_mode_masked_cat2

51924_18 masked

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.

Advertisements
  1. 2010 October 27 at 3:59 pm

    Nice! But how would it react if there were trees with no foliage, or if the photo was taken in autumn in the Northeast?

    Also, what is the size of the images in question (in square meters)? And I don’t suppose the lat/lon coordinates of the instrument are good enough to calculate the distance from the sensor to the nearest artificial cover with a size greater than x?

    Once you feel comfortable with the tool, send an artificial cover % for each station area my way and I’ll play around with it a bit.

  2. 2010 October 27 at 4:45 pm

    Give me a few laps around the track and I will run a small set of the image classifications against the NOAA Impervious and see how well they track.

    Google is mercator, varies with latitude,

    This 640×640 square is (very) approximately 320 m x 320m ~ 100,000m^2 ~ 0.1km^2

  3. Steven Mosher
    2010 November 9 at 11:29 am

    Nice,

    A while back I wanted to give grass a go because of the links to R.

    You can transform the google map out of mercator if you want to. Just use raster “project” I’ve got the code to do that if you want it.

  1. 2010 October 27 at 7:09 pm
Comments are closed.