Submitted by dylan on Fri, 2006-03-17 20:13.
Initial classification of terrain based on 9 topographic parameters extracted from a 10 meter DEM. Slope, profile curvature, maximum curvature, minimum curvature, longitudinal curvature, and cross-sectional curvature were computed from a 15x15 moving window by r.param.scale. A 9x9 moving window was used to compute the hypsometric integral as defined by Wood and Snell (1960). These 7 parameters were extracted from a given region, and classified with the clara() function from the cluster package in R. Transferal of raster data between GRASS and R was made possible by the spgrass6 package. An excellent overview of some of the features in the spgrass6 package is included in the third volume of the GRASS newsletter. Some initial results are displayed in Figures 1 and 2. Terrain data is better visualized in 3 dimensions: Figure 3 contains a view of both Temblor and Volcanic lanforms near McCabe Canyon. Ideas for this type of classification technique came from Roger Bivand (Bivand, 2000).
Initial Calculation of Terrain Parameters:
#compute hypsometric integral, approixmiated by the elevation-relief-ratio
r.neighbors in=elev_meters out=avg size=$er_wind_size method=average
r.neighbors in=elev_meters out=min size=$er_wind_size method=minimum
r.neighbors in=elev_meters out=max size=$er_wind_size method=maximum
r.mapcalc "er = (avg - min)/(max - min) "
r.param.scale in=elev_meters size=$wind_size param=slope out=slope
r.param.scale in=elev_meters size=$wind_size param=profc out=profc
r.param.scale in=elev_meters size=$wind_size param=crosc out=crosc
r.param.scale in=elev_meters size=$wind_size param=minic out=minic
r.param.scale in=elev_meters size=$wind_size param=maxic out=maxic
r.param.scale in=elev_meters size=$wind_size param=longc out=longc
Clustering in R Note that these are just examples, and not neccessarily the best approach. Detailed knowledge of your study area, the terrain in question, as well as an understanding of local landscape dynamics are all important factors which will determine the best approach.
#load required packages
#load GRASS region data
gmeta6 <- gmeta6()
#read in our 7 raster files from GRASS
x <- readFLOAT6sp(c("er","crosc","longc","slope","profc","minic","maxic"))
#assemble a matrix of our terrain variables
morph <- data.frame(cbind(x$er, x$crosc, x$longc, x$slope, x$profc, x$minic, x$maxic))
# normailize slope by dividing my max(slope) -- seems to produce interesting results
# not neccessarily the best approach
morph <- data.frame(cbind(x$er, x$crosc, x$longc, x$slope/max(x$slope), x$profc, x$minic, x$maxic))
names(morph) <- c("er","crosc","longc","slope_n","profc","minic","maxic")
# perform the clustering
# stand=T <---- data is to be standardized before clustering
# probably a good idea as each map has different units, and scales
morph.clara <- clara(morph, k=5, stand=T)
#send result back to GRASS:
x$morph_class <- morph.clara$clustering
- Bivand, R. Integrating GRASS 5.0 and R: GIS and modern statistics Computers & Geosciences, 2000, 26, 1043–1052.
- Blaszczynski, J. Landform characterization with geographical information systems Photogrammetric Engineering and Remote Sensing, 1997, 63, 183-191
- Wood, W.F. and Snell, J.B. A Quatitative system for classifying landforms (Technical Report ER-124). U.S. Quatermaster Research & Engineering Center, Natick, MS. 1960.