[GRASS-SVN] r55621 - grass-addons/grass7/raster/r.mess
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Apr 4 04:04:26 PDT 2013
Author: pvanbosgeo
Date: 2013-04-04 04:04:26 -0700 (Thu, 04 Apr 2013)
New Revision: 55621
Added:
grass-addons/grass7/raster/r.mess/r.mess
Log:
Copy from grass 6 addons
Copied: grass-addons/grass7/raster/r.mess/r.mess (from rev 55620, grass-addons/grass6/raster/r.mess/r.mess)
===================================================================
--- grass-addons/grass7/raster/r.mess/r.mess (rev 0)
+++ grass-addons/grass7/raster/r.mess/r.mess 2013-04-04 11:04:26 UTC (rev 55621)
@@ -0,0 +1,475 @@
+#!/bin/sh
+#
+#set -x
+########################################################################
+#
+# MODULE: r.mess
+# AUTHOR(S): Paulo van Breugel <p.vanbreugel AT gmail.com>
+# PURPOSE: Calculate the multivariate environmental similarity
+# surface (MESS) as proposed by Elith et al., 2010,
+# Methods in Ecology & Evolution, 1(330–342).
+#
+# NOTES: This GRASS script should give results very close (but
+# due to rounding differences probably not exactly the
+# same) to those calculated in Maxent. If you want to
+# compare, make sure to use the same points as used in
+# Maxent to calculate the MESS, which are, by default,
+# all background points if I understand it correctly. See
+# [fill in url] for a discussion about the assumptions
+# made and the use of the input parameter 'digits'
+#
+# Disclaimer: I use it in my work, but I am aware that it needs
+# improvements. Suggestions for improvements are most
+# welcome. In the meantime, use it at your own risk
+#
+# COPYRIGHT: (C) 2012 Paulo van Breugel
+# http://ecodiv.org
+# http://pvanb.wordpress.com/
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+########################################################################
+#
+#%Module
+#% description: Computes multivariate environmental similarity surface
+#%End
+
+#%option
+#% key: ref_rast
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Reference distribution as raster
+#% key_desc: name
+#% required: no
+#% multiple: no
+#% guisection: reference distribution
+#%end
+
+#%option
+#% key: ref_vect
+#% type: string
+#% gisprompt: old,vector
+#% description: Reference distribution as point vector layer
+#% key_desc: name
+#% required: no
+#% multiple: no
+#% guisection: reference distribution
+#%end
+
+#%option
+#% key: env_var
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input (predictor) raster map(s)
+#% key_desc: names
+#% required: yes
+#% multiple: yes
+#%end
+
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Root name of the output MESS data layers
+#% key_desc: name
+#% required: yes
+#%end
+
+#%option
+#% key: digits
+#% type: string
+#% description: Precision of your input layers values
+#% key_desc: string
+#% answer: 0.001
+#% required: yes
+#%end
+
+##%flag # not implemented yet
+##% key: n
+##% description: Keep individual environmental similarity layers (IES)
+##% guisection: Additional statistics
+##%end
+
+#%flag
+#% key: m
+#% description: Calculate Most dissimilar variable (MoD)
+#%end
+
+#%flag
+#% key: k
+#% description: Calculate mean of IES layers
+#%end
+
+#%flag
+#% key: l
+#% description: Calculate median of IES layers
+#%end
+
+## Set easier variable names
+OUTMAPS="${GIS_OPT_OUTPUT}"
+INMAPS="${GIS_OPT_ENV_VAR}"
+
+#=======================================================================
+## GRASS team recommandations
+#=======================================================================
+
+## Check if in GRASS
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." 1>&2
+ exit 1
+fi
+
+## check for awk
+if [ ! -x "$(which awk)" ] ; then
+ g.message -e "<awk> required, please install <awk> or <gawk> first"
+ exit 1
+fi
+
+## check for R
+if [ ! -x "$(which R)" ] ; then
+ g.message -e "<R> required, please install <R> first"
+ exit 1
+fi
+
+## To parse the code into interactive menu
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+## set environment so that awk works properly in all languages ##
+unset LC_ALL
+export LC_NUMERIC=C
+
+
+## what to do in case of user break:
+exitprocedure()
+{
+ echo "User break!"
+ cleanup
+ exit 1
+}
+
+## shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 15
+
+#=======================================================================
+## Config and general procedures
+#=======================================================================
+
+##fix this path
+if [ -z "$PROCESSDIR" ] ; then
+ PROCESSDIR="$HOME"
+fi
+
+#fix this path
+if [ -z "$LOGDIR" ] ; then
+ LOGDIR="$HOME"
+fi
+LOGFILE="$LOGDIR/r.mess.log"
+
+echo "r.mess :" >> "$LOGFILE"
+
+#=======================================================================
+## Preparing input variables for use in R script
+#=======================================================================
+
+# Create vectors with names output maps [arrOUT]
+# Create vector with names temporary layers [tmpLY]
+# Create vector with names input layers without mapset name [arrIN]
+# The 'removeThis' elements are to initiate the vector. How can I avoid this?
+
+arrOUT="${OUTMAPS}_MESS"
+tmpLY="removeThis"
+arrIN="removeThis2"
+counter=0
+IFS=,
+for nvar in ${INMAPS} ; do
+ arrOUT="$arrOUT;${OUTMAPS}_`echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`"
+ counter=`expr $counter + 1`
+ tmpLY="${tmpLY};tmp_mess_$$_$counter"
+ arrIN="${arrIN};`echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`"
+done
+
+if [ -z "$GIS_OPT_REF_VECT" -a -z "$GIS_OPT_REF_RAST" ];
+then
+ g.message -e message='Missing input: map with the reference distribution (point vector layer -or- a raster layer)'
+ exit 1
+else
+ if [ -n "$GIS_OPT_REF_VECT" -a -n "$GIS_OPT_REF_RAST" ]; then
+ g.message message='You gave a vector and raster layer as reference distribution layer. !Using the raster layer!'
+ REF_LAY="${GIS_OPT_REF_RAST}"
+ REF_TYPE="R"
+ else
+ if [ -n "$GIS_OPT_REF_RAST" ]; then
+ REF_LAY="${GIS_OPT_REF_RAST}"
+ REF_TYPE="R"
+ else
+ REF_LAY="${GIS_OPT_REF_VECT}"
+ REF_TYPE="V"
+ fi
+ fi
+fi
+
+#=======================================================================
+## testing if exising output maps already exist
+## and for missing input maps
+#=======================================================================
+
+# test for output raster map [1]
+g.findfile element=cell file=${OUTMAPS}_MESS > /dev/null
+ if [ $? -eq 0 ] ; then
+ g.message -e 'The output map '${OUTMAPS}'_MESS already exists'
+ exit 1
+fi
+
+# test for output raster maps [2]
+oIFS=$IFS
+IFS=,
+for nvar in $INMAPS ; do
+ arrIN=${OUTMAPS}_`echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`
+ g.findfile element=cell file=${arrIN} > /dev/null
+ if [ $? -eq 0 ] ; then
+ g.message -e 'The output map '${arrIN}' already exists'
+ exit 1
+ fi
+done
+IFS=$oIFS
+unset arrIN
+
+# test for missing input raster maps
+oIFS=$IFS
+IFS=,
+for nvar in $INMAPS ; do
+ tstIN=`echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`
+ g.findfile element=cell file=${tstIN} > /dev/null
+ if [ $? -gt 0 ] ; then
+ g.message -e 'The map '${tstIN}' is missing'
+ exit 1
+ fi
+done
+IFS=$oIFS
+unset tstIN
+
+#=======================================================================
+## Creating the R script
+#=======================================================================
+
+writeScript(){
+cat > $1 << "EOF"
+
+options(echo = FALSE)
+require(spgrass6)
+
+# Check grass version
+grassversion <- as.numeric(system("g.version | cut -c7", intern=T))
+
+## Get vector with variables
+argus <- commandArgs(trailingOnly=TRUE)
+ipn <- unlist(strsplit(argus[1],";"))[-1] # variable names
+ipl <- unlist(strsplit(argus[2],",")) # environmental layers
+opl <- unlist(strsplit(argus[3],";")) # output layers
+opi <- opl[-1] # base name individual layers
+opc <- opl[1] # name of MESS layer
+tml <- unlist(strsplit(argus[4], ";"))[-1] # temporary layers
+vtl <- argus[5] # reference distribution raster layer
+rtl <- argus[10] # raster or vector layer
+digits <- as.numeric(argus[9]) # Precision
+rdigits <- nchar((1/digits)-1)
+
+#-----------------------------------------------------------------------
+# Create the r.mapcalc expressions to calculate the mess for the
+# individual layers. The main step is to define the graph function to
+# be used -
+#-----------------------------------------------------------------------
+
+if(rtl=="R"){
+
+ # Reference distribution layer is raster
+ #-----------------------------------------------------------------------
+
+ flnm1 <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".txt")
+ flnm2 <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".txt")
+ for(i in 1:length(ipl)){
+
+ # Import the environmental data ipl[i] in R. I could use r.stats to
+ # get frequency table rather then raw data. But I first need to find
+ # out how to deal with floating data in r.stats, as counts of each
+ # unique value is needed.
+ # One option is to convert the floating layers to integer layers after
+ # multiplying first with a number equal to the reciprocal of the value
+ # given under the digits option
+
+ fft <- system("g.findfile --quiet element=cell file=MASK", ignore.stdout = TRUE)
+ if(fft==0){
+ rname <- paste("MASK",paste(sample(c(0:9, letters, LETTERS), 12, replace=TRUE), collapse=""), sep="_")
+ execGRASS("g.copy", rast=paste("MASK", rname, sep=","))
+ system("g.remove -f rast=MASK")
+ }
+
+ # make compatible for both v6.4 and 7
+ if(grassversion==7){
+ system(paste("r.mask raster=", vtl, sep=""))
+ }else{
+ system(paste("r.mask input=", vtl, sep=""))
+ }
+
+ system(paste("r.out.xyz input=", ipl[i], " output=", flnm1, " --overwrite", sep=""))
+ system(paste("awk -F \"|\" '{print $3}' ", flnm1, " > ", flnm2, sep=""))
+ spld <- scan(file=flnm2)
+ #unlink(c(flnm1, flnm2))
+ spld <- round(spld, rdigits)
+
+ # Calculate the frequency distribution
+ a <- table(spld)
+ envmin <- min(spld)
+ envmax <- max(spld)
+ Drange <- system(paste("r.info -r --verbose map=", ipl[i], sep=""), intern=T)
+ Dmin <- as.numeric(sub("min=", "", Drange[1]))
+ Dmax <- as.numeric(sub("max=", "", Drange[2]))
+ system("g.remove -f rast=MASK")
+
+ # Create recode rules
+ x1 <- c(Dmin, as.numeric(rownames(a))+digits)
+ x2 <- c(envmin, as.numeric(rownames(a))[-1],Dmax)
+ y1 <- c(0,cumsum(as.numeric(a))/sum(a)*100)
+ xy1 <- format(cbind(x1,x2,y1), scientific=F, trim=T, nsmall=rdigits)
+ xy2 <- apply(xy1, 1, function(x){paste(x, collapse=":")})
+ tmp.rule <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".rules")
+ write.table(xy2, file=tmp.rule, quote=F, row.names=F, col.names=F)
+
+ # Restore mask
+ if(fft==0){
+ execGRASS("g.rename", rast=paste(rname,"MASK", sep=","), intern=TRUE)
+ }
+
+ # Create the recode layer and calculate the IES
+ system(paste("r.recode input=", ipl[i], " output=", tml[i], " rules=", tmp.rule, sep=""))
+ system(paste("r.mapcalc '", opi[i], " = if(", tml[i], "==0, (float(", ipl[i], ")-", envmin, ")/(", envmax, "-", envmin, ") *100.0, if(", tml[i], "<=50, 2 * ", tml[i], ", if(", tml[i], "<100.0, 2 * (100.0- ", tml[i], "), (", envmax, "- float(", ipl[i], "))/(", envmax, "-", envmin, ") * 100.0)))'", sep=""))
+ system(paste("g.remove rast=", tml[i], sep=""))
+ unlink(tmp.rule)
+ }
+
+}else{
+
+ # Reference distribution layer is vector
+ #-----------------------------------------------------------------------
+
+ system(paste("v.extract -t input=", vtl, " type=point output=tmpMESS976543210", sep=""))
+ system(paste("v.db.addtable tmpMESS976543210 columns='", paste(ipn, " double precision", collapse=","), "'", sep=""))
+
+ # make compatible for both v6.4 and 7
+ if(grassversion==7){
+ for(m in 1:length(ipn)){
+ system(paste("v.what.rast map=tmpMESS976543210 layer=1 raster=", ipl[m], " column=", ipn[m], sep=""))
+ }
+ }else{
+ for(m in 1:length(ipn)){
+ system(paste("v.what.rast vector=tmpMESS976543210 layer=1 raster=", ipl[m], " column=", ipn[m], sep=""))
+ }
+ }
+
+ spld <- paste(ipn, collapse=",")
+ b <- execGRASS("v.db.select", parameters=list(map="tmpMESS976543210", columns=spld), intern=TRUE)
+ con <- textConnection(b)
+ spl <- na.omit(read.table(con, header=TRUE, sep="|"))
+
+ # Check for point without values
+ clpf <- na.action(spl)
+ if(length(clpf)==1){print(paste("Please note that the point", clpf, "has no value. This is probably because it lies outside the current region"))}
+ if(length(clpf)>1){print(paste("Please note that the points", paste(clpf, collapse=" and "), "have no values. This is probably because they lie outside the current region"))}
+
+ close(con)
+ system("g.remove vect=tmpMESS976543210")
+
+ for(i in 1:length(ipl)){
+ # Calculate the frequency distribution
+ a <- table(spl[i])
+ envmin <- min(spl[i])
+ envmax <- max(spl[i])
+ Drange <- system(paste("r.info -r --verbose map=", ipl[i], sep=""), intern=T)
+ Dmin <- as.numeric(sub("min=", "", Drange[1]))
+ Dmax <- as.numeric(sub("max=", "", Drange[2]))
+
+ # Create recode rules
+ x1 <- c(Dmin, as.numeric(rownames(a))+digits)
+ x2 <- c(envmin, as.numeric(rownames(a))[-1],Dmax)
+ y1 <- c(0,cumsum(as.numeric(a))/sum(a)*100)
+ xy1 <- format(cbind(x1,x2,y1), scientific=F, trim=T, nsmall=rdigits)
+ xy2 <- apply(xy1, 1, function(x){paste(x, collapse=":")})
+ tmp.rule <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = ".rules")
+ write.table(xy2, file=tmp.rule, quote=F, row.names=F, col.names=F)
+
+ # Create the recode layer and calculate the IES
+ system(paste("r.recode input=", ipl[i], " output=", tml[i], " rules=", tmp.rule, sep=""))
+ system(paste("r.mapcalc '", opi[i], " = if(", tml[i], "==0, (float(", ipl[i], ")-", envmin, ")/(", envmax, "-", envmin, ") *100.0, if(", tml[i], "<=50, 2 * ", tml[i], ", if(", tml[i], "<100.0, 2 * (100.0- ", tml[i], "), (", envmax, "- float(", ipl[i], "))/(", envmax, "-", envmin, ") * 100.0)))'", sep=""))
+ system(paste("g.remove rast=", tml[i], sep=""))
+ unlink(tmp.rule)
+ }
+}
+
+#-----------------------------------------------------------------------
+# Calculate MESS
+#-----------------------------------------------------------------------
+
+system(paste("r.series output=", opc, " input=", paste(opi, collapse=","), " method=minimum", sep=""))
+
+#-----------------------------------------------------------------------
+# Calculate other stats
+#-----------------------------------------------------------------------
+
+if(argus[6]==1){
+ system(paste("r.series output=", opc, "_MoD input=", paste(opi, collapse=","), " method=min_raster", sep=""))
+ nuv <- cbind(seq(from=0, to=length(ipn)-1, by = 1), ipn)
+ reclvar <- apply(nuv,1,function(x){paste(x[1],x[2], sep=":")})
+ tmpclas <- tempfile(pattern = "classification.rules.")
+ sink(tmpclas)
+ for(i in 1:length(reclvar)){
+ cat(reclvar[i]); cat("\n")
+ }
+ sink()
+ system(paste("r.category map=", opc, "_MoD rules=", tmpclas, sep=""))
+ unlink(tmpclas)
+}
+if(argus[7]==1){
+ system(paste("r.series output=", opc, "_mean input=", paste(opi, collapse=","), " method=average", sep=""))
+}
+if(argus[8]==1){
+ system(paste("r.series output=", opc, "_median input=", paste(opi, collapse=","), " method=median", sep=""))
+}
+
+EOF
+}
+
+# RGrass script generation
+# --------------------------
+RGRASSSCRIPT="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$RGRASSSCRIPT" ] ; then
+ g.message -e 'ERROR: unable to create temporary file for RGrass script' 1>&2
+ exit 1
+fi
+writeScript "$RGRASSSCRIPT"
+
+
+#=======================================================================
+## RGrass call
+#=======================================================================
+
+# Message
+g.message message='Calculating MESS layers.. this may take some time'
+
+##using R to create MESS layers
+#--slave
+R --vanilla --args "$arrIN" "${INMAPS}" "$arrOUT" "$tmpLY" "${REF_LAY}" "$GIS_FLAG_M" "$GIS_FLAG_K" "$GIS_FLAG_L" "${GIS_OPT_DIGITS}" "${REF_TYPE}" < "$RGRASSSCRIPT" >> "$LOGFILE" 2>&1
+if [ $? -ne 0 ] ; then
+ echo "ERROR: an error occurred during R script execution" 1>&2
+ exit 1
+fi
+
+g.message "Finished - you can find the maps in your current mapset"
+
+#=======================================================================
+
+
+
More information about the grass-commit
mailing list