[GRASS-SVN] r66293 - in grass-addons/grass7/raster: . r.stream.variables r.stream.watersheds
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Sep 22 07:58:12 PDT 2015
Author: elselvaje
Date: 2015-09-22 07:58:12 -0700 (Tue, 22 Sep 2015)
New Revision: 66293
Added:
grass-addons/grass7/raster/r.stream.variables/
grass-addons/grass7/raster/r.stream.variables/Makefile
grass-addons/grass7/raster/r.stream.variables/r.stream.variables
grass-addons/grass7/raster/r.stream.watersheds/
grass-addons/grass7/raster/r.stream.watersheds/Makefile
grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds
Log:
initial commit
Added: grass-addons/grass7/raster/r.stream.variables/Makefile
===================================================================
--- grass-addons/grass7/raster/r.stream.variables/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.stream.variables/Makefile 2015-09-22 14:58:12 UTC (rev 66293)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM=r.stream.variables
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.stream.variables/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.stream.variables/r.stream.variables
===================================================================
--- grass-addons/grass7/raster/r.stream.variables/r.stream.variables (rev 0)
+++ grass-addons/grass7/raster/r.stream.variables/r.stream.variables 2015-09-22 14:58:12 UTC (rev 66293)
@@ -0,0 +1,365 @@
+#!/bin/bash
+#
+############################################################################
+#
+# MODULE: r.stream.watersheds
+#
+# AUTHOR(S): Giuseppe Amatulli & Sami Domisch
+# based on "Domisch, S., Amatulli, G., Jetz, W. (2015)
+# 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 delineation based on the drainage direction and a gridded stream network. Based on "Domisch, S., Amatulli, G., Jetz, W. 2015 Near-global freshwater-specific environmental variables for biodiversity analyses in 1km resolution. Scientific Data"
+#% 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:yes
+#%end
+
+#%option
+#% key: output
+#% type: string
+#% key_desc: method
+#% multiple: yes
+#% options: min,max,range,mean,stddev,coeff_var,sum
+#% description: Provide the output aggregation method for the stream-specific variable: upstream 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
+#% description: Provide a scale factor to multiply or divide the final stream-specific variable (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
+
+# 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}*
+rm -fr $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/*.txt
+rm -fr $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/*.tif
+rm -fr $GIS_OPT_FOLDER/blockfile/stat_*.txt $GIS_OPT_FOLDER/blockfile/*.tif
+
+exit 1
+}
+
+# shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 9 15 19
+
+echo ""
+echo "Streams variables calculation based on the sub-watershads and sub-streams."
+echo ""
+echo "Citation: Domisch, S., Amatulli, G., Jetz, W. (2015)"
+echo "Near-global freshwater-specific environmental variables for"
+echo "biodiversity analyses in 1km resolution. Scientific Data"
+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
+
+ 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 )
+
+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
+
+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 = "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 = "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.out.gdal -c input=${GIS_OPT_VARIABLE}_${1} nodata=-9999 output=$GIS_OPT_FOLDER/blockfile/${GIS_OPT_VARIABLE}_${1}.tif type=Int32 -c createopt="COMPRESS=LZW,ZLEVEL=9" --o --q 2> /dev/null
+
+echo "$GIS_OPT_FOLDER/blockfile/${GIS_OPT_VARIABLE}_${1}.tif has been created!"
+
+' _
+
+echo "The full stream variable calculation has been done!"
+echo "To list the stream-specific output variables, use ls $GIS_OPT_FOLDER/blockfile/*.tif"
Property changes on: grass-addons/grass7/raster/r.stream.variables/r.stream.variables
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/x-sh
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.stream.watersheds/Makefile
===================================================================
--- grass-addons/grass7/raster/r.stream.watersheds/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.stream.watersheds/Makefile 2015-09-22 14:58:12 UTC (rev 66293)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM=r.stream.watersheds
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.stream.watersheds/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds
===================================================================
--- grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds (rev 0)
+++ grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds 2015-09-22 14:58:12 UTC (rev 66293)
@@ -0,0 +1,349 @@
+#!/bin/bash
+#
+############################################################################
+#
+# MODULE: r.stream.watersheds
+#
+# AUTHOR(S): Giuseppe Amatulli & Sami Domisch
+# based on "Domisch, S., Amatulli, G., Jetz, W. (2015)
+# 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 delineation based on the drainage direction and a gridded stream network. Based on "Domisch, S., Amatulli, G., Jetz, W. (2015) Near-global freshwater-specific environmental variables for biodiversity analyses in 1km resolution. Scientific Data"
+#% 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 will 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 "Sub-watershed and sub-streams delineation based on the drainage direction and a gridded stream network."
+echo ""
+echo "Citation: Domisch, S., Amatulli, G., Jetz, W. (2015)"
+echo "Near-global freshwater-specific environmental variables for"
+echo "biodiversity analyses in 1km resolution. Scientific Data"
+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=int_1 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 ) blocks in $GIS_OPT_FOLDER/blockfile/
+
+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 block $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 uniq 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
+
+cd $GIS_OPT_FOLDER
+
+echo The full sub-watershed delineation process has been done!
+echo To list the compressed files, use "ls $GIS_OPT_FOLDER/*digit4/*digit3/*digit2/*digit1/*.tar.gz"
+echo Now you can use r.stream.variables to compute stream-specific environmental variables.
+
Property changes on: grass-addons/grass7/raster/r.stream.watersheds/r.stream.watersheds
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/x-sh
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list