[GRASS-SVN] r53928 - grass-addons/grass6/raster/r.mess
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Nov 19 14:23:22 PST 2012
Author: pvanbosgeo
Date: 2012-11-19 14:23:22 -0800 (Mon, 19 Nov 2012)
New Revision: 53928
Modified:
grass-addons/grass6/raster/r.mess/r.mess
Log:
Added option to use a raster layer as reference distribution layer. Limited testing done
Modified: grass-addons/grass6/raster/r.mess/r.mess
===================================================================
--- grass-addons/grass6/raster/r.mess/r.mess 2012-11-19 22:09:35 UTC (rev 53927)
+++ grass-addons/grass6/raster/r.mess/r.mess 2012-11-19 22:23:22 UTC (rev 53928)
@@ -36,17 +36,29 @@
#%End
#%option
-#% key: vector_input
+#% 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 points
+#% description: Reference distribution as point vector layer
#% key_desc: name
-#% required: yes
+#% required: no
#% multiple: no
+#% guisection: reference distribution
#%end
#%option
-#% key: raster_input
+#% key: env_var
#% type: string
#% gisprompt: old,cell,raster
#% description: Input (predictor) raster map(s)
@@ -61,45 +73,37 @@
#% gisprompt: new
#% description: Root name of the output MESS data layers
#% key_desc: name
-#% required: no
+#% required: yes
#%end
#%option
-#% key: columns
-#% type: string
-#% description: Columns with environmental variables
-#% key_desc: string
-#% answer: Same names as input layers
-#% required: no
-#%end
-
-#%option
#% key: digits
#% type: string
-#% description: Level of precision of your data (default to 3 digits behind the commma)
+#% description: Precision of your input layers values
#% key_desc: string
#% answer: 0.001
-#% required: no
+#% required: yes
#%end
-##%flag
+##%flag # not implemented yet
##% key: n
-##% description: Keep individual environmental dissimilarity layers (IED)
+##% description: Keep individual environmental similarity layers (IES)
+##% guisection: Additional statistics
##%end
#%flag
#% key: m
-#% description: Most dissimilar variable (MoD)
+#% description: Calculate Most dissimilar variable (MoD)
#%end
#%flag
#% key: k
-#% description: Calculate MESS_mean
+#% description: Calculate mean of IES layers
#%end
#%flag
#% key: l
-#% description: Calculate MESS_median
+#% description: Calculate median of IES layers
#%end
#=======================================================================
@@ -145,7 +149,7 @@
## Set easier variable names
OUTMAPS="${GIS_OPT_OUTPUT}"
-INMAPS="${GIS_OPT_RASTER_INPUT}"
+INMAPS="${GIS_OPT_ENV_VAR}"
# test for output raster map
g.findfile element=cell file=${OUTMAPS}_combined > /dev/null
@@ -202,10 +206,30 @@
for nvar in ${INMAPS} ; do
arrOUT="$arrOUT;${OUTMAPS}_`echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`"
counter=`expr $counter + 1`
- tmpLY="${tmpLY};tmp.mess_$$_$counter"
+ 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
+
#=======================================================================
## Creating the R script
#=======================================================================
@@ -224,79 +248,112 @@
opi <- opl[-1] # base name individual layers
opc <- opl[1] # name of MESS layer
tml <- unlist(strsplit(args[4], ";"))[-1] # temporary layers
-vtl <- args[5] # vector layer
-cln <- args[6] # Column names if not the same as layer names
-digits <- as.numeric(args[10]) # Precision
+vtl <- args[5] # reference distribution raster layer
+rtl <- args[10] # raster or vector layer
+digits <- as.numeric(args[9]) # Precision
-# 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"
-#===============================================
-# 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 <- paste(ipn, collapse=",")
-}else{
- cln <- unlist(strsplit(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
+# be used -
+#-----------------------------------------------------------------------
-for(i in 1:length(ipl)){
- a <- table(spl[,i])
- envmin <- min(spl[,i])
- envmax <- max(spl[,i])
-
- x1 <- as.numeric(rownames(a))
- 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=",")
+if(rtl=="R"){
- system(paste("r.mapcalc '", tml[i], " = graph(", ipl[i], ",", a3, ")'", 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=""))
+ # Reference distribution layer is raster
+ #-----------------------------------------------------------------------
- # 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()
+ flnm1 <- system("g.tempfile -d pid=$$", intern=TRUE)
+ flnm2 <- paste(flnm1, ".txt", sep="")
+ 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
+
+ system("g.remove -f rast=MASK") # should I backup possible existing MASK?
+ 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))
+
+ # Calculate the frequency distribution
+ a <- table(spld)
+ envmin <- min(spld)
+ envmax <- max(spld)
+ rm(spld)
+ system("g.remove -f rast=MASK")
+
+ # Create the string with pairs of x-y values
+ x1 <- as.numeric(rownames(a))
+ x2 <- format(c(x1[-1],envmax+digits), scientific=FALSE, trim=TRUE)
+ x3 <- format(x1 + digits, scientific = FALSE, trim=TRUE)
+ a1 <- paste(x1,x3,sep=",", collapse=",")
+ b <- format(cumsum(as.numeric(a))/sum(a)*100, 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], ")-", 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=""))
+ }
+
+}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=""))
+ for(m in 1:length(ipn)){
+ system(paste("v.what.rast map=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 <- read.table(con, header=TRUE, sep="|")
+ 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])
+
+ # Create the string with pairs of x-y values
+ x1 <- as.numeric(rownames(a))
+ x2 <- format(c(x1[-1],envmax+digits), scientific=FALSE, trim=TRUE)
+ x3 <- format(x1 + digits, scientific = FALSE, trim=TRUE)
+ a1 <- paste(x1,x3,sep=",", collapse=",")
+ b <- format(cumsum(as.numeric(a))/sum(a)*100, 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], ")-", 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=""))
+ }
+
}
+#-----------------------------------------------------------------------
+# Calculate MESS
+#-----------------------------------------------------------------------
+
system(paste("r.series output=", opc, " input=", paste(opi, collapse=","), " method=minimum", sep=""))
-# Optionally, calculating extra layers (min_raster, mean and median)
-# Defines categories for min_raster
+#-----------------------------------------------------------------------
+# Calculate other stats
+#-----------------------------------------------------------------------
-if(args[7]==1){
+if(args[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=":")})
@@ -309,10 +366,10 @@
system(paste("r.category map=", opc, "_MoD rules=", tmpclas, sep=""))
unlink(tmpclas)
}
-if(args[8]==1){
+if(args[7]==1){
system(paste("r.series output=", opc, "_mean input=", paste(opi, collapse=","), " method=minimum", sep=""))
}
-if(args[9]==1){
+if(args[8]==1){
system(paste("r.series output=", opc, "_median input=", paste(opi, collapse=","), " method=median", sep=""))
}
@@ -338,13 +395,13 @@
##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" "${GIS_OPT_DIGITS}" < "$RGRASSSCRIPT" >> "$LOGFILE" 2>&1
+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"
+g.message "Finished - you can find the maps in your current mapset"
#=======================================================================
More information about the grass-commit
mailing list