Recently I’ve been doing work on determining oriented landscapes on Eastern Victoria Island. If you’re familiar with the Canadian NTS index I’m working primarily in NTS 77E for testing. Surficial geologists are very interested in processes that assist in explaining genesis. The landscape from a regional perspective can be broken down into two broad categories: 1)oriented & 2) non-oriented. Oriented really means stream line like landforms (think drumlins, flutes, eskers) that have an azimuth & can be grouped into fields. All other terrain, primarily hummocky are grouped into non-oriented.
Using the artefact rich CDED dem, I have a work-flow that identifies larger oriented features. Now the trick is to figure out the orientation that each feature exhibits & further begin to group fields that share similar orientations. I have not mentioned direction because that implies an angle & “flow” direction. By finding the streamline forms, the surficial geologist can determining from the drumlin fields what the direction of ice travel was. So the orientations I calculate will need to be adjusted after consulting with them. Not a big issue really, once you have the algorithm down for this.
I have been using GRASS within R. The spgrass6 package is great as it allows the creation of temporary grass location & mapsets for use. For this to work you need to have GRASS installed.
## load libraries
library(rgdal)
library(spgrass6)
library(raster)
## create tmp grass location
initGRASS(gisBase=”/usr/local/grass-6.4.3svn/”,home=tempdir(),override=T)
# gmeta<-gmeta6()
## change to permanent mapset to apply projection change for region
execGRASS(“g.mapset”,flags=c(“c”),parameters=list(mapset=”PERMANENT”))
## convert default coordinate system x,y to epsg code 3347 (can be altered later)
execGRASS(“g.proj”,flags=c(“c”),epsg=as.integer(3347))
## import clusters majority. “o” flag uses location proj. location & raster epsg must be same.
execGRASS(“r.in.gdal”,flags=c(“overwrite”,”o”),parameters=list(input=”../Tmp/saga_walk_through/Clusters_gdal.sdat”,output=”clustersInput”))
## review env settings
execGRASS(“g.gisenv”,flags=”s”)
# ## list dems
execGRASS(“g.list”,parameters=list(type=”rast”))
## set region to import clusters majority raster
execGRASS(“g.region”,flags=c(“p”,”a”),parameters=list(rast=”clustersInput”))
## apply mode 3*3 to clusters
execGRASS(“r.neighbors”,flags=c(“overwrite”),parameters=list(input=”clustersInput”,output=”clustersMode3″,method=”mode”))
### write vect to R obj for viewing
## read org clusters mode 3*3 into R
rastMode3Clusters<-readRAST6(“clustersMode3″)
## convert the spatialgriddataframe to raster object
rMode3Clusters<-raster(rastMode3Clusters)
## plot
plot(rMode3Clusters)
## write to disk for future reference. INT1U == BYTE
writeRaster(rMode3Clusters,filename=”Tmp/orgClustersMode3.tif”,format=”GTiff”,datatype=”INT1U”,NAFlag=as.integer(255),overwrite=TRUE)
I hope this gives you a quick idea of using GRASS within R.
I will be posting more on this matter in the coming days & will get a bit more technical with it.
Cheers,
Richard
