[GRASS-SVN] r59102 - grass-addons/grass7/raster/r.mess
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Feb 20 02:16:15 PST 2014
Author: pvanbosgeo
Date: 2014-02-20 02:16:14 -0800 (Thu, 20 Feb 2014)
New Revision: 59102
Modified:
grass-addons/grass7/raster/r.mess/r.mess
Log:
Added option to define vector with current and new set of environmental conditions. This will allow to calculate the MES between to time steps rather than between locations. If no set of environmental variables is given, the old ones are used.
Modified: grass-addons/grass7/raster/r.mess/r.mess
===================================================================
--- grass-addons/grass7/raster/r.mess/r.mess 2014-02-20 10:15:30 UTC (rev 59101)
+++ grass-addons/grass7/raster/r.mess/r.mess 2014-02-20 10:16:14 UTC (rev 59102)
@@ -59,7 +59,7 @@
#%end
#%option
-#% key: env_var
+#% key: env_old
#% type: string
#% gisprompt: old,cell,raster
#% description: Input (predictor) raster map(s)
@@ -69,6 +69,16 @@
#%end
#%option
+#% key: env_new
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input (predictor) raster map(s)
+#% key_desc: names
+#% required: no
+#% multiple: yes
+#%end
+
+#%option
#% key: output
#% type: string
#% gisprompt: new,cell,raster
@@ -123,7 +133,8 @@
## Set easier variable names
OUTMAPS="${GIS_OPT_OUTPUT}"
-export INMAPS="${GIS_OPT_ENV_VAR}"
+export INMAPS1="${GIS_OPT_ENV_OLD}"
+export INMAPS2="${GIS_OPT_ENV_NEW}"
#=======================================================================
## GRASS team recommandations
@@ -191,7 +202,7 @@
# Create vectors with names output maps [arrOUT]
# Create vector with names temporary layers [tmpLY]
-# Create vector with names input layers without mapset name [arrIN]
+# Create vector with names input layers without mapset name [arrIN1 and arrIN1]
# The 'removeThis' elements are to initiate the vector. How can I avoid this?
arrOUT="${OUTMAPS}_MESS"
@@ -199,7 +210,7 @@
arrIN="removeThis2"
counter=0
IFS=,
-for nvar in ${INMAPS} ; do
+for nvar in ${INMAPS1} ; do
export arrOUT="$arrOUT;${OUTMAPS}_`echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`"
counter=`expr $counter + 1`
export tmpLY="${tmpLY};tmp_mess_$$_$counter"
@@ -241,7 +252,7 @@
# test for output raster maps [2]
oIFS=$IFS
IFS=,
-for nvar in $INMAPS ; do
+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
@@ -255,7 +266,7 @@
# test for missing input raster maps
oIFS=$IFS
IFS=,
-for nvar in $INMAPS ; do
+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
@@ -306,8 +317,15 @@
## Get vector with variables
ipn <- Sys.getenv("arrIN") # variable names
ipn <- unlist(strsplit(ipn,";"))[-1]
-ipl <- Sys.getenv("INMAPS") # environmental layers
+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
@@ -327,7 +345,7 @@
#-----------------------------------------------------------------------
# Create the r.mapcalc expressions to calculate the mess for the
-# individual layers. The main step is to define the graph function to
+# individual layers. The main step is to define the recode rules to
# be used -
#-----------------------------------------------------------------------
@@ -340,10 +358,10 @@
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 the environmental data for the reference points 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.
+ # unique value is needed. Other option is to use freq() function in raster package
# 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. This may deal better with rounding errors too.
@@ -362,8 +380,12 @@
execGRASS("r.mask", input=vtl)
}
- execGRASS("r.out.xyz", input=ipl[i], output=flnm1, flags="overwrite")
- system(paste("awk -F \"|\" '{print $3}' ", flnm1, " > ", flnm2, sep=""))
+ execGRASS("r.out.xyz", input=ipl[i], output=flnm1, separator=",", flags="overwrite")
+ system(paste("awk -F \",\" '{print $3}' ", flnm1, " > ", flnm2, sep=""))
+
+ # Todo: see http://stackoverflow.com/questions/1727772/quickly-reading-very-large-tables-as-dataframes-in-r
+ # for alternative ways to import tables
+
spld <- scan(file=flnm2)
unlink(c(flnm1, flnm2))
spld <- round(spld, rdigits)
@@ -399,7 +421,7 @@
}
# Create the recode layer and calculate the IES
- execGRASS("r.recode", input=ipl[i], output=tml[i], rules=tmp.rule)
+ execGRASS("r.recode", input=ipl2[i], output=tml[i], rules=tmp.rule)
calcc <- paste(opi[i], " = if(", tml[i], "==0, (round(", ipl[i], "/",digits2,")*1.0-", round(envmin/digits2), ")/(", round(envmax/digits2)*1.0, "-", round(envmin/digits2)*1.0, ") *100.0, if(", tml[i], "<=50, 2.0 * ", tml[i], ", if(", tml[i], "<100, 2.0 * (100.0- ", tml[i], "), (", round(envmax/digits2)*1.0, "- round(", ipl[i], "/", digits2, ")*1.0)/(", round(envmax/digits2)*1.0, "-", round(envmin/digits2)*1.0, ") * 100.0)))", sep="")
execGRASS("r.mapcalc", expression=calcc, flags="overwrite")
execGRASS("g.remove", rast=tml[i])
@@ -464,7 +486,7 @@
digits2 <- 10^(-rdigits2)
# Create the recode layer and calculate the IES
- execGRASS("r.recode", input=ipl[i], output=tml[i], rules=tmp.rule)
+ execGRASS("r.recode", input=ipl2[i], output=tml[i], rules=tmp.rule)
calcc <- paste(opi[i], " = if(", tml[i], "==0, (round(", ipl[i], "/",digits2,")*1.0-", round(envmin/digits2), ")/(", round(envmax/digits2)*1.0, "-", round(envmin/digits2)*1.0, ") *100.0, if(", tml[i], "<=50, 2.0 * ", tml[i], ", if(", tml[i], "<100, 2.0 * (100.0- ", tml[i], "), (", round(envmax/digits2)*1.0, "- round(", ipl[i], "/", digits2, ")*1.0)/(", round(envmax/digits2)*1.0, "-", round(envmin/digits2)*1.0, ") * 100.0)))", sep="")
execGRASS("r.mapcalc", expression=calcc, flags="overwrite")
execGRASS("g.remove", rast=tml[i])
@@ -506,6 +528,7 @@
if(fll==1){
system(paste("r.series output=", opc, "_median input=", paste(opi, collapse=","), " method=median", sep=""))
}
+
EOF
}
@@ -527,7 +550,7 @@
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" "${INMAPS}" "$arrOUT" "$tmpLY" "$REF_LAY" "$GIS_FLAG_M" "$GIS_FLAG_K" "$GIS_FLAG_L" "$GIS_OPT_DIGITS" "$REF_TYPE"
+# 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
More information about the grass-commit
mailing list