[GRASS-SVN] r63969 - grass-addons/grass7/raster/r.mess
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jan 6 14:09:48 PST 2015
Author: pvanbosgeo
Date: 2015-01-06 14:09:48 -0800 (Tue, 06 Jan 2015)
New Revision: 63969
Added:
grass-addons/grass7/raster/r.mess/r.mess.py
Removed:
grass-addons/grass7/raster/r.mess/r.mess
Modified:
grass-addons/grass7/raster/r.mess/r.mess.html
Log:
Rewrite of the shell/R script to python.
Deleted: grass-addons/grass7/raster/r.mess/r.mess
===================================================================
--- grass-addons/grass7/raster/r.mess/r.mess 2015-01-06 22:03:42 UTC (rev 63968)
+++ grass-addons/grass7/raster/r.mess/r.mess 2015-01-06 22:09:48 UTC (rev 63969)
@@ -1,583 +0,0 @@
-#!/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) 2014 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_old
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Input (predictor) raster map(s)
-#% key_desc: names
-#% required: yes
-#% multiple: yes
-#% guisection: predictors
-#%end
-
-#%option
-#% key: env_new
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Input (predictor) raster map(s)
-#% key_desc: names
-#% required: no
-#% multiple: yes
-#% guisection: predictors
-#%end
-
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: Root name of the output MESS data layers
-#% key_desc: name
-#% required: yes
-#% guisection: Output
-#%end
-
-#%option
-#% key: digits
-#% type: integer
-#% description: Precision of your input layers values
-#% key_desc: string
-#% answer: 3
-#%end
-
-#%flag: IES
-#% key: i
-#% description: Keep individual environmental similarity layers (IES)
-#% guisection: Output
-#%end
-
-#%flag
-#% key: m
-#% description: Calculate Most dissimilar variable (MoD)
-#% guisection: Output
-#%end
-
-#%flag
-#% key: k
-#% description: Calculate mean of IES layers
-#% guisection: Output
-#%end
-
-#%flag
-#% key: l
-#% description: Calculate median of IES layers
-#% guisection: Output
-#%end
-
-#%flag
-#% key: n
-#% description: Area with negative MESS
-#% guisection: Output
-#%end
-
-
-#%option
-#% key: liblocs
-#% type: string
-#% description: Location R libraries
-#% key_desc: folder
-#% required: no
-#%end
-
-## Set easier variable names
-OUTMAPS="${GIS_OPT_OUTPUT}"
-export INMAPS1="${GIS_OPT_ENV_OLD}"
-export INMAPS2="${GIS_OPT_ENV_NEW}"
-
-#=======================================================================
-## 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 input layers without mapset name [arrIN1]
-# The 'removeThis' elements are to initiate the vector. How can I avoid this?
-
-arrOUT="${OUTMAPS}_MESS"
-arrIN="removeThis2"
-counter=0
-IFS=,
-for nvar in ${INMAPS1} ; do
- export arrOUT="$arrOUT;${OUTMAPS}_`echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`"
- counter=`expr $counter + 1`
- export 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!'
- export REF_LAY="${GIS_OPT_REF_RAST}"
- export REF_TYPE="R"
- else
- if [ -n "$GIS_OPT_REF_RAST" ]; then
- export REF_LAY="${GIS_OPT_REF_RAST}"
- export REF_TYPE="R"
- else
- export REF_LAY="${GIS_OPT_REF_VECT}"
- export 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 $INMAPS1 ; do
- arrTEST=${OUTMAPS}_`echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`
- g.findfile element=cell file=${arrTEST} > /dev/null
- if [ $? -eq 0 ] ; then
- g.message -e 'The output map '${arrTEST}' already exists'
- exit 1
- fi
-done
-IFS=$oIFS
-unset arrTEST
-
-# test for missing input raster maps
-oIFS=$IFS
-IFS=,
-for nvar in $INMAPS1 ; 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 = TRUE)
-
-# Install (if not present) and load required packages
-
-tmp.pack <- tempdir()
-libLocs <- Sys.getenv("GIS_OPT_LIBLOCS") # Location R packages
-if(libLocs!=""){
- libLocs <- c(.libPaths(), libLocs, tmp.pack)
-}else{
- libLocs <- c(.libPaths(), tmp.pack)
-}
-
-if (!"XML" %in% installed.packages(lib.loc=libLocs)){
- install.packages("XML", dep=TRUE, repos='http://cran.us.r-project.org', tmp.pack)
- require(XML, lib.loc=libLocs)
-}else{require(XML, lib.loc=libLocs)}
-
-if (!"sp" %in% installed.packages(lib.loc=libLocs)){
- install.packages("sp", dep=TRUE, repos='http://cran.us.r-project.org', tmp.pack)
- require(sp, lib.loc=libLocs)
-}else{require(sp, lib.loc=libLocs)}
-
-if (!"spgrass6" %in% installed.packages(lib.loc=libLocs)){
- install.packages("spgrass6", dep=TRUE, repos='http://cran.us.r-project.org', tmp.pack)
- require(spgrass6, lib.loc=libLocs)
-}else{require(spgrass6, lib.loc=libLocs)}
-
-outputGRASS <- function (x, separator = ",", h = FALSE)
-{
- con <- textConnection(x)
- MyVar <- read.table(con, header = h, sep = separator, comment.char = "")
- close(con)
- return(MyVar)
-}
-
-# Check grass version
-grassversion <- as.numeric(system("g.version | cut -c7", intern=T))
-
-## Get vector with variables
-ipn <- Sys.getenv("arrIN") # variable names
-ipn <- unlist(strsplit(ipn,";"))[-1]
-ipn <- tolower(ipn) # avoid problems with vector column names
-ipl <- Sys.getenv("INMAPS1") # old environmental layers
-ipl <- unlist(strsplit(ipl,","))
-ipl2 <- Sys.getenv("INMAPS2") # new environmental layers
-if(ipl2==""){
- ipl2 <- ipl
-}else{
- ipl2 <- unlist(strsplit(ipl2,","))
- if(length(ipl)!=length(ipl2)){
- stop("list with old and new environmental data layers is not of the same length")
- }
-}
-opl <- Sys.getenv("arrOUT") # output layers
-opl <- unlist(strsplit(opl,";"))
-opi <- opl[-1] # base name individual layers
-opc <- opl[1] # name of MESS layer
-vtl <- Sys.getenv("REF_LAY") # reference distribution raster layer
-rtl <- Sys.getenv("REF_TYPE") # raster or vector layer
-flm <- Sys.getenv("GIS_FLAG_M")
-flk <- Sys.getenv("GIS_FLAG_K")
-fll <- Sys.getenv("GIS_FLAG_L")
-fln <- Sys.getenv("GIS_FLAG_N")
-IL <- Sys.getenv("GIS_FLAG_I")
-digits <- as.numeric(Sys.getenv("GIS_OPT_DIGITS")) # Precision
-digits2 <- 10^digits
-options(echo=TRUE, digits=digits+4, scipen=digits+4)
-
-#-----------------------------------------------------------------------
-# Create the r.mapcalc expressions to calculate the mess for the
-# individual layers. The main step is to define the recode rules to
-# be used -
-#-----------------------------------------------------------------------
-
-if(rtl=="R"){
-
- # Reference distribution layer is raster
- #-----------------------------------------------------------------------
-
- for(i in 1:length(ipl)){
-
- # Get the environmental frequency distributions for the reference points in R.
-
- tmpf0 <- "rmess_tmp_fr9876543210a"
- execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmpf0, "=", vtl, "* 1"))
-
- citiam <- system("g.findfile --quiet element=cell file=MASK", ignore.stdout = TRUE)
- if(citiam==0){
- rname <- paste("MASK",paste(sample(c(0:9, letters, LETTERS), 12, replace=TRUE), collapse=""), sep="_")
- execGRASS("g.copy", raster=paste("MASK", rname, sep=","))
- system("g.remove -bf type=rast name=MASK")
- }
-
- # Create mask based on combined MASK/input layer
- # Make compatible for both v6.4 and 7
- if(grassversion==7){
- execGRASS("r.mask", raster=tmpf0)
- }else{
- execGRASS("r.mask", input=tmpf0)
- }
-
- # Calculate the frequency distribution
- tmpf1 <- "rmess_tmp_fr987654321a"
- execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmpf1, " = int(", digits2, "*", ipl[i], ")"))
- frtab <- execGRASS("r.stats", flags=c("c", "n"), input=tmpf1, separator=",", sort="asc", intern=TRUE)
- a <- outputGRASS(frtab, separator=",")
- a <- a[order(a[,1]),]
- b <- cumsum(as.numeric(a[,2]))/sum(a[,2])*100
-
- tmpf2 <- "rmess_tmp_fr987654321b"
- execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmpf2, " = int(", digits2, "*", ipl2[i], ")"))
- d <- execGRASS("r.univar", map=tmpf2, flags="g", intern=TRUE)
- d <- outputGRASS(d, separator="=")
-
- # Create recode rules
- minv <- as.numeric(subset(d, select=V2, V1=="min"))
- maxv <- as.numeric(subset(d, select=V2, V1=="max"))
- e1 <- ifelse(minv < min(a[,1]), minv, min(a[,1])-1)
- e2 <- ifelse(maxv > max(a[,1]), maxv, max(a[,1])+1)
- a1 <- c(e1, a[1,1], a[,1]+1)
- a2 <- c(a[1,1], a[,1], e2)
- b1 <- c(0, b, 101)
- b1 <- round(b1*digits2)
- c1 <- cbind(a1,a2, b1)
- xy2 <- apply(c1, 1, function(x){paste(format(x, scientific=FALSE, trim=TRUE), 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
- # ToDo: as option, set region to env_new (and do not restore MASK, or have an option to set
- # a new mask for the new scenario. This will make it possible to calculate MESS for new
- # geographic areas outside the current region
- if(citiam==0){
- execGRASS("r.mask", flags="overwrite", raster=rname)
- execGRASS("g.remove", flags="f", type="raster", name=rname)
- }else{
- execGRASS("r.mask", flags="r")
- }
-
- # Create the recode layer and calculate the IES
- execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmpf2, " = int(", digits2, "*", ipl2[i], ")"))
- tmpf3 <- "rmess_tmp_fr987654321c"
- execGRASS("r.recode", input=tmpf2, output=tmpf3, rules=tmp.rule, flags="overwrite")
- min1 <- min(a[,1])
- max1 <- max(a[,1])
- thr1 <- format(50*digits2, scientific=FALSE)
- thr2 <- format(100*digits2, scientific=FALSE)
- z1 <- paste("if(", tmpf3, "==0, ", digits2, "* 100.0 * (", tmpf2, "-", min1, ")/(", max1, "-", min1, ")", sep="")
- z2 <- paste(", if(", tmpf3, "<=", thr1, ", 2*", tmpf3, sep="")
- z3 <- paste(", if(", tmpf3, "<=", thr2, ", 2*(", thr2, "-", tmpf3, ")", sep="")
- z4 <- paste(",", digits2, " * 100.0 * (", max1, "- ", tmpf2, ")/(", max1, "-", min1, "))))", sep="")
- calcc <- paste(z1, z2, z3, z4, sep="")
- execGRASS("r.mapcalc", expression=paste(opi[i], " = ", calcc, sep=""), flags="overwrite")
- execGRASS("r.mapcalc", expression=paste(opi[i], " = ", opi[i], "/", digits2, ".0", sep=""), flags="overwrite")
- execGRASS("g.remove", flags="bf", type="raster", name=paste(tmpf0, tmpf1, tmpf2, tmpf3, sep=","))
- unlink(tmp.rule)
- }
-}else{
-
- # Reference distribution layer is vector
- #-----------------------------------------------------------------------
-
- execGRASS("v.extract", flags="t", input=vtl, type="point", output="tmpMESS976543210")
- system(paste("v.db.addtable tmpMESS976543210 columns='", paste(ipn, " double precision", collapse=","), "'", sep=""))
- execGRASS("db.execute", sql="CREATE UNIQUE INDEX tmpMESS123 ON tmpMESS976543210 (cat)")
-
- # 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 -f type=vect name=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=", ipl2[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)))
- 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=digits)
- 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
- tmpv1 <- "rmess_tmp_fr987654321v1"
- execGRASS("r.recode", flags="overwrite", input=ipl2[i], output=tmpv1, rules=tmp.rule)
- z1 <- paste("if(", tmpv1, "==0, (double(", ipl2[i], ")-", envmin, ".0)/(", envmax, ".0-", envmin, ".0) *100.0", sep="")
- z2 <- paste(", if(", tmpv1, "<=50, 2.0 * ", tmpv1)
- z3 <- paste(", if(", tmpv1, "<100, 2.0 * (100.0- ", tmpv1, ")")
- z4 <- paste(", (", envmax, ".0- double(", ipl2[i], "))/(", envmax, ".0-", envmin, ".0) * 100.0)))", sep="")
- calcc <- paste(z1, z2, z3, z4, sep="")
- execGRASS("r.mapcalc", expression=paste(opi[i], "=", calcc), flags="overwrite")
- execGRASS("g.remove", flags="f", type="raster", name=tmpv1)
- unlink(tmp.rule)
- }
-}
-#-----------------------------------------------------------------------
-# Calculate MESS
-#-----------------------------------------------------------------------
-
-system(paste("r.series output=", opc, " input=", paste(opi, collapse=","), " method=minimum", sep=""))
-
-#-----------------------------------------------------------------------
-# Calculate other stats
-#-----------------------------------------------------------------------
-
-if(fln==1){
- execGRASS("r.mapcalc", expression=paste(opc, "_neg = if(", opc, "<0,1)", sep=""))
-}
-
-if(flm==1){
- system(paste("r.series output=", opc, "_MoD input=", paste(opi, collapse=","), " method=min_raster", sep=""))
- nuvto <- length(ipn)-1
- nuv <- cbind(seq(from=0, to=nuvto, 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(flk==1){
- system(paste("r.series output=", opc, "_mean input=", paste(opi, collapse=","), " method=average", sep=""))
-}
-if(fll==1){
- system(paste("r.series output=", opc, "_median input=", paste(opi, collapse=","), " method=median", sep=""))
-}
-
-if(IL==0){
- system(paste("g.remove -f type=rast name=", paste(opi, collapse=","), 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
-# R --no-save --no-restore --no-site-file --no-init-file --args "$arrIN" "${INMAPS1}" "${INMAPS2}" "$arrOUT" "$tmpLY" "$REF_LAY" "$GIS_FLAG_M" "$GIS_FLAG_K" "$GIS_FLAG_L" "$GIS_OPT_DIGITS" "$REF_TYPE"
-R --no-save --no-restore < "$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"
-
-#=======================================================================
-
Modified: grass-addons/grass7/raster/r.mess/r.mess.html
===================================================================
--- grass-addons/grass7/raster/r.mess/r.mess.html 2015-01-06 22:03:42 UTC (rev 63968)
+++ grass-addons/grass7/raster/r.mess/r.mess.html 2015-01-06 22:09:48 UTC (rev 63969)
@@ -1,17 +1,25 @@
<h2>DESCRIPTION</h2>
-<em>r.mess</em> computes the Multivariate Environmental Similarity (MES) surfaces in GRASS using R as backend. The MES index was proposed by Elith et al (2010) and is also implemented in the Maxent software. The MES approach can be described as follows (from Elith et al 2010): "The multivariate environmental similarity surface (MESS) calculation represents how similar a point is to a reference set of points, with respect to a set of predictor variables (V1, V2, ...). The values in the MESS are influenced by the full distribution of the reference points, so that sites within the environmental range of the reference points but in relatively unusual environments will have a smaller value than those in very common environments."
+<p>The Multivariate Environmental Similarity (MES) surfaces was proposed by Elith et al (2010) and originally implemented in the Maxent software. They described the MES approach as: "The multivariate environmental similarity surface (MESS) calculation represents how similar a point is to a reference set of points, with respect to a set of predictor variables (V1, V2, ...). The values in the MESS are influenced by the full distribution of the reference points, so that sites within the environmental range of the reference points but in relatively unusual environments will have a smaller value than those in very common environments."
-<p>This module will compute the individual environmental similarity surfaces (IES), which represents how similar each raster cell in a raster layer is to a set of reference points for each of the input variable. MES is then simply calculated as the minimum of its similarity with respect to each variable. Note that if a given raster cell has a negative IES value, it means that that value falls outside the range of values found in the reference set. A negative MESS thus represents sites where at least one variable has a value that is outside the range of environments over the reference set.
+<em>r.mess</em> computes the MES as well as layers that help to further interpret the MES values
+<ul>
+ <li>the individual environmental similarity surfaces (IES)</li>
+ <li>the MES (which is simply the mean of the IES)</li>
+ <li>the area where for at least one of the variables has a value that falls outside the range of values found in the reference set</li>
+ <li>the most dissimilar variable (MoD)</li>
+ <li>the mean of the IES layers where IES < 0</li>
+ <li>the number of layers with negative values</li>
+</ul>
-<p>The IES/MES compares the similarity between all raster cells in the current region and a subset of the raster cells in the region (the reference points). Any sample of interest can be used for the reference set. For example, it can be occurrence records for the species or a random sample of a region. It could also be the whole region/study area, e.g., if you want to compare current and future environmental conditions (see next point). The input layer representing the reference distribution can be a vector point layer or a raster layer.
+<p>The latter two are useful to distinguish areas where only one or few of the variables have values outside the range of the reference set and areas where many variables have values outside the range of the reference set.
-<p>Besides the reference set of environmental variables used to calculate the IES/MES, the user can define a second set of environmental variables (env1). This allows the user to compute how similar current conditions in given area are to those under e.g., future climates. This is for example useful to identify areas in your study area with novel future climates.
+<p>The reference points can be a binary raster layer (with 1 representing presence and 0 representing absence) or a vector point layer as reference points. Any sample of interest can be used for the reference set. Examples are points representing occurrence records for the species and areas that represent protected areas.
+
+<p>To compare e.g., current and future environmental conditions the user needs to define a reference set of environmental variables (env_old) and a set of future environmnental variables (env_new). This is for example used to identify areas with novel future climates.
-<p>In addition to the MES, which is the minimum IES, the r.mess function also allows to compute mean and median of the IES layers and a MoD layer, which identifies for each raster cells the IESS with the lowest value. Note that the mean and median of the IES should be used with care as negative and positive IES values have different meanings. In many cases it is probably more useful to calculate mean (or sum) over all negative. That would help to distinguish areas where only one or few of the variables have values outside the range of the reference set and areas where many variables have values outside the range of the reference set. Intuitively, one can see that conditions in the latter area are more dissimilar than in the former. I will probably add this option later, but the user can of course create such a layer by using the individual IES layers.
+<p>One can also test for the similarity between two different areas. For this one needs to provide a set of environmental variables (env_old) and a reference layer (ref_rast) for one region and a second set of environmental variables for another region (env_new).
-<h2>SEE ALSO</h2>
-There is also a similar function implemented for R using the R package <em>raster</em>, and which is part of the dismo package in R now. See <a href="http://rossijeanpierre.wordpress.com/2012/08/13/computing-the-multivariate-environmental-similarity-surfaces-mess-index-in-r/">here</a> for more information and the script.
<h2>REFERENCES</h2>
<p>Elith, J., Kearney, M., & Phillips, S. 2010. The art of modelling range-shifting species. Methods in Ecology and Evolution 1:330–342.
Added: grass-addons/grass7/raster/r.mess/r.mess.py
===================================================================
--- grass-addons/grass7/raster/r.mess/r.mess.py (rev 0)
+++ grass-addons/grass7/raster/r.mess/r.mess.py 2015-01-06 22:09:48 UTC (rev 63969)
@@ -0,0 +1,492 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+########################################################################
+#
+# 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).
+#
+#
+# COPYRIGHT: (C) 2014 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_old
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input (predictor) raster map(s)
+#% key_desc: names
+#% required: yes
+#% multiple: yes
+#% guisection: predictors
+#%end
+
+#%option
+#% key: env_new
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input (predictor) raster map(s)
+#% key_desc: names
+#% required: no
+#% multiple: yes
+#% guisection: predictors
+#%end
+
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Root name of the output MESS data layers
+#% key_desc: name
+#% required: yes
+#% guisection: Output
+#%end
+
+#%option
+#% key: digits
+#% type: integer
+#% description: Precision of your input layers values
+#% key_desc: string
+#% answer: 3
+#%end
+
+#%flag
+#% key: m
+#% description: Calculate Most dissimilar variable (MoD)
+#% guisection: Output
+#%end
+
+#%flag
+#% key: n
+#% description: Area with negative MESS
+#% guisection: Output
+#%end
+
+#%flag
+#% key: k
+#% description: mean(IES), where IES < 0
+#% guisection: Output
+#%end
+
+#%flag
+#% key: c
+#% description: Number of IES layers with values < 0
+#% guisection: Output
+#%end
+
+#%flag: IES
+#% key: i
+#% description: Remove individual environmental similarity layers (IES)
+#% guisection: Output
+#%end
+
+#%flag
+#% key: r
+#% description: Keep current region (otherwise set to match env_new)
+#%end
+
+#----------------------------------------------------------------------------
+# Standard
+#----------------------------------------------------------------------------
+
+# import libraries
+import os
+import sys
+import numpy as np
+import uuid
+import atexit
+import tempfile
+import operator
+import string
+import grass.script as grass
+from grass.script import db as db
+
+if not os.environ.has_key("GISBASE"):
+ grass.message( "You must be in GRASS GIS to run this program." )
+ sys.exit(1)
+
+# create set to store names of temporary maps to be deleted upon exit
+clean_rast = set()
+
+def cleanup():
+ for rast in clean_rast:
+ grass.run_command("g.remove",
+ type="rast", name = rast, quiet = True)
+
+# main function
+def main():
+
+ #----------------------------------------------------------------------------
+ # Variables
+ #----------------------------------------------------------------------------
+
+ #Test
+ #options = {"ref_rast":"bradypus", "ref_vect":"", "env_old":"bio1 at bradypus,bio5,bio6,bio7,bio8,bio12,bio16,bio17", "env_new":"bio1 at bradypus,bio5,bio6", "output":"MESSR", "digits":"5"}
+ #flags = {"m":True, "k":True, "n":False, "i":True, "k":True}
+
+
+ # reference layer
+ REF_VECT = options['ref_vect']
+ REF_RAST = options['ref_rast']
+ if len(REF_VECT) == 0 and len(REF_RAST) == 0:
+ grass.fatal("You did not provide the reference distribution layer")
+ else:
+ if len(REF_VECT) > 0 and len(REF_RAST) > 0:
+ grass.run_command('g.message', flags="w", quiet=True,
+ message = 'Two reference layers defined (vector and raster). Using the raster layer!')
+ vtl = REF_RAST
+ rtl = 'R'
+ else:
+ if len(REF_RAST) > 0:
+ vtl = REF_RAST
+ rtl = 'R'
+ else:
+ vtl = REF_VECT
+ rtl = 'V'
+
+ # old environmental layers & variable names
+ ipl = options['env_old']
+ ipl = ipl.split(',')
+ ipn = [i.split('@')[0] for i in ipl]
+ ipn = [x.lower() for x in ipn]
+
+ # new environmental variables
+ ipl2 = options['env_new']
+ if ipl2 == '':
+ ipl_dif = False
+ ipl2 = ipl
+ else:
+ ipl_dif = True
+ ipl2 = ipl2.split(',')
+ if len(ipl2) != len(ipl) and len(ipl2) != 0:
+ grass.fatal('number of old and new environmental variables is not the same')
+
+ # output layers
+ opl = options['output']
+ opc = opl + '_MESS'
+ ipi = [opl + '_' + i for i in ipn]
+
+ # flags
+ flm = flags['m']
+ flk = flags['k']
+ fln = flags['n']
+ il = flags['i']
+ flr = flags['r']
+ fll = flags['c']
+
+ # digits / precision
+ digits = int(options['digits'])
+ digits2 = pow(10, digits)
+
+ # Color table
+ tmpcol = tempfile.mkstemp()
+ text_file = open(tmpcol[1], "w")
+ text_file.write("0% 244:109:67\n")
+ text_file.write("0 255:255:210\n")
+ text_file.write("100% 50:136:189\n")
+ text_file.close()
+
+ # Check if there is a MASK
+ citiam = grass.find_file('MASK', element = 'cell')
+
+ #----------------------------------------------------------------------------
+ # Create the recode table - Reference distribution is raster
+ #----------------------------------------------------------------------------
+
+ if rtl=="R":
+
+ # Copy mask to temporary layer)
+ if citiam['fullname'] != '':
+ rname = "MASK" + str(uuid.uuid4())
+ rname = string.replace(rname, '-', '_')
+ grass.run_command('g.copy', quiet=True, raster = ('MASK', rname))
+ clean_rast.add(rname)
+
+ for i in xrange(len(ipl)):
+
+ # Create temporary layer based on reference layer
+ tmpf0 = "rmess_tmp_" + str(uuid.uuid4())
+ tmpf0 = string.replace(tmpf0, '-', '_')
+ grass.mapcalc("$tmpf0 = $vtl * 1", vtl = vtl, tmpf0=tmpf0, quiet=True)
+ clean_rast.add(tmpf0)
+
+ # Remove mask
+ if citiam['fullname'] != '':
+ grass.run_command("r.mask", quiet=True, flags="r")
+
+ # Create mask based on combined MASK/input layer
+ grass.run_command("r.mask", quiet=True, raster=tmpf0)
+
+ # Calculate the frequency distribution
+ tmpf1 = "rmess_tmp2_" + str(uuid.uuid4())
+ tmpf1 = string.replace(tmpf1, '-', '_')
+ grass.mapcalc("$tmpf1 = int($dignum * $inplay)",
+ tmpf1 = tmpf1,
+ inplay = ipl[i],
+ dignum = digits2,
+ quiet=True)
+ clean_rast.add(tmpf1)
+ p = grass.pipe_command('r.stats', quiet=True, flags = 'cn', input = tmpf1, sort = 'asc', sep=';')
+ stval = {}
+ for line in p.stdout:
+ [val,count] = line.strip(os.linesep).split(';')
+ stval[int(val)] = int(count)
+ p.wait()
+ sstval = sorted(stval.items(), key=operator.itemgetter(0))
+ sstval = np.matrix(sstval)
+ a = np.cumsum(np.array(sstval), axis=0)
+ b = np.sum(np.array(sstval), axis=0)
+ c = a[:,1] / b[1] * 100
+
+ # Restore mask and set region to new env_new
+ # Note: if region env_new is different from region env_old, setting
+ # the mask will not do anything. If there is partial overlap, the mask
+ # may affect the area where the two overlap
+ if citiam['fullname'] != '':
+ grass.run_command("r.mask", quiet=True, overwrite = True, raster=rname)
+ grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=rname)
+ else:
+ grass.run_command("r.mask", quiet=True, flags="r")
+ if ipl_dif and not flr:
+ grass.run_command("g.region", quiet=True, raster=ipl2[0])
+
+ # Get min and max values for recode table (based on full map)
+ tmpf2 = "rmess_tmp2_" + str(uuid.uuid4())
+ tmpf2 = string.replace(tmpf2, '-', '_')
+ grass.mapcalc("$tmpf2 = int($dignum * $inplay)",
+ tmpf2 = tmpf2,
+ inplay = ipl2[i],
+ dignum = digits2)
+ d = grass.parse_command("r.univar",quiet=True, flags="g", map=tmpf2)
+
+ # Create recode rules
+ Dmin = int(d['min'])
+ Dmax = int(d['max'])
+ envmin = np.min(np.array(sstval), axis=0)[0]
+ envmax = np.max(np.array(sstval), axis=0)[0]
+
+ if Dmin < envmin: e1 = Dmin - 1
+ else: e1 = envmin -1
+ if Dmax > envmax: e2 = Dmax + 1
+ else: e2 = envmax + 1
+
+ a1 = np.hstack([(e1), np.array(sstval.T[0])[0,:] ])
+ a2 = np.hstack([np.array(sstval.T[0])[0,:] -1, (e2)])
+ b1 = np.hstack([(0), c])
+
+ tmprule = tempfile.mkstemp()
+ text_file = open(tmprule[1], "w")
+ for k in np.arange(0,len(b1.T)):
+ rtmp = str(int(a1[k])) + ":" + str(int(a2[k])) + ":" + str(b1[k])
+ text_file.write(rtmp + "\n")
+ text_file.close()
+
+ # Create the recode layer and calculate the IES
+ tmpf3 = "rmess_tmp3_" + str(uuid.uuid4())
+ tmpf3 = string.replace(tmpf3, '-', '_')
+ grass.run_command("r.recode", input=tmpf2, output=tmpf3, rules=tmprule[1])
+
+ z1 = ipi[i] + " = if(" + tmpf3 + "==0, (float(" + tmpf2 + ")-" + str(float(envmin)) + ")/(" + str(float(envmax)) + "-" + str(float(envmin)) + ") * 100"
+ z2 = ", if(" + tmpf3 + "<=50, 2*float(" + tmpf3 + ")"
+ z3 = ", if(" + tmpf3 + "<100, 2*(100-float(" + tmpf3 + "))"
+ z4 = ", (" + str(float(envmax)) + "- float(" + tmpf2 + "))/(" + str(float(envmax)) + "-" + str(float(envmin)) + ") * 100.0)))"
+ calcc = z1 + z2 + z3 + z4
+ grass.mapcalc(calcc, quiet=True)
+ grass.run_command("r.colors", quiet=True, map=ipi[i], rules=tmpcol[1])
+ grass.run_command("g.remove", quiet=True, flags="f", type="raster", pattern=[tmpf0,tmpf1,tmpf2,tmpf3])
+ os.remove(tmprule[1])
+
+ # Recover original mask
+ if citiam['fullname'] != '':
+ grass.run_command("r.mask", quiet=True, raster=rname)
+ grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=rname)
+
+
+ #----------------------------------------------------------------------------
+ # Create the recode table - Reference distribution is vector
+ #----------------------------------------------------------------------------
+
+ if rtl=="V":
+
+ # Copy point layer and add columns for variables
+ tmpf0 = "rmess_tmp_" + str(uuid.uuid4())
+ tmpf0 = string.replace(tmpf0, '-', '_')
+ clean_rast.add(tmpf0)
+
+ grass.run_command("v.extract", quiet=True, flags="t", input=vtl, type="point", output=tmpf0)
+ grass.run_command("v.db.addtable", quiet=True, map=tmpf0)
+ grass.run_command("v.db.addcolumn", quiet=True, map=tmpf0, columns="envvar double precision")
+
+
+ # Upload raster values and get value in python as frequency table
+ check_n = len(np.hstack(db.db_select(sql = "SELECT cat FROM " + tmpf0)))
+ for m in xrange(len(ipl)):
+
+ grass.run_command("db.execute" , quiet=True, sql = "UPDATE " + tmpf0 + " SET envvar = NULL")
+ grass.run_command("v.what.rast", quiet=True, map=tmpf0, layer=1, raster=ipl[m], column="envvar")
+ volval = np.vstack(db.db_select(sql = "SELECT envvar,count(envvar) from " + tmpf0 +
+ " WHERE envvar IS NOT NULL GROUP BY envvar ORDER BY envvar"))
+ volval = volval.astype(np.float, copy=False)
+ a = np.cumsum(volval[:,1], axis=0)
+ b = np.sum(volval[:,1], axis=0)
+ c = a / b * 100
+
+ # Check for point without values
+ if b < check_n:
+ print("Please note that there were " + str(check_n - b) + " points without value")
+ print("This is probably because they lies outside the current region or input rasters")
+
+ # Set region to env_new layers (if different from env_old)
+ if ipl_dif and not flr:
+ grass.run_command("g.region", quiet=True, raster=ipl2[0])
+
+ # Multiply env layer with dignum
+ tmpf2 = "rmess_tmp2_" + str(uuid.uuid4())
+ tmpf2 = string.replace(tmpf2, '-', '_')
+ grass.mapcalc("$tmpf2 = int($dignum * $inplay)",
+ tmpf2 = tmpf2,
+ inplay = ipl2[m],
+ dignum = digits2,
+ quiet=True)
+
+ # Calculate min and max values of sample points and raster layer
+ envmin = int(min(volval[:,0]) * digits2)
+ envmax = int(max(volval[:,0]) * digits2)
+ Drange = grass.read_command("r.info", flags="r", map=tmpf2)
+ Drange = Drange.split('\n')
+ Drange = np.hstack([i.split('=') for i in Drange])
+ Dmin = int(Drange[1])
+ Dmax = int(Drange[3])
+
+ if Dmin < envmin: e1 = Dmin - 1
+ else: e1 = envmin -1
+ if Dmax > envmax: e2 = Dmax + 1
+ else: e2 = envmax + 1
+
+ a0 = volval[:,0] * digits2
+ a0 = a0.astype(np.int, copy=False)
+ a1 = np.hstack([(e1) , a0 ])
+ a2 = np.hstack([a0 -1, (e2) ])
+ b1 = np.hstack([(0), c])
+
+ tmprule = tempfile.mkstemp(suffix=ipn[m])
+ text_file = open(tmprule[1], "w")
+ for k in np.arange(0,len(b1)):
+ rtmp = str(int(a1[k])) + ":" + str(int(a2[k])) + ":" + str(b1[k])
+ text_file.write(rtmp + "\n")
+ text_file.close()
+
+ # Create the recode layer and calculate the IES
+ tmpf3 = "rmess_tmp3_" + str(uuid.uuid4())
+ tmpf3 = string.replace(tmpf3, '-', '_')
+ grass.run_command("r.recode", quiet=True, input=tmpf2, output=tmpf3, rules=tmprule[1])
+
+ z1 = ipi[m] + " = if(" + tmpf3 + "==0, (float(" + tmpf2 + ")-" + str(float(envmin)) + ")/(" + str(float(envmax)) + "-" + str(float(envmin)) + ") * 100"
+ z2 = ", if(" + tmpf3 + "<=50, 2*float(" + tmpf3 + ")"
+ z3 = ", if(" + tmpf3 + "<100, 2*(100-float(" + tmpf3 + "))"
+ z4 = ", (" + str(float(envmax)) + "- float(" + tmpf2 + "))/(" + str(float(envmax)) + "-" + str(float(envmin)) + ") * 100.0)))"
+ calcc = z1 + z2 + z3 + z4
+ grass.mapcalc(calcc, quiet=True)
+ grass.run_command("r.colors", quiet=True, map=ipi[m], rules=tmpcol[1])
+
+ # Clean up
+ grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=[tmpf2,tmpf3])
+ os.remove(tmprule[1])
+
+ # Clean up
+ grass.run_command("g.remove", quiet=True, flags="f", type="vector", name=tmpf0)
+
+
+ #-----------------------------------------------------------------------
+ # Calculate MESS statistics
+ #-----------------------------------------------------------------------
+
+ # MESS
+ grass.run_command("r.series", quiet=True, output=opc, input=ipi, method="minimum")
+
+ # Area with negative MESS
+ if fln:
+ grass.mapcalc(opc + "_neg = if(" + opc + "<0,1,0)", quiet=True)
+
+ # Most dissimilar variable (MoD)
+ if flm:
+ mod = opc + "_MoD"
+ grass.run_command("r.series", quiet=True, output=mod, input=ipi, method="min_raster")
+ tmpcat = tempfile.mkstemp()
+ text_file = open(tmpcat[1], "w")
+ for cats in xrange(len(ipi)):
+ text_file.write(str(cats) + ":" + ipi[cats] + "\n")
+ text_file.close()
+ grass.run_command("r.category", quiet=True, map=mod, rules=tmpcat[1], separator=":")
+
+ # mean(IES), where IES < 0
+ if flk:
+ MinMes = grass.read_command("r.info", quiet=True, flags="r", map=opc)
+ MinMes = MinMes.split('\n')
+ MinMes = float(np.hstack([i.split('=') for i in MinMes])[1])
+ mod = opc + "MeanNeg"
+ grass.run_command("r.series", quiet=True, input=ipi, output=mod, method="average", range=(MinMes,0))
+ grass.run_command("r.colors", quiet=True, map=mod, col="ryg")
+
+ # Number of layers with negative values
+ if fll:
+ MinMes = grass.read_command("r.info", quiet=True, flags="r", map=opc)
+ MinMes = MinMes.split('\n')
+ MinMes = float(np.hstack([i.split('=') for i in MinMes])[1])
+ mod = opc + "CountNeg"
+ grass.run_command("r.series", flags="n", quiet=True, input=ipi, output=mod, method="count", range=(MinMes,0))
+ grass.run_command("r.colors", quiet=True, map=mod, col="ryg")
+
+ # Remove IES layers
+ if il:
+ grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=ipi)
+
+ #=======================================================================
+ ## Clean up
+ #=======================================================================
+
+ os.remove(tmpcol[1])
+ os.remove(tmpcat[1])
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ atexit.register(cleanup)
+ sys.exit(main())
+
+
Property changes on: grass-addons/grass7/raster/r.mess/r.mess.py
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list