[GRASS-SVN] r71043 - in grass-addons/grass7/raster: r.edm.eval r.maxent.lambdas r.stream.variables r.stream.watersheds
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat May 6 15:05:54 PDT 2017
Author: wenzeslaus
Date: 2017-05-06 15:05:54 -0700 (Sat, 06 May 2017)
New Revision: 71043
Added:
grass-addons/grass7/raster/r.edm.eval/r.edm.eval.sh
grass-addons/grass7/raster/r.maxent.lambdas/r.maxent.lambdas.sh
grass-addons/grass7/raster/r.stream.variables/r.stream.variables.sh
grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds.sh
Removed:
grass-addons/grass7/raster/r.edm.eval/r.edm.eval
grass-addons/grass7/raster/r.maxent.lambdas/r.maxent.lambdas
grass-addons/grass7/raster/r.stream.variables/r.stream.variables
grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds
Modified:
grass-addons/grass7/raster/r.edm.eval/Makefile
grass-addons/grass7/raster/r.maxent.lambdas/Makefile
grass-addons/grass7/raster/r.stream.variables/Makefile
grass-addons/grass7/raster/r.stream.watersheds/Makefile
Log:
use build support for shell scripts (depends on r71041)
Modified: grass-addons/grass7/raster/r.edm.eval/Makefile
===================================================================
--- grass-addons/grass7/raster/r.edm.eval/Makefile 2017-05-06 21:52:14 UTC (rev 71042)
+++ grass-addons/grass7/raster/r.edm.eval/Makefile 2017-05-06 22:05:54 UTC (rev 71043)
@@ -2,6 +2,6 @@
PGM = r.edm.eval
-include $(MODULE_TOPDIR)/include/Make/Script.make
+include $(MODULE_TOPDIR)/include/Make/ShScript.make
default: script
Deleted: grass-addons/grass7/raster/r.edm.eval/r.edm.eval
===================================================================
--- grass-addons/grass7/raster/r.edm.eval/r.edm.eval 2017-05-06 21:52:14 UTC (rev 71042)
+++ grass-addons/grass7/raster/r.edm.eval/r.edm.eval 2017-05-06 22:05:54 UTC (rev 71043)
@@ -1,484 +0,0 @@
-#!/bin/sh
-#
-#set -x
-########################################################################
-#
-# MODULE: r.model.eval
-# AUTHOR(S): Paulo van Breugel <paulo AT ecodiv.org>
-# PURPOSE: To evaluates how well a modeled distribution predicts an
-# observed distribution.
-#
-# NOTES: The observed distribution of e.g., a species, land
-# cover unit or vegetation unit should be a binary map
-# with 1 (present) and 0 (absent). The values of the
-# modeled distribution can be any map that represents a
-# probability distribution in space. This could be based
-# on X, but it doesn't need to be. You can also, for example,
-# evaluate how well the modeled distribution of a species X
-# predicts the distribution of species Y or of land cover
-# type Y.
-#
-# Disclaimer: Only limited testing has been done. Use at 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 evaluation statistics of an environmental distribution model, based on a layer with observed and a layer with predicted values
-#%End
-
-#%option
-#% key: obs
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Observed distribution
-#% key_desc: Raster name
-#% required: yes
-#% multiple: no
-#%end
-
-#%option
-#% key: mod
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Modeled distribution
-#% key_desc: Raster name
-#% required: yes
-#% multiple: no
-#%end
-
-#%flag
-#% key: b
-#% description: add presence to background?
-#%end
-
-#%option
-#% key: fstats
-#% type: string
-#% description: Name of output files (without extension)
-#% key_desc: File name
-#% required: yes
-#% multiple: no
-#%end
-
-#%option
-#% key: n_bins
-#% type: integer
-#% description: Number of bins in which to divide the modeled distribution scores
-#% key_desc: integer
-#% answer: 200
-#% required: no
-#%end
-
-#%flag
-#% key: l
-#% description: print log file to file?
-#%end
-
-#%option
-#% key: buffer_abs
-#% type: integer
-#% description: Restrict absence area to x km buffer
-#% key_desc: buffer zone (km)
-#% answer: 0
-#% required: no
-#% guisection: evaluation area
-#%end
-
-#%option
-#% key: buffer_pres
-#% type: integer
-#% description: Restrict presence area to x km buffer
-#% key_desc: buffer zone (km)
-#% answer: 0
-#% required: no
-#% guisection: evaluation area
-#%end
-
-#%option
-#% key: preval
-#% type: string
-#% description: Prevalence of presence points
-#% key_desc: <0-1>
-#% answer: 0
-#% required: no
-#% guisection: evaluation area
-#%end
-
-#%option
-#% key: num_pres
-#% type: integer
-#% description: number of presence points (this will not work if preval > 0)
-#% key_desc: integer
-#% answer: 0
-#% required: no
-#% guisection: evaluation area
-#%end
-
-#%option
-#% key: num_abs
-#% type: integer
-#% description: number of absence points (this will not work if preval > 0)
-#% key_desc: integer
-#% answer: 0
-#% required: no
-#% guisection: evaluation area
-#%end
-
-#=======================================================================
-## 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
-
-## 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
-
-#=======================================================================
-# test for missing input raster maps
-#=======================================================================
-
-oIFS=$IFS
-IFS=,
-arrIN=${GIS_OPT_OBS} `echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`
-g.findfile element=cell file=${arrIN} > /dev/null
-if [ $? -gt 0 ] ; then
- g.message -e 'The output map '${arrIN}' does not exists'
-exit 1
-fi
-unset arrIN
-
-arrIN=${GIS_OPT_MOD} `echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`
-g.findfile element=cell file=${arrIN} > /dev/null
-if [ $? -gt 0 ] ; then
- g.message -e 'The output map '${arrIN}' does not exists'
-exit 1
-fi
-IFS=$oIFS
-unset arrIN
-
-#=======================================================================
-## Creating the R script
-#=======================================================================
-
-writeScript(){
-cat > $1 << "EOF"
-
-options(echo = TRUE)
-library(rgrass7)
-library(reshape)
-library(caTools)
-
-# Define required functions
-get_output.GRASS <- function(x, separator=",", h=FALSE){
- con <- textConnection(x)
- MyVar <- read.table(con, header=h, sep=separator, comment.char="")
- close(con)
- return(MyVar)
-}
-
-rand.name <- function(base="tmp", l=6){
- paste(base,paste(sample(c(0:9, letters, LETTERS), l, replace=TRUE), collapse=""), sep="_")
-}
-
-copy_rast.GRASS <- function(rast, remove.orig=FALSE){
- mpst <- gmeta()$MAPSET
- fft <- execGRASS("g.findfile", flags="quiet", element="cell", file=rast, mapset=mpst, Sys_ignore.stdout=TRUE)
- if(fft==0){
- rname <- rand.name()
- execGRASS("g.copy", raster=paste(rast, rname, sep=","))
- if(remove.orig){execGRASS("g.remove", flags="f", type="raster", name=rast)}
- return(rname)
- }else{
- print("the selected file does not exist")
- }
-}
-
-## Get vector with variables
-obsmap <- Sys.getenv("GIS_OPT_OBS")
-distmod <- Sys.getenv("GIS_OPT_MOD")
-n.bins <- as.numeric(Sys.getenv("GIS_OPT_N_BINS"))
-bufabs <- as.numeric(Sys.getenv("GIS_OPT_BUFFER_ABS"))
-bufpres <- as.numeric(Sys.getenv("GIS_OPT_BUFFER_PRES"))
-fnames <- Sys.getenv("GIS_OPT_FSTATS")
-num.pres <- Sys.getenv("GIS_OPT_NUM_PRES")
-num.abs <- Sys.getenv("GIS_OPT_NUM_ABS")
-preval <- as.numeric(Sys.getenv("GIS_OPT_PREVAL"))
-print.log <- as.numeric(Sys.getenv("GIS_FLAG_L"))
-backgr <- as.numeric(Sys.getenv("GIS_FLAG_B"))
-
-# Create buffer map with buffer width of bufabs (outside) or bufpres (inside)
-if(bufabs>0){
- tmp.a <- "tmp_r_model_evaluation_0987654321a"
- tmp.b <- "tmp_r_model_evaluation_0987654321b"
- execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.a, obsmap, sep=" = "))
- execGRASS("r.null", map=tmp.a, setnull="0")
- execGRASS("r.buffer", input=tmp.a, output=tmp.b, distances=bufabs, units="kilometers")
- sink("reclass.rules"); cat("1:1:1:1"); cat("\n"); cat("2:2:0:0"); cat("\n"); sink()
- execGRASS("r.recode", flags="overwrite", input=tmp.b, output=tmp.a, rules="reclass.rules")
- obsmap <- tmp.a
- execGRASS("g.remove", flags="f", type="raster", name=tmp.b); unlink("reclass.rules")
-}
-
-if(bufpres>0){
- tmp.b <- "tmp_r_model_evaluation_0987654321b"
- tmp.c <- "tmp_r_model_evaluation_0987654321c"
- execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.b," = if(", obsmap,"==1,0,null())", sep=""))
- execGRASS("r.buffer", flags=c("z","overwrite"), input=tmp.b, output=tmp.c, distances=bufpres, units="kilometers")
- execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.b, " = if(", tmp.c, "==2,1,", obsmap, ")", sep=""))
- obsmap <- tmp.b
- execGRASS("g.remove", flags="f", type="raster", name=tmp.c)
-}
-
-# Prevalence
-if(preval>0 & preval < 1){
- tmp.e <- "tmp_r_model_evaluation_0987654321e"
- tmp.f <- "tmp_r_model_evaluation_0987654321f"
- tmp.g <- "tmp_r_model_evaluation_0987654321g"
- tmp.h <- "tmp_r_model_evaluation_0987654321h"
- tmp.i <- "tmp_r_model_evaluation_0987654321i"
- pa_cnt <- execGRASS("r.univar", flags=c("t","g"), map=obsmap, zones=obsmap, intern=TRUE)
- pa_cnt <- get_output.GRASS(pa_cnt, separator="|", h=TRUE)[,c(1,3)]
- pres.req <- ceiling((pa_cnt[1,2] + pa_cnt[2,2]) * preval)
- abs.req <- ceiling((pa_cnt[1,2] + pa_cnt[2,2]) - pres.req)
- if(pres.req > pa_cnt[2,2]){
- pres.req <- pa_cnt[2,2]
- abs.req <- ceiling((pa_cnt[2,2] / preval) - pa_cnt[2,2])
- }
- if(abs.req > pa_cnt[1,2]){
- pres.req <- ceiling(pa_cnt[1,2]/abs.req*pres.req)
- abs.req <- pa_cnt[1,2]
- }
- execGRASS("r.mapcalc", expression=paste(tmp.e, " = if(", obsmap, "== 1,1, null())", sep=""), Sys_wait=TRUE)
- execGRASS("r.random", input=tmp.e, n=as.character(pres.req), raster_output=tmp.f, Sys_wait=TRUE)
- execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.g, " = if(", obsmap, "==0,0,null())", sep=""), Sys_wait=TRUE)
- execGRASS("r.random", input=tmp.g, n=as.character(abs.req), raster_output=tmp.h, Sys_wait=TRUE)
- execGRASS("r.patch", input=paste(tmp.f, tmp.h, sep=","), output=tmp.i, flags="overwrite", Sys_wait=TRUE)
- obsmap <- tmp.i
- execGRASS("g.remove", flags="f", type="raster", name=paste(tmp.f, tmp.g, tmp.h, tmp.e, sep=","), Sys_wait=TRUE)
-}else{
- # Set number of presence and absence points
- if(num.pres > 0 | num.abs > 0){
- pa_cnt <- execGRASS("r.univar", flags=c("t","g"), map=obsmap, zones=obsmap, intern=TRUE)
- pa_cnt <- get_output.GRASS(pa_cnt, separator="|", h=TRUE)[,c(1,3)]
- num.pres <- ifelse(num.pres >= pa_cnt[2,2], 0, num.pres)
- num.abs <- ifelse(num.abs >= pa_cnt[1,2], 0, num.abs)
- }
- if(num.pres > 0){
- tmp.e <- "tmp_r_model_evaluation_0987654321e"
- tmp.f <- "tmp_r_model_evaluation_0987654321f"
- execGRASS("r.mapcalc", expression=paste(tmp.e, obsmap, sep=" = "), Sys_wait=TRUE)
- execGRASS("r.null", map=tmp.e, setnull="0", Sys_wait=TRUE)
- execGRASS("r.random", input=tmp.e, raster_output=tmp.f, n=as.character(num.pres), Sys_wait=TRUE)
- }
- if(num.abs > 0){
- tmp.g <- "tmp_r_model_evaluation_0987654321g"
- tmp.h <- "tmp_r_model_evaluation_0987654321h"
- execGRASS("r.mapcalc", expression=paste(tmp.g, obsmap, sep=" = "), Sys_wait=TRUE)
- execGRASS("r.null", map=tmp.g, setnull="1", Sys_wait=TRUE)
- execGRASS("r.random", input=tmp.g, raster_output=tmp.h, n=as.character(num.pres), Sys_wait=TRUE)
- }
- if(num.pres > 0 & num.abs == 0){
- tmp.h <- "tmp_r_model_evaluation_0987654321h"
- execGRASS("r.mapcalc", expression=paste(tmp.h, obsmap, sep=" = "), Sys_wait=TRUE)
- execGRASS("r.null", map=tmp.h, setnull="1", Sys_wait=TRUE)
- }
- if(num.abs > 0 & num.pres ==0){
- tmp.f <- "tmp_r_model_evaluation_0987654321f"
- execGRASS("r.mapcalc", expression=paste(tmp.f, obsmap, sep=" = "), Sys_wait=TRUE)
- execGRASS("r.null", map=tmp.f, setnull="0", Sys_wait=TRUE)
- }
- if(num.abs > 0 | num.pres > 0){
- tmp.j <- "tmp_r_model_evaluation_0987654321j"
- execGRASS("r.patch", input=paste(tmp.f, tmp.h, sep=","), output=tmp.j, flags="overwrite", Sys_wait=TRUE)
- obsmap <- tmp.j
- execGRASS("g.remove", flags="f", type="raster", name=paste(tmp.f, tmp.g, tmp.h, tmp.e, sep=","), Sys_wait=TRUE)
- }
-}
-
-# Create map with probability scores grouped in n.bins bins
-a <- execGRASS("r.univar", flags="g", map=distmod, intern=TRUE)
-maxmin <- get_output.GRASS(a, separator="=")
-stps <- (maxmin[5,2] - maxmin[4,2])/n.bins
-roc.st <- seq(maxmin[4,2], maxmin[5,2], by=stps)
-x1 <- c(0,roc.st[-length(roc.st)])
-x2 <- c(roc.st[1],roc.st[-1])
-y1 <- seq(length(roc.st)-1,0,-1)
-xy1 <- format(cbind(x1,x2,y1), scientific=F, trim=T, nsmall=12)
-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)
-tmp.d <- "tmp_r_model_evaluation_0987654321d"
-execGRASS("r.recode", flags="overwrite", input=distmod, output=tmp.d, rules=tmp.rule)
-unlink(tmp.rule)
-
-# Zonal stats
-if(backgr == 0){
- # Based on presence / absence
- a <- execGRASS("r.stats", flags=c("c","n"), input=paste(obsmap, tmp.d, sep=","), intern=TRUE)
- a <- get_output.GRASS(a, separator=" ")
- execGRASS("g.remove", flags="f", type="raster", name=tmp.d)
-}else{
- # Based on presence / background
- tmp.k <- "tmp_r_model_evaluation_0987654321k"
- execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.k, obsmap, sep=" = "), Sys_wait=TRUE)
- execGRASS("r.null", map=tmp.k, setnull="0", Sys_wait=TRUE)
- a1 <- execGRASS("r.stats", flags=c("c","n"), input=paste(tmp.k, tmp.d, sep=","), intern=TRUE)
- a1 <- get_output.GRASS(a1, separator=" ")
-
- tmp.m <- "tmp_r_model_evaluation_0987654321m"
- execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.m, obsmap, sep=" = "), Sys_wait=TRUE)
- execGRASS("r.null", map=tmp.m, setnull="1", Sys_wait=TRUE)
- a2 <- execGRASS("r.stats", flags=c("c","n"), input=paste(tmp.m, tmp.d, sep=","), intern=TRUE)
- a2 <- get_output.GRASS(a2, separator=" ")
- a3 <- cbind(a2[,-3],a1[,3] + a2[,3])
- names(a3) <- names(a1)
- a <- rbind(a1, a3)
-
- # Clean up
- execGRASS("g.remove", flags="f", type="raster", name=paste(tmp.k, tmp.m, sep=","))
-}
-
-# Clean up
-execGRASS("g.remove", flags="f", type="raster", name=tmp.d)
-if(bufabs>0){execGRASS("g.remove", flags="f", type="raster", name=tmp.a)}
-if(bufabs>0){execGRASS("g.remove", flags="f", type="raster", name=tmp.b)}
-if(preval>0){execGRASS("g.remove", flags="f", type="raster", name=tmp.i)}
-if(num.abs > 0 | num.pres > 0){execGRASS("g.remove", flags="f", type="raster", name=tmp.j)}
-
-# the temporary files do not always get cleaned up, need to check why
-# this is a temporary fix!
-aa <- execGRASS("g.list", type="raster", pattern="tmp_r_model_evaluation_*", intern=TRUE)
-if(length(aa)>0){
- execGRASS("g.remove", type="raster", pattern="tmp_r_model_evaluation_*", flags="f")
-}
-
-# Calculate evaluation stats
-a.cast <- cast(a, V2 ~ V1)
-a.cast[] <- sapply(a.cast, function(x) replace(x,is.na(x),0))
-a.sum <- apply(a.cast,2,sum)
-a.stats <- cbind(a.cast,cumsum(a.cast[,2]),cumsum(a.cast[,3]))[,-c(2,3)]
-names(a.stats) <- c("bin", "FP", "TP")
-a.stats$TN <- a.sum[1]-a.stats$FP
-a.stats$FN <- a.sum[2]-a.stats$TP
-a.stats$TPR <- apply(a.stats,1,function(x){x[3]/(x[3]+x[5])})
-a.stats$FPR <- apply(a.stats,1,function(x){x[2]/(x[2]+x[4])})
-a.stats$TNR <- apply(a.stats,1,function(x){x[4]/(x[2]+x[4])})
-a.stats$TSS <- apply(a.stats,1,function(x){x[6]+x[8]-1})
-PRa <- apply(a.stats,1,function(x){(x[3]+x[4])/sum(x)})
-Yobs <- apply(a.stats,1,function(x){(x[3]+x[5])/sum(x)})
-Ymod <- apply(a.stats,1,function(x){(x[3]+x[2])/sum(x)})
-PRe <- Yobs*Ymod + (1-Yobs)*(1-Ymod)
-a.stats$kappa <- (PRa - PRe) / (1 - PRe)
-n.presence <- a.stats$TP[2] + a.stats$FN[2]
-n.absence <- a.stats$TN[2] + a.stats$FP[2]
-prevalence <- n.presence / (n.presence + n.absence)
-
-# Calculate threshold values
-a.auc <- trapz(a.stats$FPR, a.stats$TPR)
-a.TSS <- max(a.stats$TSS)
-threshold.bin <- n.bins - a.stats[which.max(a.stats$TSS),1]
-a.TSS_threshold <- x2[threshold.bin]
-a.kappa_max <- max(a.stats$kappa)
-threshold.bin <- n.bins - a.stats[which.max(a.stats$kappa),1]
-a.kappa_threshold <- x2[threshold.bin]
-
-sink(paste(fnames, "summary.txt", sep="_"))
-cat(paste("AUC", as.character(a.auc), sep=" = ")); cat("\n")
-cat(paste("maximum TSS", a.TSS, sep=" = ")); cat("\n")
-cat(paste("maximum TSS threshold", a.TSS_threshold, sep=" = ")); cat("\n")
-cat(paste("maximum kappa", a.kappa_max, sep=" = ")); cat("\n")
-cat(paste("maximum kappa threshold", a.kappa_threshold, sep=" = ")); cat("\n")
-cat(paste("prevalence", prevalence, sep=" = ")); cat("\n")
-cat(paste("presence points", n.presence, sep=" = ")); cat("\n")
-cat(paste("absence points", n.absence, sep=" = ")); cat("\n")
-cat(paste("stats txt file =", paste(fnames, "summary.txt", sep="_"))); cat("\n")
-cat(paste("stats txt file =", paste(fnames, "stats.csv", sep="_"))); cat("\n")
-if(print.log==1){
- cat(paste("log file =", paste(fnames, "log.txt", sep="_"))); cat("\n")
-}
-sink()
-b.stats <- cbind(sort(x2[-length(x2)], decreasing=TRUE), a.stats)
-names(b.stats) <- c("threshold", names(a.stats))
-write.table(b.stats, file=paste(fnames, "stats.csv", sep="_"), append=TRUE, sep=";", row.names=FALSE)
-
-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='Working on it...'
-g.message message='-----'
-g.message message=''
-
-#using R to create MESS layers
-R CMD BATCH --slave "$RGRASSSCRIPT" ${GIS_OPT_FSTATS}_log.txt;
-
-if [ "$GIS_FLAG_L" -eq 0 ] ; then
- rm ${GIS_OPT_FSTATS}_log.txt
-fi
-
-head -11 ${GIS_OPT_FSTATS}_summary.txt
-
-g.message message='-----'
-g.message "Don't forget to check out the output files"
-
-#=======================================================================
-
Copied: grass-addons/grass7/raster/r.edm.eval/r.edm.eval.sh (from rev 71042, grass-addons/grass7/raster/r.edm.eval/r.edm.eval)
===================================================================
--- grass-addons/grass7/raster/r.edm.eval/r.edm.eval.sh (rev 0)
+++ grass-addons/grass7/raster/r.edm.eval/r.edm.eval.sh 2017-05-06 22:05:54 UTC (rev 71043)
@@ -0,0 +1,484 @@
+#!/bin/sh
+#
+#set -x
+########################################################################
+#
+# MODULE: r.model.eval
+# AUTHOR(S): Paulo van Breugel <paulo AT ecodiv.org>
+# PURPOSE: To evaluates how well a modeled distribution predicts an
+# observed distribution.
+#
+# NOTES: The observed distribution of e.g., a species, land
+# cover unit or vegetation unit should be a binary map
+# with 1 (present) and 0 (absent). The values of the
+# modeled distribution can be any map that represents a
+# probability distribution in space. This could be based
+# on X, but it doesn't need to be. You can also, for example,
+# evaluate how well the modeled distribution of a species X
+# predicts the distribution of species Y or of land cover
+# type Y.
+#
+# Disclaimer: Only limited testing has been done. Use at 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 evaluation statistics of an environmental distribution model, based on a layer with observed and a layer with predicted values
+#%End
+
+#%option
+#% key: obs
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Observed distribution
+#% key_desc: Raster name
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: mod
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Modeled distribution
+#% key_desc: Raster name
+#% required: yes
+#% multiple: no
+#%end
+
+#%flag
+#% key: b
+#% description: add presence to background?
+#%end
+
+#%option
+#% key: fstats
+#% type: string
+#% description: Name of output files (without extension)
+#% key_desc: File name
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: n_bins
+#% type: integer
+#% description: Number of bins in which to divide the modeled distribution scores
+#% key_desc: integer
+#% answer: 200
+#% required: no
+#%end
+
+#%flag
+#% key: l
+#% description: print log file to file?
+#%end
+
+#%option
+#% key: buffer_abs
+#% type: integer
+#% description: Restrict absence area to x km buffer
+#% key_desc: buffer zone (km)
+#% answer: 0
+#% required: no
+#% guisection: evaluation area
+#%end
+
+#%option
+#% key: buffer_pres
+#% type: integer
+#% description: Restrict presence area to x km buffer
+#% key_desc: buffer zone (km)
+#% answer: 0
+#% required: no
+#% guisection: evaluation area
+#%end
+
+#%option
+#% key: preval
+#% type: string
+#% description: Prevalence of presence points
+#% key_desc: <0-1>
+#% answer: 0
+#% required: no
+#% guisection: evaluation area
+#%end
+
+#%option
+#% key: num_pres
+#% type: integer
+#% description: number of presence points (this will not work if preval > 0)
+#% key_desc: integer
+#% answer: 0
+#% required: no
+#% guisection: evaluation area
+#%end
+
+#%option
+#% key: num_abs
+#% type: integer
+#% description: number of absence points (this will not work if preval > 0)
+#% key_desc: integer
+#% answer: 0
+#% required: no
+#% guisection: evaluation area
+#%end
+
+#=======================================================================
+## 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
+
+## 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
+
+#=======================================================================
+# test for missing input raster maps
+#=======================================================================
+
+oIFS=$IFS
+IFS=,
+arrIN=${GIS_OPT_OBS} `echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`
+g.findfile element=cell file=${arrIN} > /dev/null
+if [ $? -gt 0 ] ; then
+ g.message -e 'The output map '${arrIN}' does not exists'
+exit 1
+fi
+unset arrIN
+
+arrIN=${GIS_OPT_MOD} `echo $nvar | awk 'BEGIN{FS="@"}{print $1}'`
+g.findfile element=cell file=${arrIN} > /dev/null
+if [ $? -gt 0 ] ; then
+ g.message -e 'The output map '${arrIN}' does not exists'
+exit 1
+fi
+IFS=$oIFS
+unset arrIN
+
+#=======================================================================
+## Creating the R script
+#=======================================================================
+
+writeScript(){
+cat > $1 << "EOF"
+
+options(echo = TRUE)
+library(rgrass7)
+library(reshape)
+library(caTools)
+
+# Define required functions
+get_output.GRASS <- function(x, separator=",", h=FALSE){
+ con <- textConnection(x)
+ MyVar <- read.table(con, header=h, sep=separator, comment.char="")
+ close(con)
+ return(MyVar)
+}
+
+rand.name <- function(base="tmp", l=6){
+ paste(base,paste(sample(c(0:9, letters, LETTERS), l, replace=TRUE), collapse=""), sep="_")
+}
+
+copy_rast.GRASS <- function(rast, remove.orig=FALSE){
+ mpst <- gmeta()$MAPSET
+ fft <- execGRASS("g.findfile", flags="quiet", element="cell", file=rast, mapset=mpst, Sys_ignore.stdout=TRUE)
+ if(fft==0){
+ rname <- rand.name()
+ execGRASS("g.copy", raster=paste(rast, rname, sep=","))
+ if(remove.orig){execGRASS("g.remove", flags="f", type="raster", name=rast)}
+ return(rname)
+ }else{
+ print("the selected file does not exist")
+ }
+}
+
+## Get vector with variables
+obsmap <- Sys.getenv("GIS_OPT_OBS")
+distmod <- Sys.getenv("GIS_OPT_MOD")
+n.bins <- as.numeric(Sys.getenv("GIS_OPT_N_BINS"))
+bufabs <- as.numeric(Sys.getenv("GIS_OPT_BUFFER_ABS"))
+bufpres <- as.numeric(Sys.getenv("GIS_OPT_BUFFER_PRES"))
+fnames <- Sys.getenv("GIS_OPT_FSTATS")
+num.pres <- Sys.getenv("GIS_OPT_NUM_PRES")
+num.abs <- Sys.getenv("GIS_OPT_NUM_ABS")
+preval <- as.numeric(Sys.getenv("GIS_OPT_PREVAL"))
+print.log <- as.numeric(Sys.getenv("GIS_FLAG_L"))
+backgr <- as.numeric(Sys.getenv("GIS_FLAG_B"))
+
+# Create buffer map with buffer width of bufabs (outside) or bufpres (inside)
+if(bufabs>0){
+ tmp.a <- "tmp_r_model_evaluation_0987654321a"
+ tmp.b <- "tmp_r_model_evaluation_0987654321b"
+ execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.a, obsmap, sep=" = "))
+ execGRASS("r.null", map=tmp.a, setnull="0")
+ execGRASS("r.buffer", input=tmp.a, output=tmp.b, distances=bufabs, units="kilometers")
+ sink("reclass.rules"); cat("1:1:1:1"); cat("\n"); cat("2:2:0:0"); cat("\n"); sink()
+ execGRASS("r.recode", flags="overwrite", input=tmp.b, output=tmp.a, rules="reclass.rules")
+ obsmap <- tmp.a
+ execGRASS("g.remove", flags="f", type="raster", name=tmp.b); unlink("reclass.rules")
+}
+
+if(bufpres>0){
+ tmp.b <- "tmp_r_model_evaluation_0987654321b"
+ tmp.c <- "tmp_r_model_evaluation_0987654321c"
+ execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.b," = if(", obsmap,"==1,0,null())", sep=""))
+ execGRASS("r.buffer", flags=c("z","overwrite"), input=tmp.b, output=tmp.c, distances=bufpres, units="kilometers")
+ execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.b, " = if(", tmp.c, "==2,1,", obsmap, ")", sep=""))
+ obsmap <- tmp.b
+ execGRASS("g.remove", flags="f", type="raster", name=tmp.c)
+}
+
+# Prevalence
+if(preval>0 & preval < 1){
+ tmp.e <- "tmp_r_model_evaluation_0987654321e"
+ tmp.f <- "tmp_r_model_evaluation_0987654321f"
+ tmp.g <- "tmp_r_model_evaluation_0987654321g"
+ tmp.h <- "tmp_r_model_evaluation_0987654321h"
+ tmp.i <- "tmp_r_model_evaluation_0987654321i"
+ pa_cnt <- execGRASS("r.univar", flags=c("t","g"), map=obsmap, zones=obsmap, intern=TRUE)
+ pa_cnt <- get_output.GRASS(pa_cnt, separator="|", h=TRUE)[,c(1,3)]
+ pres.req <- ceiling((pa_cnt[1,2] + pa_cnt[2,2]) * preval)
+ abs.req <- ceiling((pa_cnt[1,2] + pa_cnt[2,2]) - pres.req)
+ if(pres.req > pa_cnt[2,2]){
+ pres.req <- pa_cnt[2,2]
+ abs.req <- ceiling((pa_cnt[2,2] / preval) - pa_cnt[2,2])
+ }
+ if(abs.req > pa_cnt[1,2]){
+ pres.req <- ceiling(pa_cnt[1,2]/abs.req*pres.req)
+ abs.req <- pa_cnt[1,2]
+ }
+ execGRASS("r.mapcalc", expression=paste(tmp.e, " = if(", obsmap, "== 1,1, null())", sep=""), Sys_wait=TRUE)
+ execGRASS("r.random", input=tmp.e, n=as.character(pres.req), raster_output=tmp.f, Sys_wait=TRUE)
+ execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.g, " = if(", obsmap, "==0,0,null())", sep=""), Sys_wait=TRUE)
+ execGRASS("r.random", input=tmp.g, n=as.character(abs.req), raster_output=tmp.h, Sys_wait=TRUE)
+ execGRASS("r.patch", input=paste(tmp.f, tmp.h, sep=","), output=tmp.i, flags="overwrite", Sys_wait=TRUE)
+ obsmap <- tmp.i
+ execGRASS("g.remove", flags="f", type="raster", name=paste(tmp.f, tmp.g, tmp.h, tmp.e, sep=","), Sys_wait=TRUE)
+}else{
+ # Set number of presence and absence points
+ if(num.pres > 0 | num.abs > 0){
+ pa_cnt <- execGRASS("r.univar", flags=c("t","g"), map=obsmap, zones=obsmap, intern=TRUE)
+ pa_cnt <- get_output.GRASS(pa_cnt, separator="|", h=TRUE)[,c(1,3)]
+ num.pres <- ifelse(num.pres >= pa_cnt[2,2], 0, num.pres)
+ num.abs <- ifelse(num.abs >= pa_cnt[1,2], 0, num.abs)
+ }
+ if(num.pres > 0){
+ tmp.e <- "tmp_r_model_evaluation_0987654321e"
+ tmp.f <- "tmp_r_model_evaluation_0987654321f"
+ execGRASS("r.mapcalc", expression=paste(tmp.e, obsmap, sep=" = "), Sys_wait=TRUE)
+ execGRASS("r.null", map=tmp.e, setnull="0", Sys_wait=TRUE)
+ execGRASS("r.random", input=tmp.e, raster_output=tmp.f, n=as.character(num.pres), Sys_wait=TRUE)
+ }
+ if(num.abs > 0){
+ tmp.g <- "tmp_r_model_evaluation_0987654321g"
+ tmp.h <- "tmp_r_model_evaluation_0987654321h"
+ execGRASS("r.mapcalc", expression=paste(tmp.g, obsmap, sep=" = "), Sys_wait=TRUE)
+ execGRASS("r.null", map=tmp.g, setnull="1", Sys_wait=TRUE)
+ execGRASS("r.random", input=tmp.g, raster_output=tmp.h, n=as.character(num.pres), Sys_wait=TRUE)
+ }
+ if(num.pres > 0 & num.abs == 0){
+ tmp.h <- "tmp_r_model_evaluation_0987654321h"
+ execGRASS("r.mapcalc", expression=paste(tmp.h, obsmap, sep=" = "), Sys_wait=TRUE)
+ execGRASS("r.null", map=tmp.h, setnull="1", Sys_wait=TRUE)
+ }
+ if(num.abs > 0 & num.pres ==0){
+ tmp.f <- "tmp_r_model_evaluation_0987654321f"
+ execGRASS("r.mapcalc", expression=paste(tmp.f, obsmap, sep=" = "), Sys_wait=TRUE)
+ execGRASS("r.null", map=tmp.f, setnull="0", Sys_wait=TRUE)
+ }
+ if(num.abs > 0 | num.pres > 0){
+ tmp.j <- "tmp_r_model_evaluation_0987654321j"
+ execGRASS("r.patch", input=paste(tmp.f, tmp.h, sep=","), output=tmp.j, flags="overwrite", Sys_wait=TRUE)
+ obsmap <- tmp.j
+ execGRASS("g.remove", flags="f", type="raster", name=paste(tmp.f, tmp.g, tmp.h, tmp.e, sep=","), Sys_wait=TRUE)
+ }
+}
+
+# Create map with probability scores grouped in n.bins bins
+a <- execGRASS("r.univar", flags="g", map=distmod, intern=TRUE)
+maxmin <- get_output.GRASS(a, separator="=")
+stps <- (maxmin[5,2] - maxmin[4,2])/n.bins
+roc.st <- seq(maxmin[4,2], maxmin[5,2], by=stps)
+x1 <- c(0,roc.st[-length(roc.st)])
+x2 <- c(roc.st[1],roc.st[-1])
+y1 <- seq(length(roc.st)-1,0,-1)
+xy1 <- format(cbind(x1,x2,y1), scientific=F, trim=T, nsmall=12)
+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)
+tmp.d <- "tmp_r_model_evaluation_0987654321d"
+execGRASS("r.recode", flags="overwrite", input=distmod, output=tmp.d, rules=tmp.rule)
+unlink(tmp.rule)
+
+# Zonal stats
+if(backgr == 0){
+ # Based on presence / absence
+ a <- execGRASS("r.stats", flags=c("c","n"), input=paste(obsmap, tmp.d, sep=","), intern=TRUE)
+ a <- get_output.GRASS(a, separator=" ")
+ execGRASS("g.remove", flags="f", type="raster", name=tmp.d)
+}else{
+ # Based on presence / background
+ tmp.k <- "tmp_r_model_evaluation_0987654321k"
+ execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.k, obsmap, sep=" = "), Sys_wait=TRUE)
+ execGRASS("r.null", map=tmp.k, setnull="0", Sys_wait=TRUE)
+ a1 <- execGRASS("r.stats", flags=c("c","n"), input=paste(tmp.k, tmp.d, sep=","), intern=TRUE)
+ a1 <- get_output.GRASS(a1, separator=" ")
+
+ tmp.m <- "tmp_r_model_evaluation_0987654321m"
+ execGRASS("r.mapcalc", flags="overwrite", expression=paste(tmp.m, obsmap, sep=" = "), Sys_wait=TRUE)
+ execGRASS("r.null", map=tmp.m, setnull="1", Sys_wait=TRUE)
+ a2 <- execGRASS("r.stats", flags=c("c","n"), input=paste(tmp.m, tmp.d, sep=","), intern=TRUE)
+ a2 <- get_output.GRASS(a2, separator=" ")
+ a3 <- cbind(a2[,-3],a1[,3] + a2[,3])
+ names(a3) <- names(a1)
+ a <- rbind(a1, a3)
+
+ # Clean up
+ execGRASS("g.remove", flags="f", type="raster", name=paste(tmp.k, tmp.m, sep=","))
+}
+
+# Clean up
+execGRASS("g.remove", flags="f", type="raster", name=tmp.d)
+if(bufabs>0){execGRASS("g.remove", flags="f", type="raster", name=tmp.a)}
+if(bufabs>0){execGRASS("g.remove", flags="f", type="raster", name=tmp.b)}
+if(preval>0){execGRASS("g.remove", flags="f", type="raster", name=tmp.i)}
+if(num.abs > 0 | num.pres > 0){execGRASS("g.remove", flags="f", type="raster", name=tmp.j)}
+
+# the temporary files do not always get cleaned up, need to check why
+# this is a temporary fix!
+aa <- execGRASS("g.list", type="raster", pattern="tmp_r_model_evaluation_*", intern=TRUE)
+if(length(aa)>0){
+ execGRASS("g.remove", type="raster", pattern="tmp_r_model_evaluation_*", flags="f")
+}
+
+# Calculate evaluation stats
+a.cast <- cast(a, V2 ~ V1)
+a.cast[] <- sapply(a.cast, function(x) replace(x,is.na(x),0))
+a.sum <- apply(a.cast,2,sum)
+a.stats <- cbind(a.cast,cumsum(a.cast[,2]),cumsum(a.cast[,3]))[,-c(2,3)]
+names(a.stats) <- c("bin", "FP", "TP")
+a.stats$TN <- a.sum[1]-a.stats$FP
+a.stats$FN <- a.sum[2]-a.stats$TP
+a.stats$TPR <- apply(a.stats,1,function(x){x[3]/(x[3]+x[5])})
+a.stats$FPR <- apply(a.stats,1,function(x){x[2]/(x[2]+x[4])})
+a.stats$TNR <- apply(a.stats,1,function(x){x[4]/(x[2]+x[4])})
+a.stats$TSS <- apply(a.stats,1,function(x){x[6]+x[8]-1})
+PRa <- apply(a.stats,1,function(x){(x[3]+x[4])/sum(x)})
+Yobs <- apply(a.stats,1,function(x){(x[3]+x[5])/sum(x)})
+Ymod <- apply(a.stats,1,function(x){(x[3]+x[2])/sum(x)})
+PRe <- Yobs*Ymod + (1-Yobs)*(1-Ymod)
+a.stats$kappa <- (PRa - PRe) / (1 - PRe)
+n.presence <- a.stats$TP[2] + a.stats$FN[2]
+n.absence <- a.stats$TN[2] + a.stats$FP[2]
+prevalence <- n.presence / (n.presence + n.absence)
+
+# Calculate threshold values
+a.auc <- trapz(a.stats$FPR, a.stats$TPR)
+a.TSS <- max(a.stats$TSS)
+threshold.bin <- n.bins - a.stats[which.max(a.stats$TSS),1]
+a.TSS_threshold <- x2[threshold.bin]
+a.kappa_max <- max(a.stats$kappa)
+threshold.bin <- n.bins - a.stats[which.max(a.stats$kappa),1]
+a.kappa_threshold <- x2[threshold.bin]
+
+sink(paste(fnames, "summary.txt", sep="_"))
+cat(paste("AUC", as.character(a.auc), sep=" = ")); cat("\n")
+cat(paste("maximum TSS", a.TSS, sep=" = ")); cat("\n")
+cat(paste("maximum TSS threshold", a.TSS_threshold, sep=" = ")); cat("\n")
+cat(paste("maximum kappa", a.kappa_max, sep=" = ")); cat("\n")
+cat(paste("maximum kappa threshold", a.kappa_threshold, sep=" = ")); cat("\n")
+cat(paste("prevalence", prevalence, sep=" = ")); cat("\n")
+cat(paste("presence points", n.presence, sep=" = ")); cat("\n")
+cat(paste("absence points", n.absence, sep=" = ")); cat("\n")
+cat(paste("stats txt file =", paste(fnames, "summary.txt", sep="_"))); cat("\n")
+cat(paste("stats txt file =", paste(fnames, "stats.csv", sep="_"))); cat("\n")
+if(print.log==1){
+ cat(paste("log file =", paste(fnames, "log.txt", sep="_"))); cat("\n")
+}
+sink()
+b.stats <- cbind(sort(x2[-length(x2)], decreasing=TRUE), a.stats)
+names(b.stats) <- c("threshold", names(a.stats))
+write.table(b.stats, file=paste(fnames, "stats.csv", sep="_"), append=TRUE, sep=";", row.names=FALSE)
+
+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='Working on it...'
+g.message message='-----'
+g.message message=''
+
+#using R to create MESS layers
+R CMD BATCH --slave "$RGRASSSCRIPT" ${GIS_OPT_FSTATS}_log.txt;
+
+if [ "$GIS_FLAG_L" -eq 0 ] ; then
+ rm ${GIS_OPT_FSTATS}_log.txt
+fi
+
+head -11 ${GIS_OPT_FSTATS}_summary.txt
+
+g.message message='-----'
+g.message "Don't forget to check out the output files"
+
+#=======================================================================
+
Modified: grass-addons/grass7/raster/r.maxent.lambdas/Makefile
===================================================================
--- grass-addons/grass7/raster/r.maxent.lambdas/Makefile 2017-05-06 21:52:14 UTC (rev 71042)
+++ grass-addons/grass7/raster/r.maxent.lambdas/Makefile 2017-05-06 22:05:54 UTC (rev 71043)
@@ -2,6 +2,6 @@
PGM = r.maxent.lambdas
-include $(MODULE_TOPDIR)/include/Make/Script.make
+include $(MODULE_TOPDIR)/include/Make/ShScript.make
default: script
Deleted: grass-addons/grass7/raster/r.maxent.lambdas/r.maxent.lambdas
===================================================================
--- grass-addons/grass7/raster/r.maxent.lambdas/r.maxent.lambdas 2017-05-06 21:52:14 UTC (rev 71042)
+++ grass-addons/grass7/raster/r.maxent.lambdas/r.maxent.lambdas 2017-05-06 22:05:54 UTC (rev 71043)
@@ -1,275 +0,0 @@
-#!/bin/sh
-#
-############################################################################
-#
-# MODULE: r.maxent.lambdas
-# AUTHOR(S): Stefan Blumentrath <stefan dot blumentrath at nina dot no >
-# Proposed small change in how raw features are extracted
-# from the lambdas file as this didn't work in original
-# code (Paulo van Breugel)
-# PURPOSE: Compute raw and/or logistic prediction maps from a lambdas
-# file produced with MaxEnt 3.3.3e.
-#
-# !!!This script works only if the input data to MaxEnt
-# are accessible from the current region.!!!
-#
-# This script will parse the specified lambdas-file from
-# MaxEnt 3.3.3e (see http://www.cs.princeton.edu/~schapire/maxent/)
-# and translate it into an r.mapcalc-expression which is then stored
-# in a temporary file and finally piped to r.mapcalc.
-# If alias names had been used in MaxEnt, these alias names can
-# automatically be replaced according to a CSV-like file provided
-# by the user. This file should contain alias names in the first
-# column and map names in the second column, seperated by comma,
-# without header. It should look e.g. like this:
-#
-# alias_1,map_1
-# alias_2,map_2
-# ...,...
-#
-# If such a CSV-file with alias names used in MaxEnt is provided,
-# the alias names from MaxEnt are replaced by map names.
-#
-# A raw output map is always computed from the MaxEnt model
-# as a first step. If logistic output is requested, the raw output
-# map can be deleted by the script ( using the l-flag). The
-# production of logistic output can be omitted. The logistic map
-# can be produced as an integer map. To do so the user has to specify
-# the number of digits after comma, that should be preserved in
-# integer output.
-# Optionally the map calculator expressions can be saved in a text
-# file, as especially the one for the raw output is likely to exceed
-# the space in the map history.
-# Due to conversion from double to floating-point in exp()-function, a
-# loss of precision from the 7th digit onwards is possible in the
-# logistic output.
-#
-# COPYRIGHT: (C) 2011 by the Norwegian Institute for Nature Research (NINA)
-# http://www.nina.no
-#
-# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
-# for details.
-#
-#############################################################################
-#
-# REQUIREMENTS:
-# awk
-#
-#%Module
-#% description: Computes raw and/or logistic prediction maps from MaxEnt lambdas files
-#%End
-#
-#%flag
-#% key: r
-#% description: Produce only raw output (both are computed by default).
-#%end
-#
-#%flag
-#% key: l
-#% description: Produce only logistic output (both are computed by default).
-#%end
-#
-#%option
-#% key: lambdas_file
-#% type: string
-#% description: MaxEnt lambdas-file to compute distribution-model from
-#% required : yes
-#% gisprompt: old_file,file,input
-#%end
-#
-#%option
-#% key: output_prefix
-#% type: string
-#% description: Prefix for output raster maps
-#% required : yes
-#% gisprompt: new,cell,raster
-#%end
-#
-#%option
-#% key: alias_file
-#% type: string
-#% description: CSV-file to replace alias names from MaxEnt by GRASS map names
-#% required : no
-#% gisprompt: old_file,file,input
-#%end
-#
-#%option
-#% key: integer_output
-#% type: integer
-#% description: Produce logistic integer output with this number of digits preserved
-#% required : no
-#% answer : 0
-#%end
-#
-#%option
-#% key: output_mapcalc
-#% type: string
-#% description: Save r.mapcalc expression to file
-#% required : no
-#% gisprompt: new_file,file,output
-#%end
-#
-#
-FLAG_ONLY_RAW=${GIS_FLAG_R}
-FLAG_ONLY_LOGISTIC=${GIS_FLAG_L}
-LAMBDAS_FILE="${GIS_OPT_LAMBDAS_FILE}"
-INTEGER_OUTPUT=${GIS_OPT_INTEGER_OUTPUT}
-#
-ALIAS_FILE="${GIS_OPT_ALIAS_FILE}"
-#
-OUTPUT_PREFIX="${GIS_OPT_OUTPUT_PREFIX}"
-OUTPUT_MAPCALC="${GIS_OPT_OUTPUT_MAPCALC}"
-#
-#Check if script is started from within GRASS
-if [ -z "$GISBASE" ] ; then
- echo "You must be in GRASS GIS to run this program." 1>&2
- exit 1
-fi
-#
-#Check if awk is installed
- if [ ! -x "`which awk`" ] ; then
- g.message -e "awk is required, please install awk or gawk first"
- exit 1
- fi
-#Set environment so that awk works properly in all languages
- unset LC_ALL
- LC_NUMERIC=C
- export LC_NUMERIC
-#
-#Pass evtl. command line arguments to gui and start it
-if [ "$1" != "@ARGS_PARSED@" ] ; then
- exec g.parser "$0" "$@"
-fi
-#
-#Check if input file exists
-if [ ! -r "$LAMBDAS_FILE" ] ; then
- g.message -e "MaxEnt lambdas-file could not be found or is not readable."
- exit 1
-fi
-#
-#Check if raw output files exists
-eval `g.findfile element=cell file=${OUTPUT_PREFIX}_raw`
-if [ -n "$name" ] ; then
- g.message -e "Raw output file already exists."
- exit 1
-fi
-#
-#Check if logistic output files exists when requestd
-if [ "${FLAG_ONLY_RAW}" -ne 1 ] ; then
- eval `g.findfile element=cell file=${OUTPUT_PREFIX}_logistic`
- if [ -n "$name" ] ; then
- g.message -e "Logistic output file already exists."
- exit 1
- fi
-fi
-#
-#Check if output flags are set properly
-if [ "${FLAG_ONLY_RAW}" -eq 1 -a "${FLAG_ONLY_LOGISTIC}" -eq 1 ] ; then
- g.message -e "The r and l flags are exclusive. Do not tick both!"
- rm "$temp1"
- exit 1
-fi
-#
-###Save r.mapcalc expression to file if requested, else to tmp-file
-if [ -z "$OUTPUT_MAPCALC" ] ; then
- #
- #Create tempfile
- temp1=`g.tempfile pid=$$`
- if [ $? -ne 0 ] || [ -z "$temp1" ] ; then
- g.message -e "ERROR: unable to create temporary files"
- exit 1
- fi
-else
- temp1="$OUTPUT_MAPCALC"
-fi
-#
-###Parse lambdas-file and translate it to a mapcalculator expression (safed in temporary file)
-###Get variables linearPredictorNormalizer, densityNormalizer and entropy from lambdas-file
-linearPredictorNormalizer=$(cat ${LAMBDAS_FILE} | grep linearPredictorNormalizer | sed 's/ //g' | tr -d '\r' | cut -f2 -d",")
-densityNormalizer=$(cat ${LAMBDAS_FILE} | grep densityNormalizer | sed 's/ //g' | tr -d '\r' | cut -f2 -d",")
-entropy=$(cat ${LAMBDAS_FILE} | grep entropy | sed 's/ //g' | tr -d '\r' | cut -f2 -d",")
-###Create file with initial part of mapcalc-expression
-echo "${OUTPUT_PREFIX}_raw = exp(((\\" > "$temp1"
-###Extract raw features
-cat -v "${LAMBDAS_FILE}" | awk '$0 !~ /\^|\(|\*/' | grep -vP "\`" | grep -vP "\'" | sed 's/\^M//g' | sed 's/,//g' | awk 'NF==4{print "if(isnull(" $1 "),0,(" $2 "*(" $1 "-" $3 ")/(" $4 "-" $3 ")))" "+" "\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
-###Extract quadratic features
-cat -v "${LAMBDAS_FILE}" | grep '\^2' | sed 's/\^M//g' | sed 's/,//g' | sed 's/\^2//g' | awk '{print "if(isnull(" $1 "),0,((" $2 "*(" $1 "*" $1 "-" $3 "))/" "(" $4 "-" $3 ")))" "+" "\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
-###Extract product features
-cat -v "${LAMBDAS_FILE}" | grep '*' | sed 's/\^M//g' | sed 's/,//g' | tr '*' ' ' | awk '{print "if((isnull(" $1 ")||isnull(" $2 ")),0,(" $3 "*(" $1 "*" $2 "-" $4 ")/(" $5 "-" $4 ")))+\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
-###Extract forward hinge features
-cat -v "${LAMBDAS_FILE}" | grep -P "\'" | sed 's/\^M//g' | sed 's/,//g' | cut -b 2- | tr '<' ' ' | awk '{print "if(isnull(" $1 "),0,if(" $1 ">=" $3 ",(" $2 "*(" $1 "-" $3 ")/(" $4 "-" $3 ")),0.0))+\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
-###Extract reverse hinge features
-cat -v "${LAMBDAS_FILE}" | grep -P "\`" | sed 's/\^M//g' | sed 's/,//g' | cut -b 2- | tr '<' ' ' | awk '{print "if(isnull(" $1 "),0,if(" $1 "<" $4 ",(" $2 "*(" $4 "-" $1 ")/(" $4 "-" $3 ")),0.0))+\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
-###Extract threshold features
-cat -v "${LAMBDAS_FILE}" | grep '(' | grep -v '=' | sed 's/\^M//g' | sed 's/,//g' | sed 's/(//g' | sed 's/)//g' | tr '<' ' ' | awk '{print "if(isnull(" $2 "),0,if(" $2 ">=" $1 "," $3 "))" "+\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
-###Extract categoric features
-cat -v "${LAMBDAS_FILE}" | grep '(' | grep '=' | sed 's/\^M//g' | sed 's/,//g' | sed 's/(//g' | sed 's/)//g' | tr '=' ' ' | awk '{print "if(isnull(" $1 "),0,if(" $1 "==" $2 "," $3 "))" "+\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
-###Replace last '\' by tail part of mapcalc-expression
-sed -i '$s/..$/\\/' "$temp1"
-echo ")-${linearPredictorNormalizer}))/${densityNormalizer}" >> "$temp1"
-###Replace possible repetitions of mathmatical signs
-sed -i 's/--/+/g' "$temp1"
-sed -i 's/+-/-/g' "$temp1"
-#
-#Check if alias file is provided and readable
-if [ -n "$ALIAS_FILE" ] ; then
- if [ ! -r "$ALIAS_FILE" ] ; then
- g.message -e "Alias file could not be found or is not readable."
- exit 1
- else
- #Parse alias-file to replace MaxEnt alias names by GRASS map names (including mapset specification)
- alias_names=$(cat $ALIAS_FILE | cut -f1 -d",")
- for a in $alias_names
- do
- m=$(cat $ALIAS_FILE | grep -w "${a}" | cut -f2 -d",")
- #
- g.message -v "Map name for alias $a is $m"
- #
- sed -i "s/(${a}-/(${m}-/" "$temp1"
- sed -i "s/(${a}+/(${m}+/" "$temp1"
- sed -i "s/(${a}>/(${m}>/" "$temp1"
- sed -i "s/(${a}</(${m}</" "$temp1"
- sed -i "s/(${a}=/(${m}=/" "$temp1"
- sed -i "s/(${a}\*/(${m}\*/" "$temp1"
- #
- sed -i "s/-${a})/-${m})/" "$temp1"
- sed -i "s/+${a})/+${m})/" "$temp1"
- sed -i "s/>${a})/>${m})/" "$temp1"
- sed -i "s/<${a})/<${m})/" "$temp1"
- sed -i "s/=${a})/=${m})/" "$temp1"
- sed -i "s/\*${a}-/\*${m}-/" "$temp1"
- sed -i "s/\*${a}+/\*${m}+/" "$temp1"
- done
- #
- fi
-fi
-#
-###Compute raw output map by sending expression saved in file temporary file to r.mapcalc
-cat "$temp1" | r.mapcalc
-#
-###Compute logistic output map if not suppressed
-if [ "${FLAG_ONLY_RAW}" -eq 0 ] ; then
- if [ "${INTEGER_OUTPUT}" -le 0 ] ; then
- echo "${OUTPUT_PREFIX}_logistic = (${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy})))" >> "$temp1"
- r.mapcalc "${OUTPUT_PREFIX}_logistic = (${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy})))"
- else
- if [ "${INTEGER_OUTPUT}" -lt 5 ] ; then
- echo "${OUTPUT_PREFIX}_logistic = round(((${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy}))))*(10^${INTEGER_OUTPUT}))" >> "$temp1"
- r.mapcalc "${OUTPUT_PREFIX}_logistic = round(((${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy}))))*(10^${INTEGER_OUTPUT}))"
- else
- echo "${OUTPUT_PREFIX}_logistic = round(((${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy}))))*100000.0)" >> "$temp1"
- r.mapcalc "${OUTPUT_PREFIX}_logistic = round(((${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy}))))*100000.0)"
- fi
- fi
-fi
-###Remove raw output map if requested
-if [ "${FLAG_ONLY_LOGISTIC}" -eq 1 ] ; then
- g.remove -f type='rast' name="${OUTPUT_PREFIX}"_raw
-fi
-###Remove tmp-file containing r.mapcalc expressions (if save to file not requested)
-if [ -z "$OUTPUT_MAPCALC" ] ; then
- rm -f "$temp1"
-fi
-
-g.message "Done"
Copied: grass-addons/grass7/raster/r.maxent.lambdas/r.maxent.lambdas.sh (from rev 71042, grass-addons/grass7/raster/r.maxent.lambdas/r.maxent.lambdas)
===================================================================
--- grass-addons/grass7/raster/r.maxent.lambdas/r.maxent.lambdas.sh (rev 0)
+++ grass-addons/grass7/raster/r.maxent.lambdas/r.maxent.lambdas.sh 2017-05-06 22:05:54 UTC (rev 71043)
@@ -0,0 +1,275 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: r.maxent.lambdas
+# AUTHOR(S): Stefan Blumentrath <stefan dot blumentrath at nina dot no >
+# Proposed small change in how raw features are extracted
+# from the lambdas file as this didn't work in original
+# code (Paulo van Breugel)
+# PURPOSE: Compute raw and/or logistic prediction maps from a lambdas
+# file produced with MaxEnt 3.3.3e.
+#
+# !!!This script works only if the input data to MaxEnt
+# are accessible from the current region.!!!
+#
+# This script will parse the specified lambdas-file from
+# MaxEnt 3.3.3e (see http://www.cs.princeton.edu/~schapire/maxent/)
+# and translate it into an r.mapcalc-expression which is then stored
+# in a temporary file and finally piped to r.mapcalc.
+# If alias names had been used in MaxEnt, these alias names can
+# automatically be replaced according to a CSV-like file provided
+# by the user. This file should contain alias names in the first
+# column and map names in the second column, seperated by comma,
+# without header. It should look e.g. like this:
+#
+# alias_1,map_1
+# alias_2,map_2
+# ...,...
+#
+# If such a CSV-file with alias names used in MaxEnt is provided,
+# the alias names from MaxEnt are replaced by map names.
+#
+# A raw output map is always computed from the MaxEnt model
+# as a first step. If logistic output is requested, the raw output
+# map can be deleted by the script ( using the l-flag). The
+# production of logistic output can be omitted. The logistic map
+# can be produced as an integer map. To do so the user has to specify
+# the number of digits after comma, that should be preserved in
+# integer output.
+# Optionally the map calculator expressions can be saved in a text
+# file, as especially the one for the raw output is likely to exceed
+# the space in the map history.
+# Due to conversion from double to floating-point in exp()-function, a
+# loss of precision from the 7th digit onwards is possible in the
+# logistic output.
+#
+# COPYRIGHT: (C) 2011 by the Norwegian Institute for Nature Research (NINA)
+# http://www.nina.no
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+# REQUIREMENTS:
+# awk
+#
+#%Module
+#% description: Computes raw and/or logistic prediction maps from MaxEnt lambdas files
+#%End
+#
+#%flag
+#% key: r
+#% description: Produce only raw output (both are computed by default).
+#%end
+#
+#%flag
+#% key: l
+#% description: Produce only logistic output (both are computed by default).
+#%end
+#
+#%option
+#% key: lambdas_file
+#% type: string
+#% description: MaxEnt lambdas-file to compute distribution-model from
+#% required : yes
+#% gisprompt: old_file,file,input
+#%end
+#
+#%option
+#% key: output_prefix
+#% type: string
+#% description: Prefix for output raster maps
+#% required : yes
+#% gisprompt: new,cell,raster
+#%end
+#
+#%option
+#% key: alias_file
+#% type: string
+#% description: CSV-file to replace alias names from MaxEnt by GRASS map names
+#% required : no
+#% gisprompt: old_file,file,input
+#%end
+#
+#%option
+#% key: integer_output
+#% type: integer
+#% description: Produce logistic integer output with this number of digits preserved
+#% required : no
+#% answer : 0
+#%end
+#
+#%option
+#% key: output_mapcalc
+#% type: string
+#% description: Save r.mapcalc expression to file
+#% required : no
+#% gisprompt: new_file,file,output
+#%end
+#
+#
+FLAG_ONLY_RAW=${GIS_FLAG_R}
+FLAG_ONLY_LOGISTIC=${GIS_FLAG_L}
+LAMBDAS_FILE="${GIS_OPT_LAMBDAS_FILE}"
+INTEGER_OUTPUT=${GIS_OPT_INTEGER_OUTPUT}
+#
+ALIAS_FILE="${GIS_OPT_ALIAS_FILE}"
+#
+OUTPUT_PREFIX="${GIS_OPT_OUTPUT_PREFIX}"
+OUTPUT_MAPCALC="${GIS_OPT_OUTPUT_MAPCALC}"
+#
+#Check if script is started from within GRASS
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." 1>&2
+ exit 1
+fi
+#
+#Check if awk is installed
+ if [ ! -x "`which awk`" ] ; then
+ g.message -e "awk is required, please install awk or gawk first"
+ exit 1
+ fi
+#Set environment so that awk works properly in all languages
+ unset LC_ALL
+ LC_NUMERIC=C
+ export LC_NUMERIC
+#
+#Pass evtl. command line arguments to gui and start it
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+#
+#Check if input file exists
+if [ ! -r "$LAMBDAS_FILE" ] ; then
+ g.message -e "MaxEnt lambdas-file could not be found or is not readable."
+ exit 1
+fi
+#
+#Check if raw output files exists
+eval `g.findfile element=cell file=${OUTPUT_PREFIX}_raw`
+if [ -n "$name" ] ; then
+ g.message -e "Raw output file already exists."
+ exit 1
+fi
+#
+#Check if logistic output files exists when requestd
+if [ "${FLAG_ONLY_RAW}" -ne 1 ] ; then
+ eval `g.findfile element=cell file=${OUTPUT_PREFIX}_logistic`
+ if [ -n "$name" ] ; then
+ g.message -e "Logistic output file already exists."
+ exit 1
+ fi
+fi
+#
+#Check if output flags are set properly
+if [ "${FLAG_ONLY_RAW}" -eq 1 -a "${FLAG_ONLY_LOGISTIC}" -eq 1 ] ; then
+ g.message -e "The r and l flags are exclusive. Do not tick both!"
+ rm "$temp1"
+ exit 1
+fi
+#
+###Save r.mapcalc expression to file if requested, else to tmp-file
+if [ -z "$OUTPUT_MAPCALC" ] ; then
+ #
+ #Create tempfile
+ temp1=`g.tempfile pid=$$`
+ if [ $? -ne 0 ] || [ -z "$temp1" ] ; then
+ g.message -e "ERROR: unable to create temporary files"
+ exit 1
+ fi
+else
+ temp1="$OUTPUT_MAPCALC"
+fi
+#
+###Parse lambdas-file and translate it to a mapcalculator expression (safed in temporary file)
+###Get variables linearPredictorNormalizer, densityNormalizer and entropy from lambdas-file
+linearPredictorNormalizer=$(cat ${LAMBDAS_FILE} | grep linearPredictorNormalizer | sed 's/ //g' | tr -d '\r' | cut -f2 -d",")
+densityNormalizer=$(cat ${LAMBDAS_FILE} | grep densityNormalizer | sed 's/ //g' | tr -d '\r' | cut -f2 -d",")
+entropy=$(cat ${LAMBDAS_FILE} | grep entropy | sed 's/ //g' | tr -d '\r' | cut -f2 -d",")
+###Create file with initial part of mapcalc-expression
+echo "${OUTPUT_PREFIX}_raw = exp(((\\" > "$temp1"
+###Extract raw features
+cat -v "${LAMBDAS_FILE}" | awk '$0 !~ /\^|\(|\*/' | grep -vP "\`" | grep -vP "\'" | sed 's/\^M//g' | sed 's/,//g' | awk 'NF==4{print "if(isnull(" $1 "),0,(" $2 "*(" $1 "-" $3 ")/(" $4 "-" $3 ")))" "+" "\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
+###Extract quadratic features
+cat -v "${LAMBDAS_FILE}" | grep '\^2' | sed 's/\^M//g' | sed 's/,//g' | sed 's/\^2//g' | awk '{print "if(isnull(" $1 "),0,((" $2 "*(" $1 "*" $1 "-" $3 "))/" "(" $4 "-" $3 ")))" "+" "\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
+###Extract product features
+cat -v "${LAMBDAS_FILE}" | grep '*' | sed 's/\^M//g' | sed 's/,//g' | tr '*' ' ' | awk '{print "if((isnull(" $1 ")||isnull(" $2 ")),0,(" $3 "*(" $1 "*" $2 "-" $4 ")/(" $5 "-" $4 ")))+\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
+###Extract forward hinge features
+cat -v "${LAMBDAS_FILE}" | grep -P "\'" | sed 's/\^M//g' | sed 's/,//g' | cut -b 2- | tr '<' ' ' | awk '{print "if(isnull(" $1 "),0,if(" $1 ">=" $3 ",(" $2 "*(" $1 "-" $3 ")/(" $4 "-" $3 ")),0.0))+\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
+###Extract reverse hinge features
+cat -v "${LAMBDAS_FILE}" | grep -P "\`" | sed 's/\^M//g' | sed 's/,//g' | cut -b 2- | tr '<' ' ' | awk '{print "if(isnull(" $1 "),0,if(" $1 "<" $4 ",(" $2 "*(" $4 "-" $1 ")/(" $4 "-" $3 ")),0.0))+\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
+###Extract threshold features
+cat -v "${LAMBDAS_FILE}" | grep '(' | grep -v '=' | sed 's/\^M//g' | sed 's/,//g' | sed 's/(//g' | sed 's/)//g' | tr '<' ' ' | awk '{print "if(isnull(" $2 "),0,if(" $2 ">=" $1 "," $3 "))" "+\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
+###Extract categoric features
+cat -v "${LAMBDAS_FILE}" | grep '(' | grep '=' | sed 's/\^M//g' | sed 's/,//g' | sed 's/(//g' | sed 's/)//g' | tr '=' ' ' | awk '{print "if(isnull(" $1 "),0,if(" $1 "==" $2 "," $3 "))" "+\\"}' | sed 's/--/+/g' | sed 's/+-/-/g' >> $temp1
+###Replace last '\' by tail part of mapcalc-expression
+sed -i '$s/..$/\\/' "$temp1"
+echo ")-${linearPredictorNormalizer}))/${densityNormalizer}" >> "$temp1"
+###Replace possible repetitions of mathmatical signs
+sed -i 's/--/+/g' "$temp1"
+sed -i 's/+-/-/g' "$temp1"
+#
+#Check if alias file is provided and readable
+if [ -n "$ALIAS_FILE" ] ; then
+ if [ ! -r "$ALIAS_FILE" ] ; then
+ g.message -e "Alias file could not be found or is not readable."
+ exit 1
+ else
+ #Parse alias-file to replace MaxEnt alias names by GRASS map names (including mapset specification)
+ alias_names=$(cat $ALIAS_FILE | cut -f1 -d",")
+ for a in $alias_names
+ do
+ m=$(cat $ALIAS_FILE | grep -w "${a}" | cut -f2 -d",")
+ #
+ g.message -v "Map name for alias $a is $m"
+ #
+ sed -i "s/(${a}-/(${m}-/" "$temp1"
+ sed -i "s/(${a}+/(${m}+/" "$temp1"
+ sed -i "s/(${a}>/(${m}>/" "$temp1"
+ sed -i "s/(${a}</(${m}</" "$temp1"
+ sed -i "s/(${a}=/(${m}=/" "$temp1"
+ sed -i "s/(${a}\*/(${m}\*/" "$temp1"
+ #
+ sed -i "s/-${a})/-${m})/" "$temp1"
+ sed -i "s/+${a})/+${m})/" "$temp1"
+ sed -i "s/>${a})/>${m})/" "$temp1"
+ sed -i "s/<${a})/<${m})/" "$temp1"
+ sed -i "s/=${a})/=${m})/" "$temp1"
+ sed -i "s/\*${a}-/\*${m}-/" "$temp1"
+ sed -i "s/\*${a}+/\*${m}+/" "$temp1"
+ done
+ #
+ fi
+fi
+#
+###Compute raw output map by sending expression saved in file temporary file to r.mapcalc
+cat "$temp1" | r.mapcalc
+#
+###Compute logistic output map if not suppressed
+if [ "${FLAG_ONLY_RAW}" -eq 0 ] ; then
+ if [ "${INTEGER_OUTPUT}" -le 0 ] ; then
+ echo "${OUTPUT_PREFIX}_logistic = (${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy})))" >> "$temp1"
+ r.mapcalc "${OUTPUT_PREFIX}_logistic = (${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy})))"
+ else
+ if [ "${INTEGER_OUTPUT}" -lt 5 ] ; then
+ echo "${OUTPUT_PREFIX}_logistic = round(((${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy}))))*(10^${INTEGER_OUTPUT}))" >> "$temp1"
+ r.mapcalc "${OUTPUT_PREFIX}_logistic = round(((${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy}))))*(10^${INTEGER_OUTPUT}))"
+ else
+ echo "${OUTPUT_PREFIX}_logistic = round(((${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy}))))*100000.0)" >> "$temp1"
+ r.mapcalc "${OUTPUT_PREFIX}_logistic = round(((${OUTPUT_PREFIX}_raw*exp(${entropy}))/(1.0+(${OUTPUT_PREFIX}_raw*exp(${entropy}))))*100000.0)"
+ fi
+ fi
+fi
+###Remove raw output map if requested
+if [ "${FLAG_ONLY_LOGISTIC}" -eq 1 ] ; then
+ g.remove -f type='rast' name="${OUTPUT_PREFIX}"_raw
+fi
+###Remove tmp-file containing r.mapcalc expressions (if save to file not requested)
+if [ -z "$OUTPUT_MAPCALC" ] ; then
+ rm -f "$temp1"
+fi
+
+g.message "Done"
Modified: grass-addons/grass7/raster/r.stream.variables/Makefile
===================================================================
--- grass-addons/grass7/raster/r.stream.variables/Makefile 2017-05-06 21:52:14 UTC (rev 71042)
+++ grass-addons/grass7/raster/r.stream.variables/Makefile 2017-05-06 22:05:54 UTC (rev 71043)
@@ -2,6 +2,6 @@
PGM=r.stream.variables
-include $(MODULE_TOPDIR)/include/Make/Script.make
+include $(MODULE_TOPDIR)/include/Make/ShScript.make
default: script
Deleted: grass-addons/grass7/raster/r.stream.variables/r.stream.variables
===================================================================
--- grass-addons/grass7/raster/r.stream.variables/r.stream.variables 2017-05-06 21:52:14 UTC (rev 71042)
+++ grass-addons/grass7/raster/r.stream.variables/r.stream.variables 2017-05-06 22:05:54 UTC (rev 71043)
@@ -1,395 +0,0 @@
-#!/bin/bash
-#
-############################################################################
-#
-# MODULE: r.stream.watersheds
-#
-# AUTHOR(S): Giuseppe Amatulli & Sami Domisch
-# based on "Domisch, S., Amatulli, G., Jetz, W. (in review)
-# Near-global freshwater-specific environmental variables for
-# biodiversity analyses in 1km resolution. Scientific Data"
-#
-# PURPOSE: Calculation of contiguous stream-specific variables that account
-# for the upstream environment (based on r.stream.watersheds).
-#
-# COPYRIGHT: (C) 2001-2012 by the GRASS Development Team
-#
-# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
-# for details.
-#
-#############################################################################
-
-#%Module
-#% description: Calculation of contiguous stream-specific variables that account for the upstream environment (based on r.stream.watersheds).
-#% keywords: drainage, stream, hydrology
-#%End
-
-#%option
-#% key: variable
-#% type: string
-#% key_desc: name
-#% description: Name of raster to be converted into a stream-specific variable
-#% required : yes
-#%end
-
-#%option
-#% key: area
-#% type: string
-#% key_desc: string
-#% multiple: no
-#% options: watershed,stream
-#% description: Area of aggregation: across the sub-watersheds or only across sub-streams
-#% required : yes
-#%end
-
-#%option
-#% key: folder
-#% type: string
-#% key_desc: name
-#% description: Provide the full folder path (same as for r.stream.watersheds)
-#% required:no
-#% answer: GISDBASE/folder_structure
-#%end
-
-#%option
-#% key: out_folder
-#% type: string
-#% key_desc: name
-#% description: Provide the full folder path for the output stream-specific variable
-#% required:no
-#% answer: GISDBASE/stream-specific_variables
-#%end
-
-#%option
-#% key: output
-#% type: string
-#% key_desc: method
-#% multiple: yes
-#% options: cells,min,max,range,mean,stddev,coeff_var,sum
-#% label: Provide the output aggregation method for the stream-specific variable
-#% description: upstream cells numbers, minimum, maximum, range, mean, standard deviation, coefficient of variation, sum. Output datatype is Int32
-#% required:yes
-#%end
-
-#%option
-#% key: scale
-#% type: double
-#% key_desc: value
-#% label: Provide a scale factor to multiply or divide the final stream-specific variable
-#% description: Provide it e.g. if input raster values are between -1 and 1, use scale=1000 to inicrease the number of decimals - all outputs will be rounded to integers
-#% answer: 1
-#% required:no
-#%end
-
-#%option
-#% key: cpu
-#% type: double
-#% description: Number of CPUs used for the parallel computation
-#% answer: 1
-#% required:no
-#%end
-
-if [ -z "$GISBASE" ] ; then
- echo "You must be in GRASS GIS to run this program." >&2
- exit 1
-fi
-
-if [ "$1" != "@ARGS_PARSED@" ] ; then
- exec g.parser "$0" "$@"
-fi
-
-#### check if we have awk
-if [ ! -x "`which awk`" ] ; then
- g.message -e "awk required, please install awk or gawk first"
- exit 1
-fi
-
-# setting environment, so that awk works properly in all languages
-unset LC_ALL
-LC_NUMERIC=C
-export LC_NUMERIC
-
-#test:
-
-if [ -z "$GIS_OPT_VARIABLE" ] ; then
- g.message "Please provide the name of raster to be converted into a stream-specific variable."
-exit 1
-fi
-
-if [ -z "$GIS_OPT_AREA" ] ; then
- g.message "Please provide area of aggregation: across the sub-watersheds or only across sub-streams."
-exit 1
-fi
-
-if [ -z "$GIS_OPT_OUTPUT" ] ; then
- g.message "Please provide the output aggregation method for the stream-specific variable: upstream min, max, range, mean, stddev, coeff_var, sum"
-exit 1
-fi
-
-#check if drainage and stream file exists
-eval `g.findfile element=cell file="$GIS_OPT_VARIABLE"`
-if [ ! "$file" ] ; then
- g.message -e "<$GIS_OPT_VARIABLE> does not exist! Aborting."
- exit 1
-fi
-
-
-export GISDBASE=$(g.gisenv get=GISDBASE separator=space)
-export LOCATION_NAME=$(g.gisenv get=LOCATION_NAME separator=space)
-export GIS_OPT_AREA
-export GIS_OPT_VARIABLE
-export GIS_OPT_OUTPUT
-export GIS_OPT_SCALE
-
-if [ ${GIS_OPT_FOLDER} = "GISDBASE/folder_structure" ] ; then
- export GIS_OPT_FOLDER=$GISDBASE"/folder_structure"
-else
- export GIS_OPT_FOLDER=$GIS_OPT_FOLDER
-fi
-
-if [ ${GIS_OPT_OUT_FOLDER} = "GISDBASE/stream-specific_variables" ] ; then
- export GIS_OPT_OUT_FOLDER=$GISDBASE"/stream-specific_variables"
-else
- export GIS_OPT_OUT_FOLDER=$GIS_OPT_OUT_FOLDER
-fi
-
-mkdir $GIS_OPT_OUT_FOLDER 2> /dev/null
-
-# what to do in case of user break:
-exitprocedure()
-{
-echo ""
-g.message -e 'Process aborted by user. All intermediate files have been deleted!'
-
-# reset in the permanent mapset
-g.gisenv set="MAPSET=PERMANENT"
-g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
-g.gisenv set="GISDBASE=$GISDBASE"
-
-# delete intermediate files
-rm -fr $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}*
-find $GIS_OPT_FOLDER/ -maxdepth 1 -name '*.txt' -delete
-find $GIS_OPT_FOLDER/ -maxdepth 1 -name '*.tif' -delete
-rm -fr $GIS_OPT_FOLDER/blockfile/stat_*.txt
-
-exit 1
-}
-
-# shell check for user break (signal list: trap -l)
-trap "exitprocedure" 2 3 9 15 19
-
-echo ""
-echo "##########################################################################"
-echo "Stream-specific variable calculation based on the sub-watersheds and sub-streams"
-echo ""
-echo "Citation: Domisch, S., Amatulli, G., Jetz, W. (in review)"
-echo "Near-global freshwater-specific environmental variables for"
-echo "biodiversity analyses in 1km resolution. Scientific Data"
-echo "##########################################################################"
-echo ""
-
-export GISRC_def=$GISRC
-
-rm -fr $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}*
-
-echo Using $( ls $GIS_OPT_FOLDER/blockfile/blockfile* | wc -l ) blocks in $GIS_OPT_FOLDER/blockfile/
-
-for file in $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/blockfile*_sub_${GIS_OPT_AREA}.tar.gz ; do
-
- export DIRNAME=$(dirname $file )
- filename=$(basename $file )
- cd $DIRNAME
- tar -xf $file
- cd $GISDBASE
-
- export BLKID=$(echo $filename | awk -v GIS_OPT_AREA=${GIS_OPT_AREA} '{ gsub("0000", "") ; gsub("blockfile","") ; gsub("_sub_","") ; gsub("GIS_OPT_AREA","") ; gsub(".tar.gz","") ; print }')
- export dir4d=${BLKID: -4:1} ; if [ -z $dir4d ] ; then dir4d=0 ; fi
- export dir3d=${BLKID: -3:1} ; if [ -z $dir3d ] ; then dir3d=0 ; fi
- export dir2d=${BLKID: -2:1} ; if [ -z $dir2d ] ; then dir2d=0 ; fi
- export dir1d=${BLKID: -1:1} ; if [ -z $dir1d ] ; then dir1d=0 ; fi
-
-g.gisenv set="GISDBASE=$GISDBASE"
-g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
-g.gisenv set="MAPSET=PERMANENT"
-
-echo "Start the stream-specific variable aggregation for block $file"
-
-ls $DIRNAME/sub_${GIS_OPT_AREA}ID*.tif | awk '{ if (NR>0) { print $1 , NR} }' | xargs -n 2 -P $GIS_OPT_CPU bash -c $'
-
-file=$1
-filename=$(basename $file .tif )
-if [ $GIS_OPT_AREA = "watershed" ] ; then ID=${filename:15:99}; fi
-if [ $GIS_OPT_AREA = "stream" ] ; then ID=${filename:12:99}; fi
-
-NR=$2
-
-if [ $( expr $( expr $NR + 1 ) / 100) -eq $( expr $NR / 100 + 1 ) ] ; then
-echo -en "\r$(expr $NR / 100 + 1 )% done"
-fi
-
-# replicate the start-GISRC in a uniq GISRC
-
-cp $GISRC_def $HOME/.grass7/rc$ID
-export GISRC=$HOME/.grass7/rc$ID
-
-g.mapset -c mapset=sub_${GIS_OPT_AREA}ID$ID location=$LOCATION_NAME dbase=$GISDBASE --quiet
-
-rm -f $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}ID$ID/.gislock
-
-export GISBASE=$( grep gisbase $(which grass70) | awk \'{ if(NR==2) { gsub ("\\"","" ) ; print $3 } }\' )
-export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
-export LD_LIBRARY_PATH="$GISBASE/lib"
-export GRASS_LD_LIBRARY_PATH="$LD_LIBRARY_PATH"
-export PYTHONPATH="$GISBASE/etc/python:$PYTHONPATH"
-export MANPATH=$MANPATH:$GISBASE/man
-
-# echo Load $DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif
-
-r.in.gdal input=$DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif output=sub_${GIS_OPT_AREA}ID${ID} --q
-
-g.region rast=sub_${GIS_OPT_AREA}ID${ID}@sub_${GIS_OPT_AREA}ID${ID} zoom=sub_${GIS_OPT_AREA}ID${ID}@sub_${GIS_OPT_AREA}ID${ID} --q
-r.mapcalc "sub_$GIS_OPT_VARIABLE = if ( sub_${GIS_OPT_AREA}ID${ID} == 1 , $GIS_OPT_VARIABLE at PERMANENT , null())" --o --q
-
-echo $ID"|"$( r.univar -t --q map=sub_$GIS_OPT_VARIABLE | awk \'{ if (NR==2 ) print}\' ) > $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt
-
-FULL=$(awk -F "|" \'{if (NF==13) {print 1 } else {print 0} }\' $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt)
-
-if [ $FULL -eq 0 ] ; then
-rm $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt
-else
-rm -f $DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif
-fi
-
-rm -r $HOME/.grass7/rc$ID $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}ID$ID
-
-' _
-
-g.gisenv set="MAPSET=PERMANENT"
-g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
-g.gisenv set="GISDBASE=$GISDBASE"
-
-echo -en "\r100% done"
-echo ""
-
-echo "Check for missing files due to RAM overload"
-
-if ls $DIRNAME/sub_${GIS_OPT_AREA}ID*.tif 1> /dev/null 2>&1 ; then
-
-echo $(ls $DIRNAME/sub_${GIS_OPT_AREA}ID*.tif | wc -l ) files had RAM issues
-
-ls $DIRNAME/sub_${GIS_OPT_AREA}ID*.tif | awk '{ if (NR>0) { print $1 , NR} }' | xargs -n 2 -P 1 bash -c $'
-
-file=$1
-filename=$(basename $file .tif )
-
-echo Fix missing file $file using only 1 CPU
-
-if [ $GIS_OPT_AREA = "watershed" ] ; then ID=${filename:15:99}; fi
-if [ $GIS_OPT_AREA = "stream" ] ; then ID=${filename:12:99}; fi
-
-NR=$2
-
-if [ ! -f $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt ] ; then
-
-if [ $( expr $( expr $NR + 1 ) / 100) -eq $( expr $NR / 100 + 1 ) ] ; then
-echo -en "\r$(expr $NR / 100 + 1 )% done"
-fi
-
-# replicate the start-GISRC in a uniq GISRC
-
-cp $GISRC_def $HOME/.grass7/rc$ID
-export GISRC=$HOME/.grass7/rc$ID
-
-g.mapset -c mapset=sub_${GIS_OPT_AREA}ID$ID location=$LOCATION_NAME dbase=$GISDBASE --quiet
-
-rm -f $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}ID$ID/.gislock
-
-export GISBASE=$( grep gisbase $(which grass70) | awk \'{ if(NR==2) { gsub ("\\"","" ) ; print $3 } }\' )
-export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
-export LD_LIBRARY_PATH="$GISBASE/lib"
-export GRASS_LD_LIBRARY_PATH="$LD_LIBRARY_PATH"
-export PYTHONPATH="$GISBASE/etc/python:$PYTHONPATH"
-export MANPATH=$MANPATH:$GISBASE/man
-
-# echo Load $DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif
-
-r.in.gdal input=$DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif output=sub_${GIS_OPT_AREA}ID${ID} --q
-rm -f $DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif
-
-g.region rast=sub_${GIS_OPT_AREA}ID${ID}@sub_${GIS_OPT_AREA}ID${ID} zoom=sub_${GIS_OPT_AREA}ID${ID}@sub_${GIS_OPT_AREA}ID${ID} --q
-r.mapcalc "sub_$GIS_OPT_VARIABLE = if ( sub_${GIS_OPT_AREA}ID${ID} == 1 , $GIS_OPT_VARIABLE at PERMANENT , null())" --o --q
-
-echo $ID"|"$( r.univar -t --q map=sub_$GIS_OPT_VARIABLE | awk \'{ if (NR==2 ) print}\' ) > $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt
-
-FULL=$(awk -F "|" \'{if (NF==13) {print 1 } else {print 0} }\' $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt)
-
-if [ $FULL -eq 0 ] ; then rm $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt ; fi
-
-rm -r $HOME/.grass7/rc$ID $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}ID$ID
-
-fi
-
-' _
-
-else
-echo 0 files had RAM issues
-fi
-
-
-cat $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID*.txt > $DIRNAME/stat_${GIS_OPT_VARIABLE}.txt
-rm -f $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID*.txt
-
-done
-
-echo "Aggregating the final table to reclassify the raster ID"
-
-echo "ID|non_null_cells|null_cells|min|max|range|mean|mean_of_abs|stddev|variance|coeff_var|sum|sum_abs" > $GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}.txt
-cat $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/stat_${GIS_OPT_VARIABLE}.txt >> $GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}.txt
-rm -f $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/stat_${GIS_OPT_VARIABLE}.txt
-
-echo Reclass the grid_id for the following output ${GIS_OPT_OUTPUT//","/" "}
-
-echo ${GIS_OPT_OUTPUT//","/" "} | xargs -n 1 -P $GIS_OPT_CPU bash -c $'
-
-if [ $1 = "cells" ] ; then col=2 ; fi
-if [ $1 = "min" ] ; then col=4 ; fi
-if [ $1 = "max" ] ; then col=5 ; fi
-if [ $1 = "range" ] ; then col=6 ; fi
-if [ $1 = "mean" ] ; then col=7 ; fi
-if [ $1 = "stddev" ] ; then col=9 ; fi
-if [ $1 = "coeff_var" ] ; then col=11 ; fi
-if [ $1 = "sum" ] ; then col=12 ; fi
-
-awk -F "|" -v col=$col -v scale=$GIS_OPT_SCALE \' { if (NR>1 ) {gsub("-nan","0",$col); gsub("inf","0",$col); print $1" = "int($col * scale) }}\' $GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}.txt > $GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}_$1.txt
-
-' _
-
-echo ${GIS_OPT_OUTPUT//","/" "} | xargs -n 1 -P 1 bash -c $'
-
-if [ $1 = "cells" ] ; then col=2 ; fi
-if [ $1 = "min" ] ; then col=4 ; fi
-if [ $1 = "max" ] ; then col=5 ; fi
-if [ $1 = "range" ] ; then col=6 ; fi
-if [ $1 = "mean" ] ; then col=7 ; fi
-if [ $1 = "stddev" ] ; then col=9 ; fi
-if [ $1 = "coeff_var" ] ; then col=11 ; fi
-if [ $1 = "sum" ] ; then col=12 ; fi
-
-echo Create stream-specific variable: $GIS_OPT_VARIABLE $1
-
-r.reclass input=grid_ID output=${GIS_OPT_VARIABLE}_$1 rules=$GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}_${1}.txt --overwrite --q
-# rm -f $GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}_${1}.txt
-
-r.mapcalc " ${GIS_OPT_VARIABLE}_$1 = ${GIS_OPT_VARIABLE}_$1 " --o
-
-r.out.gdal -c input=${GIS_OPT_VARIABLE}_${1} nodata=-9999 output=$GIS_OPT_OUT_FOLDER/${GIS_OPT_VARIABLE}_${1}.tif type=Int32 -c createopt="COMPRESS=LZW,ZLEVEL=9" --o --q 2> /dev/null
-if test -f "$GIS_OPT_FOLDER/blockfile/${GIS_OPT_VARIABLE}_${1}.tif"; then echo "$GIS_OPT_FOLDER/blockfile/${GIS_OPT_VARIABLE}_${1}.tif has been created!";fi
-#echo "$GIS_OPT_OUT_FOLDER/blockfile/${GIS_OPT_VARIABLE}_${1}.tif has been created!"
-
-' _
-
-exit
-
-echo "The full stream variable calculation has been done!"
-echo "To list the stream-specific output variables: ls $GIS_OPT_OUT_FOLDER/*.tif"
Copied: grass-addons/grass7/raster/r.stream.variables/r.stream.variables.sh (from rev 71042, grass-addons/grass7/raster/r.stream.variables/r.stream.variables)
===================================================================
--- grass-addons/grass7/raster/r.stream.variables/r.stream.variables.sh (rev 0)
+++ grass-addons/grass7/raster/r.stream.variables/r.stream.variables.sh 2017-05-06 22:05:54 UTC (rev 71043)
@@ -0,0 +1,395 @@
+#!/bin/bash
+#
+############################################################################
+#
+# MODULE: r.stream.watersheds
+#
+# AUTHOR(S): Giuseppe Amatulli & Sami Domisch
+# based on "Domisch, S., Amatulli, G., Jetz, W. (in review)
+# Near-global freshwater-specific environmental variables for
+# biodiversity analyses in 1km resolution. Scientific Data"
+#
+# PURPOSE: Calculation of contiguous stream-specific variables that account
+# for the upstream environment (based on r.stream.watersheds).
+#
+# COPYRIGHT: (C) 2001-2012 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+#%Module
+#% description: Calculation of contiguous stream-specific variables that account for the upstream environment (based on r.stream.watersheds).
+#% keywords: drainage, stream, hydrology
+#%End
+
+#%option
+#% key: variable
+#% type: string
+#% key_desc: name
+#% description: Name of raster to be converted into a stream-specific variable
+#% required : yes
+#%end
+
+#%option
+#% key: area
+#% type: string
+#% key_desc: string
+#% multiple: no
+#% options: watershed,stream
+#% description: Area of aggregation: across the sub-watersheds or only across sub-streams
+#% required : yes
+#%end
+
+#%option
+#% key: folder
+#% type: string
+#% key_desc: name
+#% description: Provide the full folder path (same as for r.stream.watersheds)
+#% required:no
+#% answer: GISDBASE/folder_structure
+#%end
+
+#%option
+#% key: out_folder
+#% type: string
+#% key_desc: name
+#% description: Provide the full folder path for the output stream-specific variable
+#% required:no
+#% answer: GISDBASE/stream-specific_variables
+#%end
+
+#%option
+#% key: output
+#% type: string
+#% key_desc: method
+#% multiple: yes
+#% options: cells,min,max,range,mean,stddev,coeff_var,sum
+#% label: Provide the output aggregation method for the stream-specific variable
+#% description: upstream cells numbers, minimum, maximum, range, mean, standard deviation, coefficient of variation, sum. Output datatype is Int32
+#% required:yes
+#%end
+
+#%option
+#% key: scale
+#% type: double
+#% key_desc: value
+#% label: Provide a scale factor to multiply or divide the final stream-specific variable
+#% description: Provide it e.g. if input raster values are between -1 and 1, use scale=1000 to inicrease the number of decimals - all outputs will be rounded to integers
+#% answer: 1
+#% required:no
+#%end
+
+#%option
+#% key: cpu
+#% type: double
+#% description: Number of CPUs used for the parallel computation
+#% answer: 1
+#% required:no
+#%end
+
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." >&2
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+#### check if we have awk
+if [ ! -x "`which awk`" ] ; then
+ g.message -e "awk required, please install awk or gawk first"
+ exit 1
+fi
+
+# setting environment, so that awk works properly in all languages
+unset LC_ALL
+LC_NUMERIC=C
+export LC_NUMERIC
+
+#test:
+
+if [ -z "$GIS_OPT_VARIABLE" ] ; then
+ g.message "Please provide the name of raster to be converted into a stream-specific variable."
+exit 1
+fi
+
+if [ -z "$GIS_OPT_AREA" ] ; then
+ g.message "Please provide area of aggregation: across the sub-watersheds or only across sub-streams."
+exit 1
+fi
+
+if [ -z "$GIS_OPT_OUTPUT" ] ; then
+ g.message "Please provide the output aggregation method for the stream-specific variable: upstream min, max, range, mean, stddev, coeff_var, sum"
+exit 1
+fi
+
+#check if drainage and stream file exists
+eval `g.findfile element=cell file="$GIS_OPT_VARIABLE"`
+if [ ! "$file" ] ; then
+ g.message -e "<$GIS_OPT_VARIABLE> does not exist! Aborting."
+ exit 1
+fi
+
+
+export GISDBASE=$(g.gisenv get=GISDBASE separator=space)
+export LOCATION_NAME=$(g.gisenv get=LOCATION_NAME separator=space)
+export GIS_OPT_AREA
+export GIS_OPT_VARIABLE
+export GIS_OPT_OUTPUT
+export GIS_OPT_SCALE
+
+if [ ${GIS_OPT_FOLDER} = "GISDBASE/folder_structure" ] ; then
+ export GIS_OPT_FOLDER=$GISDBASE"/folder_structure"
+else
+ export GIS_OPT_FOLDER=$GIS_OPT_FOLDER
+fi
+
+if [ ${GIS_OPT_OUT_FOLDER} = "GISDBASE/stream-specific_variables" ] ; then
+ export GIS_OPT_OUT_FOLDER=$GISDBASE"/stream-specific_variables"
+else
+ export GIS_OPT_OUT_FOLDER=$GIS_OPT_OUT_FOLDER
+fi
+
+mkdir $GIS_OPT_OUT_FOLDER 2> /dev/null
+
+# what to do in case of user break:
+exitprocedure()
+{
+echo ""
+g.message -e 'Process aborted by user. All intermediate files have been deleted!'
+
+# reset in the permanent mapset
+g.gisenv set="MAPSET=PERMANENT"
+g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
+g.gisenv set="GISDBASE=$GISDBASE"
+
+# delete intermediate files
+rm -fr $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}*
+find $GIS_OPT_FOLDER/ -maxdepth 1 -name '*.txt' -delete
+find $GIS_OPT_FOLDER/ -maxdepth 1 -name '*.tif' -delete
+rm -fr $GIS_OPT_FOLDER/blockfile/stat_*.txt
+
+exit 1
+}
+
+# shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 9 15 19
+
+echo ""
+echo "##########################################################################"
+echo "Stream-specific variable calculation based on the sub-watersheds and sub-streams"
+echo ""
+echo "Citation: Domisch, S., Amatulli, G., Jetz, W. (in review)"
+echo "Near-global freshwater-specific environmental variables for"
+echo "biodiversity analyses in 1km resolution. Scientific Data"
+echo "##########################################################################"
+echo ""
+
+export GISRC_def=$GISRC
+
+rm -fr $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}*
+
+echo Using $( ls $GIS_OPT_FOLDER/blockfile/blockfile* | wc -l ) blocks in $GIS_OPT_FOLDER/blockfile/
+
+for file in $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/blockfile*_sub_${GIS_OPT_AREA}.tar.gz ; do
+
+ export DIRNAME=$(dirname $file )
+ filename=$(basename $file )
+ cd $DIRNAME
+ tar -xf $file
+ cd $GISDBASE
+
+ export BLKID=$(echo $filename | awk -v GIS_OPT_AREA=${GIS_OPT_AREA} '{ gsub("0000", "") ; gsub("blockfile","") ; gsub("_sub_","") ; gsub("GIS_OPT_AREA","") ; gsub(".tar.gz","") ; print }')
+ export dir4d=${BLKID: -4:1} ; if [ -z $dir4d ] ; then dir4d=0 ; fi
+ export dir3d=${BLKID: -3:1} ; if [ -z $dir3d ] ; then dir3d=0 ; fi
+ export dir2d=${BLKID: -2:1} ; if [ -z $dir2d ] ; then dir2d=0 ; fi
+ export dir1d=${BLKID: -1:1} ; if [ -z $dir1d ] ; then dir1d=0 ; fi
+
+g.gisenv set="GISDBASE=$GISDBASE"
+g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
+g.gisenv set="MAPSET=PERMANENT"
+
+echo "Start the stream-specific variable aggregation for block $file"
+
+ls $DIRNAME/sub_${GIS_OPT_AREA}ID*.tif | awk '{ if (NR>0) { print $1 , NR} }' | xargs -n 2 -P $GIS_OPT_CPU bash -c $'
+
+file=$1
+filename=$(basename $file .tif )
+if [ $GIS_OPT_AREA = "watershed" ] ; then ID=${filename:15:99}; fi
+if [ $GIS_OPT_AREA = "stream" ] ; then ID=${filename:12:99}; fi
+
+NR=$2
+
+if [ $( expr $( expr $NR + 1 ) / 100) -eq $( expr $NR / 100 + 1 ) ] ; then
+echo -en "\r$(expr $NR / 100 + 1 )% done"
+fi
+
+# replicate the start-GISRC in a uniq GISRC
+
+cp $GISRC_def $HOME/.grass7/rc$ID
+export GISRC=$HOME/.grass7/rc$ID
+
+g.mapset -c mapset=sub_${GIS_OPT_AREA}ID$ID location=$LOCATION_NAME dbase=$GISDBASE --quiet
+
+rm -f $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}ID$ID/.gislock
+
+export GISBASE=$( grep gisbase $(which grass70) | awk \'{ if(NR==2) { gsub ("\\"","" ) ; print $3 } }\' )
+export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
+export LD_LIBRARY_PATH="$GISBASE/lib"
+export GRASS_LD_LIBRARY_PATH="$LD_LIBRARY_PATH"
+export PYTHONPATH="$GISBASE/etc/python:$PYTHONPATH"
+export MANPATH=$MANPATH:$GISBASE/man
+
+# echo Load $DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif
+
+r.in.gdal input=$DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif output=sub_${GIS_OPT_AREA}ID${ID} --q
+
+g.region rast=sub_${GIS_OPT_AREA}ID${ID}@sub_${GIS_OPT_AREA}ID${ID} zoom=sub_${GIS_OPT_AREA}ID${ID}@sub_${GIS_OPT_AREA}ID${ID} --q
+r.mapcalc "sub_$GIS_OPT_VARIABLE = if ( sub_${GIS_OPT_AREA}ID${ID} == 1 , $GIS_OPT_VARIABLE at PERMANENT , null())" --o --q
+
+echo $ID"|"$( r.univar -t --q map=sub_$GIS_OPT_VARIABLE | awk \'{ if (NR==2 ) print}\' ) > $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt
+
+FULL=$(awk -F "|" \'{if (NF==13) {print 1 } else {print 0} }\' $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt)
+
+if [ $FULL -eq 0 ] ; then
+rm $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt
+else
+rm -f $DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif
+fi
+
+rm -r $HOME/.grass7/rc$ID $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}ID$ID
+
+' _
+
+g.gisenv set="MAPSET=PERMANENT"
+g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
+g.gisenv set="GISDBASE=$GISDBASE"
+
+echo -en "\r100% done"
+echo ""
+
+echo "Check for missing files due to RAM overload"
+
+if ls $DIRNAME/sub_${GIS_OPT_AREA}ID*.tif 1> /dev/null 2>&1 ; then
+
+echo $(ls $DIRNAME/sub_${GIS_OPT_AREA}ID*.tif | wc -l ) files had RAM issues
+
+ls $DIRNAME/sub_${GIS_OPT_AREA}ID*.tif | awk '{ if (NR>0) { print $1 , NR} }' | xargs -n 2 -P 1 bash -c $'
+
+file=$1
+filename=$(basename $file .tif )
+
+echo Fix missing file $file using only 1 CPU
+
+if [ $GIS_OPT_AREA = "watershed" ] ; then ID=${filename:15:99}; fi
+if [ $GIS_OPT_AREA = "stream" ] ; then ID=${filename:12:99}; fi
+
+NR=$2
+
+if [ ! -f $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt ] ; then
+
+if [ $( expr $( expr $NR + 1 ) / 100) -eq $( expr $NR / 100 + 1 ) ] ; then
+echo -en "\r$(expr $NR / 100 + 1 )% done"
+fi
+
+# replicate the start-GISRC in a uniq GISRC
+
+cp $GISRC_def $HOME/.grass7/rc$ID
+export GISRC=$HOME/.grass7/rc$ID
+
+g.mapset -c mapset=sub_${GIS_OPT_AREA}ID$ID location=$LOCATION_NAME dbase=$GISDBASE --quiet
+
+rm -f $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}ID$ID/.gislock
+
+export GISBASE=$( grep gisbase $(which grass70) | awk \'{ if(NR==2) { gsub ("\\"","" ) ; print $3 } }\' )
+export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
+export LD_LIBRARY_PATH="$GISBASE/lib"
+export GRASS_LD_LIBRARY_PATH="$LD_LIBRARY_PATH"
+export PYTHONPATH="$GISBASE/etc/python:$PYTHONPATH"
+export MANPATH=$MANPATH:$GISBASE/man
+
+# echo Load $DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif
+
+r.in.gdal input=$DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif output=sub_${GIS_OPT_AREA}ID${ID} --q
+rm -f $DIRNAME/sub_${GIS_OPT_AREA}ID${ID}.tif
+
+g.region rast=sub_${GIS_OPT_AREA}ID${ID}@sub_${GIS_OPT_AREA}ID${ID} zoom=sub_${GIS_OPT_AREA}ID${ID}@sub_${GIS_OPT_AREA}ID${ID} --q
+r.mapcalc "sub_$GIS_OPT_VARIABLE = if ( sub_${GIS_OPT_AREA}ID${ID} == 1 , $GIS_OPT_VARIABLE at PERMANENT , null())" --o --q
+
+echo $ID"|"$( r.univar -t --q map=sub_$GIS_OPT_VARIABLE | awk \'{ if (NR==2 ) print}\' ) > $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt
+
+FULL=$(awk -F "|" \'{if (NF==13) {print 1 } else {print 0} }\' $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt)
+
+if [ $FULL -eq 0 ] ; then rm $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID$ID.txt ; fi
+
+rm -r $HOME/.grass7/rc$ID $GISDBASE/$LOCATION_NAME/sub_${GIS_OPT_AREA}ID$ID
+
+fi
+
+' _
+
+else
+echo 0 files had RAM issues
+fi
+
+
+cat $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID*.txt > $DIRNAME/stat_${GIS_OPT_VARIABLE}.txt
+rm -f $DIRNAME/stat_${GIS_OPT_VARIABLE}_ID*.txt
+
+done
+
+echo "Aggregating the final table to reclassify the raster ID"
+
+echo "ID|non_null_cells|null_cells|min|max|range|mean|mean_of_abs|stddev|variance|coeff_var|sum|sum_abs" > $GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}.txt
+cat $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/stat_${GIS_OPT_VARIABLE}.txt >> $GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}.txt
+rm -f $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/stat_${GIS_OPT_VARIABLE}.txt
+
+echo Reclass the grid_id for the following output ${GIS_OPT_OUTPUT//","/" "}
+
+echo ${GIS_OPT_OUTPUT//","/" "} | xargs -n 1 -P $GIS_OPT_CPU bash -c $'
+
+if [ $1 = "cells" ] ; then col=2 ; fi
+if [ $1 = "min" ] ; then col=4 ; fi
+if [ $1 = "max" ] ; then col=5 ; fi
+if [ $1 = "range" ] ; then col=6 ; fi
+if [ $1 = "mean" ] ; then col=7 ; fi
+if [ $1 = "stddev" ] ; then col=9 ; fi
+if [ $1 = "coeff_var" ] ; then col=11 ; fi
+if [ $1 = "sum" ] ; then col=12 ; fi
+
+awk -F "|" -v col=$col -v scale=$GIS_OPT_SCALE \' { if (NR>1 ) {gsub("-nan","0",$col); gsub("inf","0",$col); print $1" = "int($col * scale) }}\' $GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}.txt > $GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}_$1.txt
+
+' _
+
+echo ${GIS_OPT_OUTPUT//","/" "} | xargs -n 1 -P 1 bash -c $'
+
+if [ $1 = "cells" ] ; then col=2 ; fi
+if [ $1 = "min" ] ; then col=4 ; fi
+if [ $1 = "max" ] ; then col=5 ; fi
+if [ $1 = "range" ] ; then col=6 ; fi
+if [ $1 = "mean" ] ; then col=7 ; fi
+if [ $1 = "stddev" ] ; then col=9 ; fi
+if [ $1 = "coeff_var" ] ; then col=11 ; fi
+if [ $1 = "sum" ] ; then col=12 ; fi
+
+echo Create stream-specific variable: $GIS_OPT_VARIABLE $1
+
+r.reclass input=grid_ID output=${GIS_OPT_VARIABLE}_$1 rules=$GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}_${1}.txt --overwrite --q
+# rm -f $GIS_OPT_FOLDER/blockfile/stat_${GIS_OPT_VARIABLE}_${1}.txt
+
+r.mapcalc " ${GIS_OPT_VARIABLE}_$1 = ${GIS_OPT_VARIABLE}_$1 " --o
+
+r.out.gdal -c input=${GIS_OPT_VARIABLE}_${1} nodata=-9999 output=$GIS_OPT_OUT_FOLDER/${GIS_OPT_VARIABLE}_${1}.tif type=Int32 -c createopt="COMPRESS=LZW,ZLEVEL=9" --o --q 2> /dev/null
+if test -f "$GIS_OPT_FOLDER/blockfile/${GIS_OPT_VARIABLE}_${1}.tif"; then echo "$GIS_OPT_FOLDER/blockfile/${GIS_OPT_VARIABLE}_${1}.tif has been created!";fi
+#echo "$GIS_OPT_OUT_FOLDER/blockfile/${GIS_OPT_VARIABLE}_${1}.tif has been created!"
+
+' _
+
+exit
+
+echo "The full stream variable calculation has been done!"
+echo "To list the stream-specific output variables: ls $GIS_OPT_OUT_FOLDER/*.tif"
Modified: grass-addons/grass7/raster/r.stream.watersheds/Makefile
===================================================================
--- grass-addons/grass7/raster/r.stream.watersheds/Makefile 2017-05-06 21:52:14 UTC (rev 71042)
+++ grass-addons/grass7/raster/r.stream.watersheds/Makefile 2017-05-06 22:05:54 UTC (rev 71043)
@@ -2,6 +2,6 @@
PGM=r.stream.watersheds
-include $(MODULE_TOPDIR)/include/Make/Script.make
+include $(MODULE_TOPDIR)/include/Make/ShScript.make
default: script
Deleted: grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds
===================================================================
--- grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds 2017-05-06 21:52:14 UTC (rev 71042)
+++ grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds 2017-05-06 22:05:54 UTC (rev 71043)
@@ -1,348 +0,0 @@
-#!/bin/bash
-#
-############################################################################
-#
-# MODULE: r.stream.watersheds
-#
-# AUTHOR(S): Giuseppe Amatulli & Sami Domisch
-# based on "Domisch, S., Amatulli, G., Jetz, W. (in review)
-# Near-global freshwater-specific environmental variables for
-# biodiversity analyses in 1km resolution. Scientific Data"
-#
-# PURPOSE: Delineate the upstream contributing area ('sub-watershed') and
-# stream sections ('sub-stream') for each grid cell of a
-# stream network
-#
-# COPYRIGHT: (C) 2001-2012 by the GRASS Development Team
-#
-# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
-# for details.
-#
-#############################################################################
-
-#%Module
-#% description: Sub-watershed and sub-stream delineation based on the drainage direction and a gridded stream network.
-#% keywords: drainage, stream, hydrology
-#%End
-
-#%option
-#% key: drainage
-#% type: string
-#% key_desc: name
-#% description: Name of the drainage direction raster (generated with r.watershed)
-#% required : yes
-#%end
-
-#%option
-#% key: stream
-#% type: string
-#% key_desc: name
-#% description: Name of the stream network raster
-#% required : yes
-#%end
-
-#%option
-#% key: folder
-#% type: string
-#% key_desc: name
-#% description: Provide the full folder path and name where the sub-watersheds and sub-streams should be stored
-#% required:no
-#% answer: GISDBASE/folder_structure
-#%end
-
-#%option
-#% key: cpu
-#% type: double
-#% description: Number of CPUs used for the parallel computation
-#% answer: 1
-#% required:no
-#%end
-
-if [ -z "$GISBASE" ] ; then
- echo "You must be in GRASS GIS to run this program." >&2
- exit 1
-fi
-
-if [ "$1" != "@ARGS_PARSED@" ] ; then
- exec g.parser "$0" "$@"
-fi
-
-#### check if we have awk
-if [ ! -x "`which awk`" ] ; then
- g.message -e "awk required, please install awk or gawk first"
- exit 1
-fi
-
-# setting environment, so that awk works properly in all languages
-unset LC_ALL
-LC_NUMERIC=C
-export LC_NUMERIC
-
-#test:
-
-if [ -z "$GIS_OPT_DRAINAGE" ] ; then
- g.message "Please provide the name of the drainage direction raster."
-exit 1
-fi
-
-if [ -z "$GIS_OPT_STREAM" ] ; then
- g.message "Please provide the name of the stream network raster."
-exit 1
-fi
-
-#check if drainage and stream file exists
-eval `g.findfile element=cell file="$GIS_OPT_STREAM"`
-if [ ! "$file" ] ; then
- g.message -e "<$GIS_OPT_STREAM> does not exist! Aborting."
- exit 1
-fi
-
-eval `g.findfile element=cell file="$GIS_OPT_DRAINAGE"`
-if [ ! "$file" ] ; then
- g.message -e "<$GIS_OPT_DRAINAGE> does not exist! Aborting."
- exit 1
-fi
-
-export GISDBASE=$(g.gisenv get=GISDBASE separator=space)
-export LOCATION_NAME=$(g.gisenv get=LOCATION_NAME separator=space)
-export GIS_OPT_STREAM
-export GIS_OPT_DRAINAGE
-export GIS_OPT_CPU
-
-if [ ${GIS_OPT_FOLDER} = "GISDBASE/folder_structure" ] ; then
- export GIS_OPT_FOLDER=$GISDBASE"/folder_structure"
-else
- export GIS_OPT_FOLDER=$GIS_OPT_FOLDER
-fi
-
-# what to do in case of user break:
-exitprocedure()
-{
-echo ""
-g.message -e 'Process aborted by user. All intermediate files have been deleted!'
-
-# reset in the permanent mapset
-g.gisenv set="MAPSET=PERMANENT"
-g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
-g.gisenv set="GISDBASE=$GISDBASE"
-
-# delete folder structure and sub-watershed mapset
-rm -fr $GIS_OPT_FOLDER $GISDBASE/$LOCATION_NAME/sub_watershedID*
-
-exit 1
-}
-
-# shell check for user break (signal list: trap -l)
-trap "exitprocedure" 2 3 15
-
-echo ""
-echo "#######################################################################################################"
-echo "Sub-watershed and sub-stream delineation based on the drainage direction and a gridded stream network."
-echo ""
-echo "Citation: Domisch, S., Amatulli, G., Jetz, W. (in review)"
-echo "Near-global freshwater-specific environmental variables for"
-echo "biodiversity analyses in 1km resolution. Scientific Data"
-echo "#######################################################################################################"
-echo ""
-echo Exporting stream grid cell coordinates and saving in $GIS_OPT_FOLDER/stream_coord_lines.txt
-
-rm -fr $GIS_OPT_FOLDER 2> /dev/null
-mkdir $GIS_OPT_FOLDER 2> /dev/null
-
-r.out.xyz --o input=$GIS_OPT_STREAM separator=space output=$GIS_OPT_FOLDER/stream_coord.txt --q
-
-awk 'BEGIN {print "X","Y","ID"} { print $1, $2 , NR }' $GIS_OPT_FOLDER/stream_coord.txt > $GIS_OPT_FOLDER/stream_coord_lines.txt
-
-echo Create the raster ID file
-
-v.in.ascii --overwrite input=$GIS_OPT_FOLDER/stream_coord_lines.txt output=vector_ID x=1 y=2 separator=' ' skip=1 --q 2> /dev/null
-
-echo Rasterize the stream network point coordinates
-
-v.to.rast --overwrite input=vector_ID output=grid_ID layer=vector_ID use=attr attrcolumn=cat type=point --q
-
-echo Create the folder structure in $GIS_OPT_FOLDER/ to save results
-
-# build up folder structure for all the digit bigger than 9999.
-# from 0 to 9999 goes in the /0digit4/0digit3/0digit2/0digit1/
-# from 10000 to 19999 goes in the /0digit4/0digit3/0digit2/1digit1/
-
-max_line=$(wc -l $GIS_OPT_FOLDER/stream_coord.txt | awk '{ print $1 }' )
-max_seq1d=${max_line: -5:1} ; if [ -z $max_seq1d ] ; then max_seq1d=0 ; fi
-max_seq2d=${max_line: -6:1} ; if [ -z $max_seq2d ] ; then max_seq2d=0 ; else max_seq1d=9 ; fi
-max_seq3d=${max_line: -7:1} ; if [ -z $max_seq3d ] ; then max_seq3d=0 ; else max_seq2d=9 ; max_seq1d=9 ; fi
-max_seq4d=${max_line: -8:1} ; if [ -z $max_seq4d ] ; then max_seq4d=0 ; else max_seq3d=9 ; max_seq2d=9 ; max_seq1d=9 ; fi
-
-for dir4d in $(seq 0 $max_seq4d) ; do
- for dir3d in $(seq 0 $max_seq3d) ; do
- for dir2d in $(seq 0 $max_seq2d) ; do
- for dir1d in $(seq 0 $max_seq1d) ; do
- mkdir -p $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1
- done
- done
- done
-done
-
-mkdir $GIS_OPT_FOLDER/blockfile
-
-echo "Split the stream grid cell coordinate file with $( wc -l $GIS_OPT_FOLDER/stream_coord.txt | awk '{ print $1 }' ) grid cells into blocks of 10000 cells"
-
-awk -v PATH=$GIS_OPT_FOLDER 'NR%10000==1{x=PATH"/blockfile/blockfile"(++i*10000-10000)}{print > x}' $GIS_OPT_FOLDER/stream_coord_lines.txt
-
-awk '{ if(NR>1) print }' $GIS_OPT_FOLDER/blockfile/blockfile0 > $GIS_OPT_FOLDER/blockfile/blockfile0_tmp
-mv $GIS_OPT_FOLDER/blockfile/blockfile0_tmp $GIS_OPT_FOLDER/blockfile/blockfile0
-
-echo Create $( ls $GIS_OPT_FOLDER/blockfile/blockfile* | wc -l ) subsets: $GIS_OPT_FOLDER/blockfile/blockfile*.tif
-
-rm -fr $GISDBASE/$LOCATION_NAME/sub_watershedID*
-
-export GISRC_def=$GISRC
-
-for file in $GIS_OPT_FOLDER/blockfile/blockfile* ; do
-
- filename=$(basename $file )
-
- export BLKID=$(echo $filename | awk '{ gsub("0000", "") ; gsub("blockfile","") ; print }')
- export dir4d=${BLKID: -4:1} ; if [ -z $dir4d ] ; then dir4d=0 ; fi
- export dir3d=${BLKID: -3:1} ; if [ -z $dir3d ] ; then dir3d=0 ; fi
- export dir2d=${BLKID: -2:1} ; if [ -z $dir2d ] ; then dir2d=0 ; fi
- export dir1d=${BLKID: -1:1} ; if [ -z $dir1d ] ; then dir1d=0 ; fi
- export file=$file
-
-rm -f $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/*.tif
-
-g.region -d rast=$GIS_OPT_DRAINAGE at PERMANENT
-
-g.gisenv set="GISDBASE=$GISDBASE"
-g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
-g.gisenv set="MAPSET=PERMANENT"
-
-echo "Start the sub-watershed delineation for subset $file"
-
-awk '{ print NR , $1 , $2 ,$3 }' $file | xargs -n 4 -P $GIS_OPT_CPU bash -c $'
-NR=$1
-X_coord=$2
-Y_coord=$3
-ID=$4
-
-
-if [ $( expr $( expr $NR + 1 ) / 100) -eq $( expr $NR / 100 + 1 ) ] ; then
-echo -en "\r$(expr $NR / 100 + 1 )% done"
-fi
-
-
-# replicate the start-GISRC in a unique GISRC
-
-cp $GISRC_def $HOME/.grass7/rc$ID
-export GISRC=$HOME/.grass7/rc$ID
-
-g.mapset -c mapset=sub_watershedID$ID location=$LOCATION_NAME dbase=$GISDBASE --quiet
-
-rm -f $GISDBASE/$LOCATION_NAME/sub_watershedID$ID/.gislock
-
-export GISBASE=$( grep gisbase $(which grass70) | awk \'{ if(NR==2) { gsub ("\\"","" ) ; print $3 } }\' )
-export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
-export LD_LIBRARY_PATH="$GISBASE/lib"
-export GRASS_LD_LIBRARY_PATH="$LD_LIBRARY_PATH"
-export PYTHONPATH="$GISBASE/etc/python:$PYTHONPATH"
-export MANPATH=$MANPATH:$GISBASE/man
-
-r.water.outlet --o --q input=$GIS_OPT_DRAINAGE at PERMANENT output=sub_watershedID$ID coordinates=$X_coord,$Y_coord
-
-g.region rast=sub_watershedID${ID}@sub_watershedID${ID} zoom=sub_watershedID${ID}@sub_watershedID${ID}
-r.out.gdal -c type=Byte --o --q nodata=255 createopt="COMPRESS=LZW,ZLEVEL=9" input=sub_watershedID$ID at sub_watershedID${ID} output=$GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_watershedID$ID.tif
-
-### Crop this basin template --> for temperature only for the stream cells in this basin
-
-r.mapcalc "sub_streamID${ID} = if ( sub_watershedID${ID}@sub_watershedID${ID} == 1 & ${GIS_OPT_STREAM}@PERMANENT >= 1 , 1 , null() )" --o --q
-r.null map=sub_streamID${ID} setnull=0 --q
-r.out.gdal -c type=Byte --o --q nodata=255 createopt="COMPRESS=LZW,ZLEVEL=9" input=sub_streamID${ID}@sub_watershedID${ID} output=$GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_streamID${ID}.tif
-
-if [ ! -f $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_streamID${ID}.tif ] ; then the file $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_streamID${ID}.tif dose not exist ; fi
-
-rm -r $HOME/.grass7/rc$ID $GISDBASE/$LOCATION_NAME/sub_watershedID$ID
-
-' _
-
-echo -en "\r100% done"
-
-rm -fr $GISDBASE/$LOCATION_NAME/sub_watershedID*
-
-echo ""
-echo $(ls -l $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_watershedID*.tif | wc -l ) sub-watersheds have been processed
-
-# reset mapset to PERMANENT
-
-g.gisenv set="MAPSET=PERMANENT"
-g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
-g.gisenv set="GISDBASE=$GISDBASE"
-
-echo Check for missing files due to RAM overload
-
-cat $file | xargs -n 3 -P 1 bash -c $'
-
-X_coord=$1
-Y_coord=$2
-ID=$3
-
-if [ ! -f $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_watershedID$ID.tif ] ; then
-
-echo Fix missing file $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_watershedID$ID.tif using only 1 CPU
-
-# replicate the start-GISRC in a unique GISRC
-
-cp $GISRC_def $HOME/.grass7/rc$ID
-export GISRC=$HOME/.grass7/rc$ID
-
-g.mapset -c mapset=sub_watershedID$ID location=$LOCATION_NAME dbase=$GISDBASE --quiet
-
-rm -f $GISDBASE/$LOCATION_NAME/sub_watershedID$ID/.gislock
-
-export GISBASE=$( grep gisbase $(which grass70) | awk \'{ if(NR==2) { gsub ("\\"","" ) ; print $3 } }\' )
-export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
-export LD_LIBRARY_PATH="$GISBASE/lib"
-export GRASS_LD_LIBRARY_PATH="$LD_LIBRARY_PATH"
-export PYTHONPATH="$GISBASE/etc/python:$PYTHONPATH"
-export MANPATH=$MANPATH:$GISBASE/man
-
-r.water.outlet --o --q input=$GIS_OPT_DRAINAGE at PERMANENT output=sub_watershedID$ID coordinates=$X_coord,$Y_coord
-
-g.region rast=sub_watershedID${ID}@sub_watershedID${ID} zoom=sub_watershedID${ID}@sub_watershedID${ID}
-r.out.gdal -c type=Byte --o --q nodata=255 createopt="COMPRESS=LZW,ZLEVEL=9" input=sub_watershedID$ID at sub_watershedID${ID} output=$GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_watershedID$ID.tif
-
-### Crop this basin template --> for temperature only for the stream cells in this basin
-
-r.mapcalc "sub_streamID${ID} = if ( sub_watershedID${ID}@sub_watershedID${ID} == 1 & ${GIS_OPT_STREAM}@PERMANENT >= 1 , 1 , null() )" --o --q
-r.null map=sub_streamID${ID} setnull=0 --q
-r.out.gdal -c type=Byte --o --q nodata=255 createopt="COMPRESS=LZW,ZLEVEL=9" input=sub_streamID${ID}@sub_watershedID${ID} output=$GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_streamID${ID}.tif
-
-rm -r $HOME/.grass7/rc$ID $GISDBASE/$LOCATION_NAME/sub_watershedID$ID
-
-fi
-
-' _
-
-rm -fr $GISDBASE/$LOCATION_NAME/sub_watershedID*
-
-echo "Compress the sub-watersheds and sub-streams to reduce the number of inodes filling up the hard disk"
-
-
-cd $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/
-
-tar -zcPf ${filename}_sub_watershed.tar.gz sub_watershedID*.tif
-tar -zcPf ${filename}_sub_stream.tar.gz sub_streamID*.tif
-
-cd $GIS_OPT_FOLDER
-
-g.gisenv set="MAPSET=PERMANENT"
-g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
-g.gisenv set="GISDBASE=$GISDBASE"
-
-rm -fr $GISDBASE/$LOCATION_NAME/sub_watershedID* $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/*.tif
-
-done
-
-echo "The full sub-watershed delineation process has been done!"
-echo "To list the compressed files: ls $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/*.tar.gz"
-echo "Now you can use r.stream.variables to compute stream-specific environmental variables."
Copied: grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds.sh (from rev 71042, grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds)
===================================================================
--- grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds.sh (rev 0)
+++ grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds.sh 2017-05-06 22:05:54 UTC (rev 71043)
@@ -0,0 +1,348 @@
+#!/bin/bash
+#
+############################################################################
+#
+# MODULE: r.stream.watersheds
+#
+# AUTHOR(S): Giuseppe Amatulli & Sami Domisch
+# based on "Domisch, S., Amatulli, G., Jetz, W. (in review)
+# Near-global freshwater-specific environmental variables for
+# biodiversity analyses in 1km resolution. Scientific Data"
+#
+# PURPOSE: Delineate the upstream contributing area ('sub-watershed') and
+# stream sections ('sub-stream') for each grid cell of a
+# stream network
+#
+# COPYRIGHT: (C) 2001-2012 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+#%Module
+#% description: Sub-watershed and sub-stream delineation based on the drainage direction and a gridded stream network.
+#% keywords: drainage, stream, hydrology
+#%End
+
+#%option
+#% key: drainage
+#% type: string
+#% key_desc: name
+#% description: Name of the drainage direction raster (generated with r.watershed)
+#% required : yes
+#%end
+
+#%option
+#% key: stream
+#% type: string
+#% key_desc: name
+#% description: Name of the stream network raster
+#% required : yes
+#%end
+
+#%option
+#% key: folder
+#% type: string
+#% key_desc: name
+#% description: Provide the full folder path and name where the sub-watersheds and sub-streams should be stored
+#% required:no
+#% answer: GISDBASE/folder_structure
+#%end
+
+#%option
+#% key: cpu
+#% type: double
+#% description: Number of CPUs used for the parallel computation
+#% answer: 1
+#% required:no
+#%end
+
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." >&2
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+#### check if we have awk
+if [ ! -x "`which awk`" ] ; then
+ g.message -e "awk required, please install awk or gawk first"
+ exit 1
+fi
+
+# setting environment, so that awk works properly in all languages
+unset LC_ALL
+LC_NUMERIC=C
+export LC_NUMERIC
+
+#test:
+
+if [ -z "$GIS_OPT_DRAINAGE" ] ; then
+ g.message "Please provide the name of the drainage direction raster."
+exit 1
+fi
+
+if [ -z "$GIS_OPT_STREAM" ] ; then
+ g.message "Please provide the name of the stream network raster."
+exit 1
+fi
+
+#check if drainage and stream file exists
+eval `g.findfile element=cell file="$GIS_OPT_STREAM"`
+if [ ! "$file" ] ; then
+ g.message -e "<$GIS_OPT_STREAM> does not exist! Aborting."
+ exit 1
+fi
+
+eval `g.findfile element=cell file="$GIS_OPT_DRAINAGE"`
+if [ ! "$file" ] ; then
+ g.message -e "<$GIS_OPT_DRAINAGE> does not exist! Aborting."
+ exit 1
+fi
+
+export GISDBASE=$(g.gisenv get=GISDBASE separator=space)
+export LOCATION_NAME=$(g.gisenv get=LOCATION_NAME separator=space)
+export GIS_OPT_STREAM
+export GIS_OPT_DRAINAGE
+export GIS_OPT_CPU
+
+if [ ${GIS_OPT_FOLDER} = "GISDBASE/folder_structure" ] ; then
+ export GIS_OPT_FOLDER=$GISDBASE"/folder_structure"
+else
+ export GIS_OPT_FOLDER=$GIS_OPT_FOLDER
+fi
+
+# what to do in case of user break:
+exitprocedure()
+{
+echo ""
+g.message -e 'Process aborted by user. All intermediate files have been deleted!'
+
+# reset in the permanent mapset
+g.gisenv set="MAPSET=PERMANENT"
+g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
+g.gisenv set="GISDBASE=$GISDBASE"
+
+# delete folder structure and sub-watershed mapset
+rm -fr $GIS_OPT_FOLDER $GISDBASE/$LOCATION_NAME/sub_watershedID*
+
+exit 1
+}
+
+# shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 15
+
+echo ""
+echo "#######################################################################################################"
+echo "Sub-watershed and sub-stream delineation based on the drainage direction and a gridded stream network."
+echo ""
+echo "Citation: Domisch, S., Amatulli, G., Jetz, W. (in review)"
+echo "Near-global freshwater-specific environmental variables for"
+echo "biodiversity analyses in 1km resolution. Scientific Data"
+echo "#######################################################################################################"
+echo ""
+echo Exporting stream grid cell coordinates and saving in $GIS_OPT_FOLDER/stream_coord_lines.txt
+
+rm -fr $GIS_OPT_FOLDER 2> /dev/null
+mkdir $GIS_OPT_FOLDER 2> /dev/null
+
+r.out.xyz --o input=$GIS_OPT_STREAM separator=space output=$GIS_OPT_FOLDER/stream_coord.txt --q
+
+awk 'BEGIN {print "X","Y","ID"} { print $1, $2 , NR }' $GIS_OPT_FOLDER/stream_coord.txt > $GIS_OPT_FOLDER/stream_coord_lines.txt
+
+echo Create the raster ID file
+
+v.in.ascii --overwrite input=$GIS_OPT_FOLDER/stream_coord_lines.txt output=vector_ID x=1 y=2 separator=' ' skip=1 --q 2> /dev/null
+
+echo Rasterize the stream network point coordinates
+
+v.to.rast --overwrite input=vector_ID output=grid_ID layer=vector_ID use=attr attrcolumn=cat type=point --q
+
+echo Create the folder structure in $GIS_OPT_FOLDER/ to save results
+
+# build up folder structure for all the digit bigger than 9999.
+# from 0 to 9999 goes in the /0digit4/0digit3/0digit2/0digit1/
+# from 10000 to 19999 goes in the /0digit4/0digit3/0digit2/1digit1/
+
+max_line=$(wc -l $GIS_OPT_FOLDER/stream_coord.txt | awk '{ print $1 }' )
+max_seq1d=${max_line: -5:1} ; if [ -z $max_seq1d ] ; then max_seq1d=0 ; fi
+max_seq2d=${max_line: -6:1} ; if [ -z $max_seq2d ] ; then max_seq2d=0 ; else max_seq1d=9 ; fi
+max_seq3d=${max_line: -7:1} ; if [ -z $max_seq3d ] ; then max_seq3d=0 ; else max_seq2d=9 ; max_seq1d=9 ; fi
+max_seq4d=${max_line: -8:1} ; if [ -z $max_seq4d ] ; then max_seq4d=0 ; else max_seq3d=9 ; max_seq2d=9 ; max_seq1d=9 ; fi
+
+for dir4d in $(seq 0 $max_seq4d) ; do
+ for dir3d in $(seq 0 $max_seq3d) ; do
+ for dir2d in $(seq 0 $max_seq2d) ; do
+ for dir1d in $(seq 0 $max_seq1d) ; do
+ mkdir -p $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1
+ done
+ done
+ done
+done
+
+mkdir $GIS_OPT_FOLDER/blockfile
+
+echo "Split the stream grid cell coordinate file with $( wc -l $GIS_OPT_FOLDER/stream_coord.txt | awk '{ print $1 }' ) grid cells into blocks of 10000 cells"
+
+awk -v PATH=$GIS_OPT_FOLDER 'NR%10000==1{x=PATH"/blockfile/blockfile"(++i*10000-10000)}{print > x}' $GIS_OPT_FOLDER/stream_coord_lines.txt
+
+awk '{ if(NR>1) print }' $GIS_OPT_FOLDER/blockfile/blockfile0 > $GIS_OPT_FOLDER/blockfile/blockfile0_tmp
+mv $GIS_OPT_FOLDER/blockfile/blockfile0_tmp $GIS_OPT_FOLDER/blockfile/blockfile0
+
+echo Create $( ls $GIS_OPT_FOLDER/blockfile/blockfile* | wc -l ) subsets: $GIS_OPT_FOLDER/blockfile/blockfile*.tif
+
+rm -fr $GISDBASE/$LOCATION_NAME/sub_watershedID*
+
+export GISRC_def=$GISRC
+
+for file in $GIS_OPT_FOLDER/blockfile/blockfile* ; do
+
+ filename=$(basename $file )
+
+ export BLKID=$(echo $filename | awk '{ gsub("0000", "") ; gsub("blockfile","") ; print }')
+ export dir4d=${BLKID: -4:1} ; if [ -z $dir4d ] ; then dir4d=0 ; fi
+ export dir3d=${BLKID: -3:1} ; if [ -z $dir3d ] ; then dir3d=0 ; fi
+ export dir2d=${BLKID: -2:1} ; if [ -z $dir2d ] ; then dir2d=0 ; fi
+ export dir1d=${BLKID: -1:1} ; if [ -z $dir1d ] ; then dir1d=0 ; fi
+ export file=$file
+
+rm -f $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/*.tif
+
+g.region -d rast=$GIS_OPT_DRAINAGE at PERMANENT
+
+g.gisenv set="GISDBASE=$GISDBASE"
+g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
+g.gisenv set="MAPSET=PERMANENT"
+
+echo "Start the sub-watershed delineation for subset $file"
+
+awk '{ print NR , $1 , $2 ,$3 }' $file | xargs -n 4 -P $GIS_OPT_CPU bash -c $'
+NR=$1
+X_coord=$2
+Y_coord=$3
+ID=$4
+
+
+if [ $( expr $( expr $NR + 1 ) / 100) -eq $( expr $NR / 100 + 1 ) ] ; then
+echo -en "\r$(expr $NR / 100 + 1 )% done"
+fi
+
+
+# replicate the start-GISRC in a unique GISRC
+
+cp $GISRC_def $HOME/.grass7/rc$ID
+export GISRC=$HOME/.grass7/rc$ID
+
+g.mapset -c mapset=sub_watershedID$ID location=$LOCATION_NAME dbase=$GISDBASE --quiet
+
+rm -f $GISDBASE/$LOCATION_NAME/sub_watershedID$ID/.gislock
+
+export GISBASE=$( grep gisbase $(which grass70) | awk \'{ if(NR==2) { gsub ("\\"","" ) ; print $3 } }\' )
+export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
+export LD_LIBRARY_PATH="$GISBASE/lib"
+export GRASS_LD_LIBRARY_PATH="$LD_LIBRARY_PATH"
+export PYTHONPATH="$GISBASE/etc/python:$PYTHONPATH"
+export MANPATH=$MANPATH:$GISBASE/man
+
+r.water.outlet --o --q input=$GIS_OPT_DRAINAGE at PERMANENT output=sub_watershedID$ID coordinates=$X_coord,$Y_coord
+
+g.region rast=sub_watershedID${ID}@sub_watershedID${ID} zoom=sub_watershedID${ID}@sub_watershedID${ID}
+r.out.gdal -c type=Byte --o --q nodata=255 createopt="COMPRESS=LZW,ZLEVEL=9" input=sub_watershedID$ID at sub_watershedID${ID} output=$GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_watershedID$ID.tif
+
+### Crop this basin template --> for temperature only for the stream cells in this basin
+
+r.mapcalc "sub_streamID${ID} = if ( sub_watershedID${ID}@sub_watershedID${ID} == 1 & ${GIS_OPT_STREAM}@PERMANENT >= 1 , 1 , null() )" --o --q
+r.null map=sub_streamID${ID} setnull=0 --q
+r.out.gdal -c type=Byte --o --q nodata=255 createopt="COMPRESS=LZW,ZLEVEL=9" input=sub_streamID${ID}@sub_watershedID${ID} output=$GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_streamID${ID}.tif
+
+if [ ! -f $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_streamID${ID}.tif ] ; then the file $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_streamID${ID}.tif dose not exist ; fi
+
+rm -r $HOME/.grass7/rc$ID $GISDBASE/$LOCATION_NAME/sub_watershedID$ID
+
+' _
+
+echo -en "\r100% done"
+
+rm -fr $GISDBASE/$LOCATION_NAME/sub_watershedID*
+
+echo ""
+echo $(ls -l $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_watershedID*.tif | wc -l ) sub-watersheds have been processed
+
+# reset mapset to PERMANENT
+
+g.gisenv set="MAPSET=PERMANENT"
+g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
+g.gisenv set="GISDBASE=$GISDBASE"
+
+echo Check for missing files due to RAM overload
+
+cat $file | xargs -n 3 -P 1 bash -c $'
+
+X_coord=$1
+Y_coord=$2
+ID=$3
+
+if [ ! -f $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_watershedID$ID.tif ] ; then
+
+echo Fix missing file $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_watershedID$ID.tif using only 1 CPU
+
+# replicate the start-GISRC in a unique GISRC
+
+cp $GISRC_def $HOME/.grass7/rc$ID
+export GISRC=$HOME/.grass7/rc$ID
+
+g.mapset -c mapset=sub_watershedID$ID location=$LOCATION_NAME dbase=$GISDBASE --quiet
+
+rm -f $GISDBASE/$LOCATION_NAME/sub_watershedID$ID/.gislock
+
+export GISBASE=$( grep gisbase $(which grass70) | awk \'{ if(NR==2) { gsub ("\\"","" ) ; print $3 } }\' )
+export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
+export LD_LIBRARY_PATH="$GISBASE/lib"
+export GRASS_LD_LIBRARY_PATH="$LD_LIBRARY_PATH"
+export PYTHONPATH="$GISBASE/etc/python:$PYTHONPATH"
+export MANPATH=$MANPATH:$GISBASE/man
+
+r.water.outlet --o --q input=$GIS_OPT_DRAINAGE at PERMANENT output=sub_watershedID$ID coordinates=$X_coord,$Y_coord
+
+g.region rast=sub_watershedID${ID}@sub_watershedID${ID} zoom=sub_watershedID${ID}@sub_watershedID${ID}
+r.out.gdal -c type=Byte --o --q nodata=255 createopt="COMPRESS=LZW,ZLEVEL=9" input=sub_watershedID$ID at sub_watershedID${ID} output=$GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_watershedID$ID.tif
+
+### Crop this basin template --> for temperature only for the stream cells in this basin
+
+r.mapcalc "sub_streamID${ID} = if ( sub_watershedID${ID}@sub_watershedID${ID} == 1 & ${GIS_OPT_STREAM}@PERMANENT >= 1 , 1 , null() )" --o --q
+r.null map=sub_streamID${ID} setnull=0 --q
+r.out.gdal -c type=Byte --o --q nodata=255 createopt="COMPRESS=LZW,ZLEVEL=9" input=sub_streamID${ID}@sub_watershedID${ID} output=$GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/sub_streamID${ID}.tif
+
+rm -r $HOME/.grass7/rc$ID $GISDBASE/$LOCATION_NAME/sub_watershedID$ID
+
+fi
+
+' _
+
+rm -fr $GISDBASE/$LOCATION_NAME/sub_watershedID*
+
+echo "Compress the sub-watersheds and sub-streams to reduce the number of inodes filling up the hard disk"
+
+
+cd $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/
+
+tar -zcPf ${filename}_sub_watershed.tar.gz sub_watershedID*.tif
+tar -zcPf ${filename}_sub_stream.tar.gz sub_streamID*.tif
+
+cd $GIS_OPT_FOLDER
+
+g.gisenv set="MAPSET=PERMANENT"
+g.gisenv set="LOCATION_NAME=$LOCATION_NAME"
+g.gisenv set="GISDBASE=$GISDBASE"
+
+rm -fr $GISDBASE/$LOCATION_NAME/sub_watershedID* $GIS_OPT_FOLDER/${dir4d}digit4/${dir3d}digit3/${dir2d}digit2/${dir1d}digit1/*.tif
+
+done
+
+echo "The full sub-watershed delineation process has been done!"
+echo "To list the compressed files: ls $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/*.tar.gz"
+echo "Now you can use r.stream.variables to compute stream-specific environmental variables."
More information about the grass-commit
mailing list