[GRASS-SVN] r53849 - grass-addons/grass6/raster/r.mess
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Nov 16 09:45:18 PST 2012
Author: pvanbosgeo
Date: 2012-11-16 09:45:17 -0800 (Fri, 16 Nov 2012)
New Revision: 53849
Modified:
grass-addons/grass6/raster/r.mess/r.mess
Log:
Clean up of the code, can now be used in GRASS 7 too
Modified: grass-addons/grass6/raster/r.mess/r.mess
===================================================================
--- grass-addons/grass6/raster/r.mess/r.mess 2012-11-16 17:44:19 UTC (rev 53848)
+++ grass-addons/grass6/raster/r.mess/r.mess 2012-11-16 17:45:17 UTC (rev 53849)
@@ -13,7 +13,9 @@
# 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.
+# 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
@@ -47,7 +49,7 @@
#% key: raster_input
#% type: string
#% gisprompt: old,cell,raster
-#% description: Input raster map(s)
+#% description: Input (predictor) raster map(s)
#% key_desc: names
#% required: yes
#% multiple: yes
@@ -59,7 +61,7 @@
#% gisprompt: new
#% description: Root name of the output MESS data layers
#% key_desc: name
-#% required: yes
+#% required: no
#%end
#%option
@@ -71,6 +73,15 @@
#% required: no
#%end
+#%option
+#% key: digits
+#% type: string
+#% description: Level of precision of your data (default to 3 digits behind the commma)
+#% key_desc: string
+#% answer: 0.001
+#% required: no
+#%end
+
##%flag
##% key: n
##% description: Keep individual environmental dissimilarity layers (IED)
@@ -214,69 +225,70 @@
opc <- opl[1] # name of MESS layer
tml <- unlist(strsplit(args[4], ";"))[-1] # temporary layers
vtl <- args[5] # vector layer
-cln <- args[6]
+cln <- args[6] # Column names if not the same as layer names
+digits <- as.numeric(args[10]) # Precision
-# Import the vector point layer
-#spl <- readVECT6(vtl) # replace with below
-# To do: add option to select layer
-b <- execGRASS("v.db.select", parameters=list(map="vtl"), intern=TRUE)
-con <- textConnection(b)
-spl <- read.table(con, header=TRUE, sep="|")
-close(con)
+# for testing
+#==============================================
+#ipn <- c("bio_1", "bio_2")
+#ipl <- c("bio_1 at ConsStat", "bio_2 at ConsStat")
+#opl <- "AAA"
+#opi <- "AAA"
+#opc <- "AAAmess"
+#tml <- "tmp999999991"
+#vtl <- "AAtest at ConsStat"
+#digits <- 0.001
+#cln <- "Same names as input layers"
+#===============================================
-# Extract columns and create table
-# Use column names if given, otherwise,
-# use names input layers (ipl)
+# Import the vector attribute table. Use column names if given,
+# otherwise, # use names input layers (ipl)
+# To do: add option to select layer
if(cln == "Same names as input layers"){
- spld <- spl at data[,ipn]
+ spld <- paste(ipn, collapse=",")
}else{
cln <- unlist(strsplit(cln,","))
- spld <- spl at data
- spld <- spld[,cln]
+ spld <- cln
}
+b <- execGRASS("v.db.select", parameters=list(map=vtl, columns=spld), intern=TRUE)
+con <- textConnection(b)
+spl <- read.table(con, header=TRUE, sep="|")
+close(con)
#-----------------------------------------------------------------------
# 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 (check out: possible alternative would be to use reclass)
-if(is.null(dim(spld))){
- minp <- min(spld)
- maxp <- max(spld)
-}else{
- minp <- apply(spld, 2, min)
- maxp <- apply(spld, 2, max)
-}
+# be used
for(i in 1:length(ipl)){
- envir <- system(paste("r.univar -t map=", ipl[i], " fs=,", sep=""), intern=TRUE) # for use in grass7, change 'fs' to 'separator
- envir <- as.numeric(unlist(strsplit(envir[2], split=",")))
- envmax <- ifelse(envir[4]==maxp[i],maxp[i]+(maxp[i]-minp[i])*0.001,envir[4])
- envmin <- ifelse(envir[3]==minp[i],minp[i]-(maxp[i]-minp[i])*0.001,envir[3])
- if(is.null(dim(spld))){
- a <- table(spld)
- }else{a <- table(spld[,i])
- }
+ a <- table(spl[,i])
+ envmin <- min(spl[,i])
+ envmax <- max(spl[,i])
+
x1 <- as.numeric(rownames(a))
- x2 <- c(x1[-1],envmax)
- x3 <- x1 + (x2-x1)/1000
- a1 <- paste(envmin, paste(x1,x3,sep=",", collapse=","), envmax, sep=",")
- b <- round(cumsum(as.numeric(a))/sum(a)*100,3)
- a2 <- paste(0, 0, paste(b, b, sep=",", collapse=","), sep=",")
+ x2 <- format(c(x1[-1],envmax+digits), scientific=FALSE, trim=TRUE)
+ # Alternative: use r.mapcalc to create a reclass layer with all
+ # cells set to envmax if the environmental raster varsel[j] >= envmax.
+ # Next, use r.univar to find the minimum value of the reclass layer.
+ x3 <- format(x1 + digits, scientific = FALSE, trim=TRUE)
+ a1 <- paste(x1,x3,sep=",", collapse=",")
+ b <- format(round(cumsum(as.numeric(a))/sum(a)*100,3), scientific=FALSE, trim=TRUE)
+ a2 <- paste(c(0,b[-length(b)]), b, sep=",", collapse=",")
a3 <- paste(strsplit(a1,",")[[1]],strsplit(a2,",")[[1]],sep=",", collapse=",")
system(paste("r.mapcalc '", tml[i], " = graph(", ipl[i], ",", a3, ")'", sep=""))
- system(paste("r.mapcalc '", opi[i], " = if(", tml[i], "==0, (float(", ipl[i], ")-", minp[i], ")/(", maxp[i], "-", minp[i], ") *100.0, if(", tml[i], "<=50, 2 * ", tml[i], ", if(", tml[i], "<100.0, 2 * (100.0- ", tml[i], "), (", maxp[i], "- float(", ipl[i], "))/(", maxp[i], "-", minp[i], ") * 100.0)))'", 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=""))
# plotting the accumulative frequency plots for each bioclim variable
# need to add option in menu. Idea is to save the images, if this
# option is enabled, to default working directory
- #d1 <- paste(envir[3], paste(x1,x3,sep=",", collapse=","), envir[4], sep=",")
- #d2 <- paste(0, 0, paste(b, b, sep=",", collapse=","), sep=",")
- #png(paste(ipn, ".png", sep=""))
- #plot(strsplit(d1,",")[[1]],strsplit(d2,",")[[1]],type="l", xlab=ipn[i], ylab="% plots")
- #dev.off()
+ #d1 <- paste(envir[3], paste(x1,x3,sep=",", collapse=","), envir[4], sep=",")
+ #d2 <- paste(0, 0, paste(b, b, sep=",", collapse=","), sep=",")
+ #png(paste(ipn, ".png", sep=""))
+ #plot(strsplit(d1,",")[[1]],strsplit(d2,",")[[1]],type="l", xlab=ipn[i], ylab="% plots")
+ #dev.off()
}
system(paste("r.series output=", opc, " input=", paste(opi, collapse=","), " method=minimum", sep=""))
@@ -326,7 +338,7 @@
##using R to create MESS layers
#--slave
-R --vanilla --args "$arrIN" "${INMAPS}" "$arrOUT" "$tmpLY" "${GIS_OPT_VECTOR_INPUT}" "${GIS_OPT_COLUMNS}" "$GIS_FLAG_M" "$GIS_FLAG_K" "$GIS_FLAG_L" < "$RGRASSSCRIPT" >> "$LOGFILE" 2>&1
+R --vanilla --args "$arrIN" "${INMAPS}" "$arrOUT" "$tmpLY" "${GIS_OPT_VECTOR_INPUT}" "${GIS_OPT_COLUMNS}" "$GIS_FLAG_M" "$GIS_FLAG_K" "$GIS_FLAG_L" "${GIS_OPT_DIGITS}" < "$RGRASSSCRIPT" >> "$LOGFILE" 2>&1
if [ $? -ne 0 ] ; then
echo "ERROR: an error occurred during R script execution" 1>&2
exit 1
@@ -337,19 +349,3 @@
#=======================================================================
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
More information about the grass-commit
mailing list