[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