[GRASS-SVN] r52030 - in grass-addons/grass6/vector/v.strahler: . r.broscoe r.strahler
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jun 11 01:33:17 PDT 2012
Author: neteler
Date: 2012-06-11 01:33:16 -0700 (Mon, 11 Jun 2012)
New Revision: 52030
Added:
grass-addons/grass6/vector/v.strahler/r.broscoe/
grass-addons/grass6/vector/v.strahler/r.broscoe/Makefile
grass-addons/grass6/vector/v.strahler/r.broscoe/description.html
grass-addons/grass6/vector/v.strahler/r.broscoe/r.broscoe
grass-addons/grass6/vector/v.strahler/r.strahler/
grass-addons/grass6/vector/v.strahler/r.strahler/Makefile
grass-addons/grass6/vector/v.strahler/r.strahler/description.html
grass-addons/grass6/vector/v.strahler/r.strahler/r.strahler
grass-addons/grass6/vector/v.strahler/v.mainchannel.sh
Removed:
grass-addons/grass6/vector/v.strahler/r.broscoe.sh
grass-addons/grass6/vector/v.strahler/r.broscoe.sh.html
grass-addons/grass6/vector/v.strahler/r.strahler.sh
grass-addons/grass6/vector/v.strahler/r.strahler.sh.html
grass-addons/grass6/vector/v.strahler/v.mainchannel
Log:
cleanup part 1
Added: grass-addons/grass6/vector/v.strahler/r.broscoe/Makefile
===================================================================
--- grass-addons/grass6/vector/v.strahler/r.broscoe/Makefile (rev 0)
+++ grass-addons/grass6/vector/v.strahler/r.broscoe/Makefile 2012-06-11 08:33:16 UTC (rev 52030)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.broscoe
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Copied: grass-addons/grass6/vector/v.strahler/r.broscoe/description.html (from rev 51972, grass-addons/grass6/vector/v.strahler/r.broscoe.sh.html)
===================================================================
--- grass-addons/grass6/vector/v.strahler/r.broscoe/description.html (rev 0)
+++ grass-addons/grass6/vector/v.strahler/r.broscoe/description.html 2012-06-11 08:33:16 UTC (rev 52030)
@@ -0,0 +1,92 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.broscoe.sh</em> Calculates waerden test and t test statistics for some values of threshold area on a single basin, according to A.J.Broscoe theory (1959).<br />
+The program uses some <em><a href="http://www.r-project.org/">R</a></em> commands for statistical analisys and graphic rapresentation. In particoular the R package <em>"agricolae"</em> is required.<br />
+The A.J.Broscoe theory is well known as the theory of the "Mean Stream Drop" and it says that, for the extraction by DEM of a stream network, exists a threshold value wich makes <em>drop</em> constant, and this is the <em>right</em> one extraction threshold. <br />
+By definig the <em>drop</em> (H) as:<br />
+<br />
+<em>H = S L</em><br />
+<br />
+for streams in each Strahler order. <br />
+The <em>drop</em> is calculated as the quote difference by ending and starting point of a stream; S is the slope and L is the lenght of the same stream.<br />
+Using the Leopold and Miller relation (1964):<br />
+<br />
+<em>S=CA<sup>t</sup></em><br />
+<br />
+where C and t are constat values; the area we are searching (A) is the lowest that gives S to find H costant for each Strahler order (w).<br />
+<br />
+<em>H<sub>w</sub> = H<sub>w+1</sub> = H<sub>w+2</sub> = ...</em><br />
+<br />
+where H<sub>w</sub> is the <em>mean</em> of the drops related to the streams in the same Strahler order (w).<br />
+The area can be found by making some attempts for different area thresholds, doing some statistical tests (Van der Waerden test and linear regression), and choosing the <em>right</em> threshold from the output of the tests.<br />
+<br />
+<em>r.broscoe.sh</em> takes in input the DEM, the threshold values on wich calculate statistics, the outlet coords of the basin you want to study; it returns a table (text file) with the output of the Van der Waerden test and linear regression (t test) for each threshold value.<br />
+For the Van der Waerden test the parameter <em>Pvalue</em> is taken. It has to be greater than the possible, it represents the possibility of success of the test (the <em>Mean Stream Drop</em> is the same for all Strahler orders).<br />
+For the linear regression the parameters <em>t, Pr, R_squared_adj</em> are taken. <em>t</em> is the t statistic value, <em>Pr</em> is the possibility of success of the t test, <em>R_squared_adj</em> measures the dispersion of data around the mean value (for each order) for given degrees of freedom.<br />
+Three graphics called "linear_regression", "waerden_test" and "all_tests" are also generated as PDF in the home folder.<br />
+<br />
+Preferably let's take the threshold value wich gives <em>Pvalue</em> (or <em>Pr</em>) greater than 0.95, but is not granted that you can reach that result because it depends of the well-graduation (by Horton-Strahler) of the basin, its geomorphological maturity, so it is not rare that you have to take threshold where <em>Pvalue</em> is simlpy the greatest.<br />
+At the end of the calculation, at first <em>Pvalue</em> is examinated, then, only if Van der Waerden test gives no good results (low <em>Pvalue</em>), the linear regression output (<em>Pr</em>) is examinated; in fact the Van der Waerden test is preferred to linear regression because it allows you to consider the real dispersion of data around the mean: this makes you able to know the real significance of the probability (e.g. the significance is low for few data in the sample) considering an unique parameter.<br />
+<br />
+<h2>EXAMPLE</h2>
+
+An example on Menotre stream (Umbria, Italy):<br />
+The syntax:
+
+<div class="code"><pre>
+ r.broscoe.sh dem=dtm20_regione at AB 'thresholds=400 600 800 1000 1200 1400 1800 2000' xcoor=2291350.34 ycoor=4765192.22 lt=4 result=menotre_txt
+</pre></div>
+
+The results:
+<div class="code"><pre>
+threshold t Pr Radj Pvalue
+400 0.5713518 0.568486 -0.003798402 0.6085511
+600 0.8791352 0.3810997 -0.001896266 0.2798474
+800 1.053110 0.2948033 0.001067895 0.29454
+1000 0.02578308 0.9794938 -0.01233737 0.8535388
+1200 0.3985548 0.69147 -0.01234108 0.6340721
+1400 -1.024254 0.3100425 0.0008457844 0.256408
+1800 -0.6368832 0.5274277 -0.01309044 0.5764749
+2000 -0.4003206 0.6908575 -0.01901582 0.814699
+</pre></div>
+<br />
+<img src="images/wt_rbroscoe.jpg"> <img src="images/lr_rbroscoe.jpg"> <img src="images/at_rbroscoe.jpg"> <br />
+<br />
+By the report and graphics, you can see that the Van der Werden test gives not-so-good results (<em>Pvalue_max</em>=0.85 for threshold=1000 cells) but, if you consider the linear regression output (<em>Pr</em>), you can see that for the same threshold value (1000 cells) <em>Pr</em> is 97%.<br />
+So the threshold=1000 cells is chosen. Moreover the program returns a set of vector map called <em>"orderd_thresholdvalue"</em> from wich you can extract the right one orderd-network (in this case the right one is <em>"orderd_1000"</em>), you can rename and use it as well as you want.<br />
+<br />
+<img src="images/menotre.jpg"><br />
+<br />
+<h2>NOTES</h2>
+The <em>lt</em> value requested in input is a parameter that prevents eventual errors in the DEM; it considers the presence of pits and represents the height difference <em>lesserthan</em> a drop is not considered as a drop but as a pit, and extracted from <em>Mean Stream Drop</em> analysis.<br />
+<br />
+The program uses the module <em><a href="r.strahler.html">r.strahler</a></em>, so it presents the same conditions about the selection of a threshold value range.<br />
+
+<h2>SEE ALSO</h2>
+<em><a href="v.strahler.html">v.strahler</a></em><br>
+<em><a href="r.strahler.sh.html">r.strahler.sh</a></em><br>
+
+<h2>REFERENCES</h2>
+NIST, (2006). <i>Van Der Waerden.</i><br />
+URL: <em><a href="http://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/vanderwa.htm">http://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/vanderwa.htm</a></em><br />
+<p>
+D. G. Tarboton and D. P. Ames, (2001). <i>Advances in the mapping of flow networks from digital elevation data.</i><b> World Water and Environment Resources Congress</b>, presentation (2001).<br />
+<p>
+J. J. Flint, (1974). <i>Stream gradient as a function of order, magnitude, and discharge.</i><b> Water Resources Research</b>, vol.10, n.5, p.969-973. <br />
+<p>
+NIST, (2006). <i>Engineering statistical handbook: confidence limits for the mean.</i><br />
+URL: <em><a href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm">http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm</a></em><br />
+<p>
+J. C. Davis, (1990). <i>Statistics and Data Analysis in Geology</i>. John Wiley \& Sons editors (New York, NY, USA).<br />
+<p>
+A. J. Broscoe, (1959). <i>Quantitative analysis of longitudinal stream profiles of small watersheds</i>. Department of Geology, Columbia University, NY.<br />
+<p>
+F. De Mendiburu, (2006). <i>Statistical Procedures for Agricultural Research.</i><br />
+URL: <em><a href="http://rss.acs.unt.edu/Rdoc/library/agricolae/html/agricolae.package.html">http://rss.acs.unt.edu/Rdoc/library/agricolae/html/agricolae.package.html</a></em><br />
+
+<h2>AUTHORS</h2>
+
+Ivan Marchesini and Annalisa Minelli, Univ. Perugia. <br>
+
+<p>
+<i>Last changed: $Date$</i>
Copied: grass-addons/grass6/vector/v.strahler/r.broscoe/r.broscoe (from rev 51972, grass-addons/grass6/vector/v.strahler/r.broscoe.sh)
===================================================================
--- grass-addons/grass6/vector/v.strahler/r.broscoe/r.broscoe (rev 0)
+++ grass-addons/grass6/vector/v.strahler/r.broscoe/r.broscoe 2012-06-11 08:33:16 UTC (rev 52030)
@@ -0,0 +1,303 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: r.broscoe
+# AUTHOR(S): Annalisa Minelli, Ivan Marchesini
+#
+# PURPOSE: Calculates waerden test and t test statistics
+# for some values of threshold, on a single basin
+# according to A. J. Broscoe theory (1959)
+#
+# COPYRIGHT: (C) 2008 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.
+#
+# REQUIREMENTS: you need R installed with "agricolae" package
+# You can install the package by starting R as root and typing
+# install.packages("agricolae")
+# This will install the package from a CRAN mirror on the internet.
+# For more information on R see the documentation available at
+# http://www.r-project.org/
+#
+# TODO: solve stability problems for low area threshold (cf. v.strahler)
+#############################################################################
+#%Module
+#% description: Calculates waerden test and t test statistics for some values of threshold, on a single basin according to A. J. Broscoe theory (1959)
+#% keywords: waerden.test,t.test,mean stream drop
+#%End
+#%option
+#% key: dem
+#% type: string
+#% key_desc: dem
+#% gisprompt: old,cell,raster
+#% description: Name of DEM raster map
+#% required : yes
+#%END
+#%option
+#% key: thresholds
+#% type: integer
+#% description: Threshold values to calculate statistics on (separated by <space>)
+#% required : yes
+#%end
+#%option
+#% key: xcoor
+#% type: double
+#% description: x coord of outlet
+#% required : yes
+#%end
+#%option
+#% key: ycoor
+#% type: double
+#% description: y coord of outlet
+#% required : yes
+#%end
+#%option
+#% key: lt
+#% type: integer
+#% description: Lesser than (in meters), the program doesn't consider stream drops lesser than this
+#% required : yes
+#%END
+#%option
+#% key: result
+#% type: string
+#% gisprompt: new_file,file,output
+#% key_desc: name
+#% description: Resultant text file to write statistics into
+#% required: yes
+#%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
+
+dem=$GIS_OPT_DEM
+thresholds=$GIS_OPT_THRESHOLDS
+xcoor=$GIS_OPT_XCOOR
+ycoor=$GIS_OPT_YCOOR
+lt=$GIS_OPT_LT
+result=$GIS_OPT_RESULT
+
+### setup enviro vars ###
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+
+echo $LOCATION_NAME
+
+#g.region rast=$dem
+
+res=`g.region -p | grep res | sed 1d | cut -f2 -d':' | tr -d ' ' `
+
+#g.region vect=$input
+
+#g.region n=n+$res s=s-$res e=e+$res w=w-$res
+
+for j in $thresholds
+do
+
+g.remove vect=ordered_$j
+rm orderedtxt_$j
+
+echo "initializing statistics for threshold value=$j"
+
+r.strahler.sh dem=$dem xcoor=$xcoor ycoor=$ycoor thr=$j output=ordered_$j textoutput=orderedtxt_$j --overwrite
+
+rm one
+rm two
+rm three
+rm maxord
+
+`cat orderedtxt_$j | sed 2d |sed 1d > one
+sort -r -k 4,4 one > two
+head -1 two > three
+awk 'NR==1{print $4}' three > "maxord"`
+
+maxord=`cat maxord`
+
+rm thrVSmaxord_$j
+
+echo "$j $maxord" > thrVSmaxord_$j
+
+ a=1
+ until [ "$a" -gt "$maxord" ]
+ do
+ g.remove vect=ord1,ord1_pl,ord1_pl_3d
+ g.remove vect=ord1_pl_3d_cat
+ rm tmp
+ rm ordcol
+ rm textout
+
+ v.extract input=ordered_$j output=ord1 type=line layer=1 new=-1 list=$a
+
+ v.build.polylines input=ord1 output=ord1_pl cats=no
+
+ v.drape input=ord1_pl type=line rast=$dem method=nearest output=ord1_pl_3d
+
+ v.category input=ord1_pl_3d output=ord1_pl_3d_cat type=line option=add cat=1 layer=1 step=1
+
+ v.db.addtable map=ord1_pl_3d_cat table=ord1_pl_3d_cat layer=1 'columns=cat integer, sx double, sy double, sz double, ex double, ey double, ez double'
+
+ v.to.db map=ord1_pl_3d_cat type=line layer=1 qlayer=1 option=start units=meters 'column=sx, sy, sz'
+
+ v.to.db map=ord1_pl_3d_cat type=line layer=1 qlayer=1 option=end units=meters 'column=ex, ey, ez'
+
+ v.rast.stats vector=ord1_pl_3d_cat raster=$dem colprefix=stats percentile=90 --verbose -c
+
+
+ db_driver=`db.connect -p | grep driver | cut -f2 -d':'`
+ db_database=`db.connect -p | grep database | cut -f2 -d':'`
+ for i in `db.select table=ord1_pl_3d_cat database=$db_database driver=$db_driver 'sql=select cat from ord1_pl_3d_cat' | sed 1d`
+ do
+
+ min=`echo "select stats_min from ord1_pl_3d_cat where cat=$i" | db.select | sed 1d`
+ sz=`echo "select sz from ord1_pl_3d_cat where cat=$i" | db.select | sed 1d`
+ ez=`echo "select ez from ord1_pl_3d_cat where cat=$i" | db.select | sed 1d`
+ if [ `echo "$ez - $min" | bc -l | cut -f1 -d'.'` -lt "$lt" ] || [ `echo "$sz - $min" | bc -l | cut -f1 -d'.'` -lt "$lt" ]
+ then
+ if [ "$ez" -gt "$sz" ]
+ then
+ echo "$ez - $min" | bc -l | cut -f1 -d'.' >> tmp
+ echo $a >> ordcol
+ else
+ echo "$sz - $min" | bc -l | cut -f1 -d'.' >> tmp
+ echo $a >> ordcol
+ fi
+ else
+ echo "$ez - $min" | bc -l | cut -f1 -d'.' >> tmp
+ echo $a >> ordcol
+ echo "$sz - $min" | bc -l | cut -f1 -d'.' >> tmp
+ echo $a >> ordcol
+ fi
+ done
+ paste ordcol tmp > textout\_$a
+ #read
+
+ if [ "$a" -gt 1 ]
+ then
+ b=`echo "$a-1" | bc -l`
+ cat textout\_$b textout\_$a > textout_temp
+ rm textout\_$a textout\_$b
+ mv textout_temp textout\_$a
+ fi
+
+ a=$(($a+1))
+ done
+
+ mv textout\_$maxord textout
+ echo "plotting textout..."
+ cat textout
+
+ echo "
+ imported=read.csv('textout', sep = '\t', dec='.', header = F)
+ dframe=data.frame(imported)
+ str(dframe)
+ nrow_maxord=nrow(dframe[dframe\$V1==max(dframe\$V1),])
+ maxord=max(dframe\$V1)
+ ifelse(nrow_maxord==1, (dframe=na.omit(subset(dframe[dframe\$V1 < maxord,]))), (dframe))
+ cat('plotted dframe \n')
+ nord=c(1:max(dframe\$V1))
+ #search for negative values of drop, replacing them with 0.1 values
+ dframe\$V2[(dframe\$V2 < 0)]<-0.1
+ tmp<-dframe
+ #search for zero values of drop, replacing them with 0.1 values
+ tmp\$V2[dframe\$V2==0]<-0.1
+ dframe<-tmp
+ delta=c()
+ for(i in nord)
+ {
+ offset=(qt(0.025,nrow(dframe[dframe\$V1==i,])-1,lower.tail=FALSE)*(sd(log(dframe\$V2[dframe\$V1==i])))/sqrt(nrow(dframe[dframe\$V1==i,])))/(mean (log(dframe\$V2[dframe\$V1==i])))
+ delta=c(delta,offset)
+
+ }
+ delta=delta[delta<0.5]
+ delta<-data.frame(na.omit(delta))
+ dframe<-data.frame(na.omit(dframe[dframe\$V1 <= (nrow(delta)),]))
+ print(dframe)
+ cat('plotted input dataframe \n')
+ linear_regr=summary(lm(dframe\$V2~dframe\$V1))
+ print(linear_regr)
+ cat('plotted linear regression output \n')
+ t=linear_regr\$coefficients[2,3]
+ Pr=linear_regr\$coefficients[2,4]
+ Radj=linear_regr\$adj.r.squared
+ library(agricolae)
+ wer_t=waerden.test(dframe\$V2,dframe\$V1,group=F)
+ print(wer_t)
+ cat('plotted Van der Waerden test output \n')
+ sink(file=\"tmpfile\")
+ waerden.test(dframe\$V2,dframe\$V1,group=F)
+ sink()
+ Pval=try(system(\"cat tmpfile | grep Pvalue | cut -f2 -d' '\", intern=TRUE))
+ system(\"rm tmpfile\")
+ Pval<-as.numeric(Pval)
+ out_row<-as.numeric(c(\"$j\",t,Pr,Radj,Pval))
+ cat('plotting values: threshold, t, Pr, Radj, Pval \n')
+ print(out_row)
+ write(out_row,file='results\_$j',ncolumns=5,sep='\t')
+ " > R_temp
+
+ echo 'source ("R_temp")' | R --vanilla --slave
+
+ rm R_temp
+
+done
+
+rm $result
+rm thrVSmaxord
+
+echo "threshold t Pr Radj Pvalue">$result
+echo "threshold maxord">thrVSmaxord
+
+for k in $thresholds
+do
+ cat results_$k >> $result
+ rm results_$k
+
+ cat thrVSmaxord_$k >> thrVSmaxord
+ rm thrVSmaxord_$k
+done
+
+echo "plotting results..."
+cat $result
+
+echo "
+thrVSmaxord=read.csv(\"thrVSmaxord\", sep = '\t', dec='.', header = F)
+names(thrVSmaxord)[1]=\"Threshold\"
+names(thrVSmaxord)[2]=\"Maxord\"
+xgraph=read.csv(\"$result\", sep = '\t', dec='.', header = F)
+names(xgraph)[1]=\"Threshold\"
+names(xgraph)[2]=\"t_statistic\"
+names(xgraph)[3]=\"Pr\"
+names(xgraph)[4]=\"R_squared_adj\"
+names(xgraph)[5]=\"Pvalue\"
+
+pdf('waerden_test.pdf')
+print(matplot(xgraph\$Threshold, xgraph[, c(\"Pvalue\")], type=\"l\", col=\"red\", lwd=3, lty=1, ylab=\"Pvalue\", pch=1))
+dev.off()
+
+pdf('linear_regression.pdf')
+print(matplot(xgraph\$Threshold, xgraph[, c(\"Pr\")], type=\"l\", col=\"green\", lwd=3, lty=1, ylab=\"Pr\", pch=1))
+dev.off()
+
+pdf('all_tests.pdf')
+print(matplot(xgraph\$Threshold, xgraph[, c(\"Pr\",\"Pvalue\")], type=\"b\", lwd=3, lty=1, ylab=\"(1)Pr, (2)Pvalue\"))
+dev.off()
+
+pdf('thrVSmaxord.pdf')
+print(matplot(thrVSmaxord\$Threshold, thrVSmaxord[, c(\"Maxord\")], type=\"l\", col=\"blue\", lwd=3, lty=1, ylab=\"Maxord\", pch=1))
+dev.off()
+
+" > R_temp
+
+echo 'source ("R_temp")' | R --vanilla --slave
+
+rm R_temp
+
+
Deleted: grass-addons/grass6/vector/v.strahler/r.broscoe.sh
===================================================================
--- grass-addons/grass6/vector/v.strahler/r.broscoe.sh 2012-06-10 23:21:13 UTC (rev 52029)
+++ grass-addons/grass6/vector/v.strahler/r.broscoe.sh 2012-06-11 08:33:16 UTC (rev 52030)
@@ -1,303 +0,0 @@
-#!/bin/sh
-#
-############################################################################
-#
-# MODULE: r.broscoe
-# AUTHOR(S): Annalisa Minelli, Ivan Marchesini
-#
-# PURPOSE: Calculates waerden test and t test statistics
-# for some values of threshold, on a single basin
-# according to A. J. Broscoe theory (1959)
-#
-# COPYRIGHT: (C) 2008 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.
-#
-# REQUIREMENTS: you need R installed with "agricolae" package
-# You can install the package by starting R as root and typing
-# install.packages("agricolae")
-# This will install the package from a CRAN mirror on the internet.
-# For more information on R see the documentation available at
-# http://www.r-project.org/
-#
-# TODO: solve stability problems for low area threshold (cf. v.strahler)
-#############################################################################
-#%Module
-#% description: Calculates waerden test and t test statistics for some values of threshold, on a single basin according to A. J. Broscoe theory (1959)
-#% keywords: waerden.test,t.test,mean stream drop
-#%End
-#%option
-#% key: dem
-#% type: string
-#% key_desc: dem
-#% gisprompt: old,cell,raster
-#% description: Name of DEM raster map
-#% required : yes
-#%END
-#%option
-#% key: thresholds
-#% type: integer
-#% description: Threshold values to calculate statistics on (separated by <space>)
-#% required : yes
-#%end
-#%option
-#% key: xcoor
-#% type: double
-#% description: x coord of outlet
-#% required : yes
-#%end
-#%option
-#% key: ycoor
-#% type: double
-#% description: y coord of outlet
-#% required : yes
-#%end
-#%option
-#% key: lt
-#% type: integer
-#% description: Lesser than (in meters), the program doesn't consider stream drops lesser than this
-#% required : yes
-#%END
-#%option
-#% key: result
-#% type: string
-#% gisprompt: new_file,file,output
-#% key_desc: name
-#% description: Resultant text file to write statistics into
-#% required: yes
-#%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
-
-dem=$GIS_OPT_DEM
-thresholds=$GIS_OPT_THRESHOLDS
-xcoor=$GIS_OPT_XCOOR
-ycoor=$GIS_OPT_YCOOR
-lt=$GIS_OPT_LT
-result=$GIS_OPT_RESULT
-
-### setup enviro vars ###
-eval `g.gisenv`
-: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
-
-echo $LOCATION_NAME
-
-#g.region rast=$dem
-
-res=`g.region -p | grep res | sed 1d | cut -f2 -d':' | tr -d ' ' `
-
-#g.region vect=$input
-
-#g.region n=n+$res s=s-$res e=e+$res w=w-$res
-
-for j in $thresholds
-do
-
-g.remove vect=ordered_$j
-rm orderedtxt_$j
-
-echo "initializing statistics for threshold value=$j"
-
-r.strahler.sh dem=$dem xcoor=$xcoor ycoor=$ycoor thr=$j output=ordered_$j textoutput=orderedtxt_$j --overwrite
-
-rm one
-rm two
-rm three
-rm maxord
-
-`cat orderedtxt_$j | sed 2d |sed 1d > one
-sort -r -k 4,4 one > two
-head -1 two > three
-awk 'NR==1{print $4}' three > "maxord"`
-
-maxord=`cat maxord`
-
-rm thrVSmaxord_$j
-
-echo "$j $maxord" > thrVSmaxord_$j
-
- a=1
- until [ "$a" -gt "$maxord" ]
- do
- g.remove vect=ord1,ord1_pl,ord1_pl_3d
- g.remove vect=ord1_pl_3d_cat
- rm tmp
- rm ordcol
- rm textout
-
- v.extract input=ordered_$j output=ord1 type=line layer=1 new=-1 list=$a
-
- v.build.polylines input=ord1 output=ord1_pl cats=no
-
- v.drape input=ord1_pl type=line rast=$dem method=nearest output=ord1_pl_3d
-
- v.category input=ord1_pl_3d output=ord1_pl_3d_cat type=line option=add cat=1 layer=1 step=1
-
- v.db.addtable map=ord1_pl_3d_cat table=ord1_pl_3d_cat layer=1 'columns=cat integer, sx double, sy double, sz double, ex double, ey double, ez double'
-
- v.to.db map=ord1_pl_3d_cat type=line layer=1 qlayer=1 option=start units=meters 'column=sx, sy, sz'
-
- v.to.db map=ord1_pl_3d_cat type=line layer=1 qlayer=1 option=end units=meters 'column=ex, ey, ez'
-
- v.rast.stats vector=ord1_pl_3d_cat raster=$dem colprefix=stats percentile=90 --verbose -c
-
-
- db_driver=`db.connect -p | grep driver | cut -f2 -d':'`
- db_database=`db.connect -p | grep database | cut -f2 -d':'`
- for i in `db.select table=ord1_pl_3d_cat database=$db_database driver=$db_driver 'sql=select cat from ord1_pl_3d_cat' | sed 1d`
- do
-
- min=`echo "select stats_min from ord1_pl_3d_cat where cat=$i" | db.select | sed 1d`
- sz=`echo "select sz from ord1_pl_3d_cat where cat=$i" | db.select | sed 1d`
- ez=`echo "select ez from ord1_pl_3d_cat where cat=$i" | db.select | sed 1d`
- if [ `echo "$ez - $min" | bc -l | cut -f1 -d'.'` -lt "$lt" ] || [ `echo "$sz - $min" | bc -l | cut -f1 -d'.'` -lt "$lt" ]
- then
- if [ "$ez" -gt "$sz" ]
- then
- echo "$ez - $min" | bc -l | cut -f1 -d'.' >> tmp
- echo $a >> ordcol
- else
- echo "$sz - $min" | bc -l | cut -f1 -d'.' >> tmp
- echo $a >> ordcol
- fi
- else
- echo "$ez - $min" | bc -l | cut -f1 -d'.' >> tmp
- echo $a >> ordcol
- echo "$sz - $min" | bc -l | cut -f1 -d'.' >> tmp
- echo $a >> ordcol
- fi
- done
- paste ordcol tmp > textout\_$a
- #read
-
- if [ "$a" -gt 1 ]
- then
- b=`echo "$a-1" | bc -l`
- cat textout\_$b textout\_$a > textout_temp
- rm textout\_$a textout\_$b
- mv textout_temp textout\_$a
- fi
-
- a=$(($a+1))
- done
-
- mv textout\_$maxord textout
- echo "plotting textout..."
- cat textout
-
- echo "
- imported=read.csv('textout', sep = '\t', dec='.', header = F)
- dframe=data.frame(imported)
- str(dframe)
- nrow_maxord=nrow(dframe[dframe\$V1==max(dframe\$V1),])
- maxord=max(dframe\$V1)
- ifelse(nrow_maxord==1, (dframe=na.omit(subset(dframe[dframe\$V1 < maxord,]))), (dframe))
- cat('plotted dframe \n')
- nord=c(1:max(dframe\$V1))
- #search for negative values of drop, replacing them with 0.1 values
- dframe\$V2[(dframe\$V2 < 0)]<-0.1
- tmp<-dframe
- #search for zero values of drop, replacing them with 0.1 values
- tmp\$V2[dframe\$V2==0]<-0.1
- dframe<-tmp
- delta=c()
- for(i in nord)
- {
- offset=(qt(0.025,nrow(dframe[dframe\$V1==i,])-1,lower.tail=FALSE)*(sd(log(dframe\$V2[dframe\$V1==i])))/sqrt(nrow(dframe[dframe\$V1==i,])))/(mean (log(dframe\$V2[dframe\$V1==i])))
- delta=c(delta,offset)
-
- }
- delta=delta[delta<0.5]
- delta<-data.frame(na.omit(delta))
- dframe<-data.frame(na.omit(dframe[dframe\$V1 <= (nrow(delta)),]))
- print(dframe)
- cat('plotted input dataframe \n')
- linear_regr=summary(lm(dframe\$V2~dframe\$V1))
- print(linear_regr)
- cat('plotted linear regression output \n')
- t=linear_regr\$coefficients[2,3]
- Pr=linear_regr\$coefficients[2,4]
- Radj=linear_regr\$adj.r.squared
- library(agricolae)
- wer_t=waerden.test(dframe\$V2,dframe\$V1,group=F)
- print(wer_t)
- cat('plotted Van der Waerden test output \n')
- sink(file=\"tmpfile\")
- waerden.test(dframe\$V2,dframe\$V1,group=F)
- sink()
- Pval=try(system(\"cat tmpfile | grep Pvalue | cut -f2 -d' '\", intern=TRUE))
- system(\"rm tmpfile\")
- Pval<-as.numeric(Pval)
- out_row<-as.numeric(c(\"$j\",t,Pr,Radj,Pval))
- cat('plotting values: threshold, t, Pr, Radj, Pval \n')
- print(out_row)
- write(out_row,file='results\_$j',ncolumns=5,sep='\t')
- " > R_temp
-
- echo 'source ("R_temp")' | R --vanilla --slave
-
- rm R_temp
-
-done
-
-rm $result
-rm thrVSmaxord
-
-echo "threshold t Pr Radj Pvalue">$result
-echo "threshold maxord">thrVSmaxord
-
-for k in $thresholds
-do
- cat results_$k >> $result
- rm results_$k
-
- cat thrVSmaxord_$k >> thrVSmaxord
- rm thrVSmaxord_$k
-done
-
-echo "plotting results..."
-cat $result
-
-echo "
-thrVSmaxord=read.csv(\"thrVSmaxord\", sep = '\t', dec='.', header = F)
-names(thrVSmaxord)[1]=\"Threshold\"
-names(thrVSmaxord)[2]=\"Maxord\"
-xgraph=read.csv(\"$result\", sep = '\t', dec='.', header = F)
-names(xgraph)[1]=\"Threshold\"
-names(xgraph)[2]=\"t_statistic\"
-names(xgraph)[3]=\"Pr\"
-names(xgraph)[4]=\"R_squared_adj\"
-names(xgraph)[5]=\"Pvalue\"
-
-pdf('waerden_test.pdf')
-print(matplot(xgraph\$Threshold, xgraph[, c(\"Pvalue\")], type=\"l\", col=\"red\", lwd=3, lty=1, ylab=\"Pvalue\", pch=1))
-dev.off()
-
-pdf('linear_regression.pdf')
-print(matplot(xgraph\$Threshold, xgraph[, c(\"Pr\")], type=\"l\", col=\"green\", lwd=3, lty=1, ylab=\"Pr\", pch=1))
-dev.off()
-
-pdf('all_tests.pdf')
-print(matplot(xgraph\$Threshold, xgraph[, c(\"Pr\",\"Pvalue\")], type=\"b\", lwd=3, lty=1, ylab=\"(1)Pr, (2)Pvalue\"))
-dev.off()
-
-pdf('thrVSmaxord.pdf')
-print(matplot(thrVSmaxord\$Threshold, thrVSmaxord[, c(\"Maxord\")], type=\"l\", col=\"blue\", lwd=3, lty=1, ylab=\"Maxord\", pch=1))
-dev.off()
-
-" > R_temp
-
-echo 'source ("R_temp")' | R --vanilla --slave
-
-rm R_temp
-
-
Deleted: grass-addons/grass6/vector/v.strahler/r.broscoe.sh.html
===================================================================
--- grass-addons/grass6/vector/v.strahler/r.broscoe.sh.html 2012-06-10 23:21:13 UTC (rev 52029)
+++ grass-addons/grass6/vector/v.strahler/r.broscoe.sh.html 2012-06-11 08:33:16 UTC (rev 52030)
@@ -1,108 +0,0 @@
-<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
-<html>
-<head>
-<title>GRASS GIS: r.broscoe.sh</title>
-<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
-<link rel="stylesheet" href="grassdocs.css" type="text/css">
-</head>
-<body bgcolor="white">
-<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
-<h2>NAME</h2>
-<em><b>r.broscoe.sh</b></em>
-
-<h2>DESCRIPTION</h2>
-
-<em>r.broscoe.sh</em> Calculates waerden test and t test statistics for some values of threshold area on a single basin, according to A.J.Broscoe theory (1959).<br />
-The program uses some <em><a href="http://www.r-project.org/">R</a></em> commands for statistical analisys and graphic rapresentation. In particoular the R package <em>"agricolae"</em> is required.<br />
-The A.J.Broscoe theory is well known as the theory of the "Mean Stream Drop" and it says that, for the extraction by DEM of a stream network, exists a threshold value wich makes <em>drop</em> constant, and this is the <em>right</em> one extraction threshold. <br />
-By definig the <em>drop</em> (H) as:<br />
-<br />
-<em>H = S L</em><br />
-<br />
-for streams in each Strahler order. <br />
-The <em>drop</em> is calculated as the quote difference by ending and starting point of a stream; S is the slope and L is the lenght of the same stream.<br />
-Using the Leopold and Miller relation (1964):<br />
-<br />
-<em>S=CA<sup>t</sup></em><br />
-<br />
-where C and t are constat values; the area we are searching (A) is the lowest that gives S to find H costant for each Strahler order (w).<br />
-<br />
-<em>H<sub>w</sub> = H<sub>w+1</sub> = H<sub>w+2</sub> = ...</em><br />
-<br />
-where H<sub>w</sub> is the <em>mean</em> of the drops related to the streams in the same Strahler order (w).<br />
-The area can be found by making some attempts for different area thresholds, doing some statistical tests (Van der Waerden test and linear regression), and choosing the <em>right</em> threshold from the output of the tests.<br />
-<br />
-<em>r.broscoe.sh</em> takes in input the DEM, the threshold values on wich calculate statistics, the outlet coords of the basin you want to study; it returns a table (text file) with the output of the Van der Waerden test and linear regression (t test) for each threshold value.<br />
-For the Van der Waerden test the parameter <em>Pvalue</em> is taken. It has to be greater than the possible, it represents the possibility of success of the test (the <em>Mean Stream Drop</em> is the same for all Strahler orders).<br />
-For the linear regression the parameters <em>t, Pr, R_squared_adj</em> are taken. <em>t</em> is the t statistic value, <em>Pr</em> is the possibility of success of the t test, <em>R_squared_adj</em> measures the dispersion of data around the mean value (for each order) for given degrees of freedom.<br />
-Three graphics called "linear_regression", "waerden_test" and "all_tests" are also generated as PDF in the home folder.<br />
-<br />
-Preferably let's take the threshold value wich gives <em>Pvalue</em> (or <em>Pr</em>) greater than 0.95, but is not granted that you can reach that result because it depends of the well-graduation (by Horton-Strahler) of the basin, its geomorphological maturity, so it is not rare that you have to take threshold where <em>Pvalue</em> is simlpy the greatest.<br />
-At the end of the calculation, at first <em>Pvalue</em> is examinated, then, only if Van der Waerden test gives no good results (low <em>Pvalue</em>), the linear regression output (<em>Pr</em>) is examinated; in fact the Van der Waerden test is preferred to linear regression because it allows you to consider the real dispersion of data around the mean: this makes you able to know the real significance of the probability (e.g. the significance is low for few data in the sample) considering an unique parameter.<br />
-<br />
-<h2>EXAMPLE</h2>
-
-An example on Menotre stream (Umbria, Italy):<br />
-The syntax:
-
-<div class="code"><pre>
- r.broscoe.sh dem=dtm20_regione at AB 'thresholds=400 600 800 1000 1200 1400 1800 2000' xcoor=2291350.34 ycoor=4765192.22 lt=4 result=menotre_txt
-</pre></div>
-
-The results:
-<div class="code"><pre>
-threshold t Pr Radj Pvalue
-400 0.5713518 0.568486 -0.003798402 0.6085511
-600 0.8791352 0.3810997 -0.001896266 0.2798474
-800 1.053110 0.2948033 0.001067895 0.29454
-1000 0.02578308 0.9794938 -0.01233737 0.8535388
-1200 0.3985548 0.69147 -0.01234108 0.6340721
-1400 -1.024254 0.3100425 0.0008457844 0.256408
-1800 -0.6368832 0.5274277 -0.01309044 0.5764749
-2000 -0.4003206 0.6908575 -0.01901582 0.814699
-</pre></div>
-<br />
-<img src="images/wt_rbroscoe.jpg"> <img src="images/lr_rbroscoe.jpg"> <img src="images/at_rbroscoe.jpg"> <br />
-<br />
-By the report and graphics, you can see that the Van der Werden test gives not-so-good results (<em>Pvalue_max</em>=0.85 for threshold=1000 cells) but, if you consider the linear regression output (<em>Pr</em>), you can see that for the same threshold value (1000 cells) <em>Pr</em> is 97%.<br />
-So the threshold=1000 cells is chosen. Moreover the program returns a set of vector map called <em>"orderd_thresholdvalue"</em> from wich you can extract the right one orderd-network (in this case the right one is <em>"orderd_1000"</em>), you can rename and use it as well as you want.<br />
-<br />
-<img src="images/menotre.jpg"><br />
-<br />
-<h2>NOTES</h2>
-The <em>lt</em> value requested in input is a parameter that prevents eventual errors in the DEM; it considers the presence of pits and represents the height difference <em>lesserthan</em> a drop is not considered as a drop but as a pit, and extracted from <em>Mean Stream Drop</em> analysis.<br />
-<br />
-The program uses the module <em><a href="r.strahler.html">r.strahler</a></em>, so it presents the same conditions about the selection of a threshold value range.<br />
-
-<h2>SEE ALSO</h2>
-<em><a href="v.strahler.html">v.strahler</a></em><br>
-<em><a href="r.strahler.sh.html">r.strahler.sh</a></em><br>
-
-<h2>REFERENCES</h2>
-NIST, (2006). <i>Van Der Waerden.</i><br />
-URL: <em><a href="http://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/vanderwa.htm">http://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/vanderwa.htm</a></em><br />
-<p>
-D. G. Tarboton and D. P. Ames, (2001). <i>Advances in the mapping of flow networks from digital elevation data.</i><b> World Water and Environment Resources Congress</b>, presentation (2001).<br />
-<p>
-J. J. Flint, (1974). <i>Stream gradient as a function of order, magnitude, and discharge.</i><b> Water Resources Research</b>, vol.10, n.5, p.969-973. <br />
-<p>
-NIST, (2006). <i>Engineering statistical handbook: confidence limits for the mean.</i><br />
-URL: <em><a href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm">http://www.itl.nist.gov/div898/handbook/eda/section3/eda352.htm</a></em><br />
-<p>
-J. C. Davis, (1990). <i>Statistics and Data Analysis in Geology</i>. John Wiley \& Sons editors (New York, NY, USA).<br />
-<p>
-A. J. Broscoe, (1959). <i>Quantitative analysis of longitudinal stream profiles of small watersheds</i>. Department of Geology, Columbia University, NY.<br />
-<p>
-F. De Mendiburu, (2006). <i>Statistical Procedures for Agricultural Research.</i><br />
-URL: <em><a href="http://rss.acs.unt.edu/Rdoc/library/agricolae/html/agricolae.package.html">http://rss.acs.unt.edu/Rdoc/library/agricolae/html/agricolae.package.html</a></em><br />
-
-<h2>AUTHORS</h2>
-
-Ivan Marchesini and Annalisa Minelli, Univ. Perugia. <br>
-
-<p>
-<i>Last changed: $Date$</i>
-<HR>
-<p><a href="index.html">Main index</a> - <a href="raster.html">raster index</a> - <a href="full_index.html">Full index</a></P>
-</body>
-</html>
Added: grass-addons/grass6/vector/v.strahler/r.strahler/Makefile
===================================================================
--- grass-addons/grass6/vector/v.strahler/r.strahler/Makefile (rev 0)
+++ grass-addons/grass6/vector/v.strahler/r.strahler/Makefile 2012-06-11 08:33:16 UTC (rev 52030)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.strahler
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Copied: grass-addons/grass6/vector/v.strahler/r.strahler/description.html (from rev 51972, grass-addons/grass6/vector/v.strahler/r.strahler.sh.html)
===================================================================
--- grass-addons/grass6/vector/v.strahler/r.strahler/description.html (rev 0)
+++ grass-addons/grass6/vector/v.strahler/r.strahler/description.html 2012-06-11 08:33:16 UTC (rev 52030)
@@ -0,0 +1,39 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.strahler.sh</em> Creates a vector map of Strahler ordered streams of a single basin starting from a DEM.
+Once selected the area where you want to work by visualizing it in the GRASS monitor, the programs extracts lines form the DEM by using <em><a href="r.watershed.html">r.watershed</a></em> , makes you able to select an outlet for the basin (or takes the outlet coords manually inserted at the beginning), individuates the basin by using <em><a href="r.water.outlet.html">r.water.outlet</a></em> , cleans the topology, executes <em><a href="v.strahler.html">v.strahler</a></em> for the extracted network. <br />
+The program returns an ouput text file (see <em><a href="v.strahler.html">v.strahler</a></em> ) and the basin network ordered by Horton-Strahler.
+The program is also able to plot a background map during the elaboration, so, if you have some aerial photos or the same DEM, it's easier to select the right outlet cross-section in the map.<br />
+If the <b>-i</b> interactive flag is given then you are working in interactive mode and the program makes you able to choose the outlet cross section directly from the monitor, instead the program works in un-interactive mode but it needs the outlet coords manually inserted.
+
+<h2>EXAMPLE</h2>
+
+The interactive-mode syntax:
+
+<div class="code"><pre>
+ r.strahler.sh -i dem=dem_20 thr=1000
+ output=ordered_net textoutput=text_net bkgrmap=dem_20
+</pre></div>
+
+The un-interactive-mode syntax:
+
+<div class="code"><pre>
+ r.strahler.sh dem=dem_20 xcoor=2291350.34 ycoor=4765192.22 thr=1000
+ output=ordered_net textoutput=text_net bkgrmap=dem_20
+</pre></div>
+
+<h2>BUGS</h2>
+The program has some stability problems for low area threshold; from sperimental analysis the program works well for thresholds greater than the ones in graphic (below) for given cells of the studied-area in the DEM.<br />
+<img src="images/thrVScell.jpg"> <br />
+
+<h2>SEE ALSO</h2>
+<em><a href="v.strahler.html">v.strahler</a></em><br>
+<em><a href="r.broscoe.sh.html">r.broscoe.sh</a></em><br>
+
+<h2>AUTHORS</h2>
+
+Ivan Marchesini and Annalisa Minelli, Univ. Perugia. <br>
+
+
+<p>
+<i>Last changed: $Date$</i>
Copied: grass-addons/grass6/vector/v.strahler/r.strahler/r.strahler (from rev 51972, grass-addons/grass6/vector/v.strahler/r.strahler.sh)
===================================================================
--- grass-addons/grass6/vector/v.strahler/r.strahler/r.strahler (rev 0)
+++ grass-addons/grass6/vector/v.strahler/r.strahler/r.strahler 2012-06-11 08:33:16 UTC (rev 52030)
@@ -0,0 +1,222 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: r.strahler.sh
+# AUTHOR(S): Annalisa Minelli, Ivan Marchesini
+# PURPOSE: Create a vector map of strahler ordered streems of a single
+# basin, starting from a DEM.
+# The script extracts from a DEM the network (by using
+# r.watershed), converts the network to vector lines,
+# cleans the topology, lets you decide an outlet for your basin,
+# and executes v.strahler for the upstream basin.
+#
+# COPYRIGHT: (C) 2008 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.
+#
+# TODO: solve stability problems for low area threshold
+#############################################################################
+#%Module
+#% description: Create a vector map of strahler ordered streams of a single basin starting from a DEM
+#% keywords: strahler, streams, vector
+#%End
+#%flag
+#% key: i
+#% description: find interactively the outlet coords
+#%END
+#%option
+#% key: dem
+#% type: string
+#% key_desc: dem
+#% gisprompt: old,cell,raster
+#% description: Name of DEM raster map
+#% required : yes
+#%END
+#%option
+#% key: xcoor
+#% type: double
+#% description: x coord of outlet
+#% required : no
+#%end
+#%option
+#% key: ycoor
+#% type: double
+#% description: y coord of outlet
+#% required : no
+#%end
+#%option
+#% key: thr
+#% type: integer
+#% description: minimum size of exterior watershed basin
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Name of streams ordered map created
+#% required : yes
+#%end
+#%option
+#% key: textoutput
+#% type: string
+#% gisprompt: new_file,file,output
+#% key_desc: name
+#% description: Name for output text file
+#% required: yes
+#%end
+#%option
+#% key: bkgrmap
+#% type: string
+#% key_desc: bkgrmap
+#% gisprompt: old,cell,raster
+#% description: Name of background map to plot
+#% 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
+
+dem=$GIS_OPT_DEM
+xcoor=$GIS_OPT_XCOOR
+ycoor=$GIS_OPT_YCOOR
+thr=$GIS_OPT_THR
+output=$GIS_OPT_OUTPUT
+textoutput=$GIS_OPT_TEXTOUTPUT
+bkgrmap=$GIS_OPT_BKGRMAP
+
+#check presence of x and y coordinates if you want to use the un-interactive mode
+
+if [ $GIS_FLAG_I -eq 0 ] && [ -z $xcoor ] && [ -z $ycoor ]; then
+echo "
+
+
+If you don't want to use the interactive-mode you must select the outlets' coodinates
+
+
+"
+exit 0
+
+fi
+
+#check presence of raster MASK, put it aside
+MASKFOUND=0
+eval `g.findfile element=cell file=MASK`
+if [ "$file" ] ; then
+ g.message "Raster MASK found, temporarily disabled"
+ g.rename MASK,"${TMPNAME}_origmask" --quiet
+ MASKFOUND=1
+fi
+
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+LOCATION=$GISDBASE/$LOCATION_NAME/$MAPSET
+
+#uncomment this line if you need to remove some maps remaining from previous runs
+#g.remove rast=rnetwork,rnetwork_thin,rnetwork_b,rnetwork_l,rnetwork_patch,rnetwork_thin2,drainage,MASK,basin vect=vnetwork,vnetwork_b,vnetwork_b_add,vnetwork2,vnetwork_dangle
+
+#get resolution of DEM
+res=`g.region -p | grep nsres | cut -f2 -d':' | tr -d ' '`
+
+#find raster streams and plot them
+r.watershed elevation=$dem threshold=$thr stream=rnetwork drainage=drainage --overwrite
+r.null map=rnetwork setnull=0
+r.thin input=rnetwork output=rnetwork_th iterations=200 --overwrite
+r.mapcalc "rnetwork_thin=rnetwork_th/rnetwork_th"
+d.mon stop=x0
+d.mon start=x0
+
+if [ "$bkgrmap" ]; then
+d.rast map=$bkgrmap
+fi
+
+
+d.rast -o map=rnetwork_thin
+
+#use the interactive mode
+if [ $GIS_FLAG_I -eq 1 ] ; then
+
+ #ask the user for zooming and choosing the outlet cell
+ #-------
+ echo "
+
+
+ Now, please zoom to the area where you want to put the outlet cross section
+
+
+ "
+
+ d.zoom
+
+ echo "
+
+
+ Now, please click on the cell representing the cross section
+
+
+ "
+
+ #prevent from choosing a wrong cell outside from the streams
+ cat=2
+ while [ "$cat" != "1" ]
+ do
+ result=`d.what.rast -t -1 rnetwork_thin`
+ coor=`echo $result | cut -f1 -d ' '`
+ x=`echo $coor | cut -f1 -d':'`
+ y=`echo $coor | cut -f2 -d':'`
+ cat=`echo $result | cut -f3 -d ' ' | tr -d :`
+ done
+
+ #come back to the previous zoom
+ d.zoom -r
+
+#do not use the interactive mode
+else
+ x=$xcoor
+ y=$ycoor
+fi
+
+
+#find the basin and set the mask
+r.water.outlet drainage=drainage basin=basin easting=$x northing=$y
+d.rast basin
+r.null map=basin setnull=0
+r.mapcalc "MASK=basin"
+d.rast MASK
+
+#clean raster map and convert it to vector
+g.region rast=MASK
+r.to.vect input=rnetwork_thin output=vnetwork feature=line --overwrite
+d.vect map=vnetwork
+v.type input=vnetwork output=vnetwork_b type=line,boundary --overwrite
+v.centroids input=vnetwork_b output=vnetwork_b_add option=add layer=1 cat=1 step=1 --overwrite
+newres=`echo "$res/4" | bc -l`
+g.region res=$newres
+v.to.rast input=vnetwork_b_add output=rnetwork_b use=val type=area layer=1 value=1 rows=4096 --overwrite
+v.to.rast input=vnetwork output=rnetwork_l use=val type=line layer=1 value=1 rows=4096 --overwrite
+r.patch input=rnetwork_b,rnetwork_l output=rnetwork_patch --overwrite
+r.thin input=rnetwork_patch output=rnetwork_thin2 iterations=200 --overwrite
+r.to.vect input=rnetwork_thin2 output=vnetwork2 feature=line --overwrite
+
+#remove dangles and create polylines
+soglia=`echo "sqrt(2*($res^2))+1" | bc -l`
+v.clean input=vnetwork2 output=vnetwork_dangle type=line tool=rmdangle thresh=$soglia --overwrite
+v.build.polylines input=vnetwork_dangle output=vnetwork_poly cats=first --overwrite
+d.vect vnetwork_poly
+
+#enlarge region to ensure v.strahler get raster values at streams ends and use v.strahler
+g.region -a res=$res n=n+$soglia s=s-$soglia e=e+$soglia w=w-$soglia
+v.strahler input=vnetwork_poly output=$output dem=$dem sloppy=0 layer=1 txout=$textoutput
+
+#remove temporary files
+g.remove rast=rnetwork,rnetwork_thin,rnetwork_b,rnetwork_l,rnetwork_patch,rnetwork_thin2,drainage,MASK,basin vect=vnetwork,vnetwork_b,vnetwork_b_add,vnetwork2,vnetwork_dangle
+
Deleted: grass-addons/grass6/vector/v.strahler/r.strahler.sh
===================================================================
--- grass-addons/grass6/vector/v.strahler/r.strahler.sh 2012-06-10 23:21:13 UTC (rev 52029)
+++ grass-addons/grass6/vector/v.strahler/r.strahler.sh 2012-06-11 08:33:16 UTC (rev 52030)
@@ -1,222 +0,0 @@
-#!/bin/sh
-#
-############################################################################
-#
-# MODULE: r.strahler.sh
-# AUTHOR(S): Annalisa Minelli, Ivan Marchesini
-# PURPOSE: Create a vector map of strahler ordered streems of a single
-# basin, starting from a DEM.
-# The script extracts from a DEM the network (by using
-# r.watershed), converts the network to vector lines,
-# cleans the topology, lets you decide an outlet for your basin,
-# and executes v.strahler for the upstream basin.
-#
-# COPYRIGHT: (C) 2008 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.
-#
-# TODO: solve stability problems for low area threshold
-#############################################################################
-#%Module
-#% description: Create a vector map of strahler ordered streams of a single basin starting from a DEM
-#% keywords: strahler, streams, vector
-#%End
-#%flag
-#% key: i
-#% description: find interactively the outlet coords
-#%END
-#%option
-#% key: dem
-#% type: string
-#% key_desc: dem
-#% gisprompt: old,cell,raster
-#% description: Name of DEM raster map
-#% required : yes
-#%END
-#%option
-#% key: xcoor
-#% type: double
-#% description: x coord of outlet
-#% required : no
-#%end
-#%option
-#% key: ycoor
-#% type: double
-#% description: y coord of outlet
-#% required : no
-#%end
-#%option
-#% key: thr
-#% type: integer
-#% description: minimum size of exterior watershed basin
-#% required : yes
-#%end
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,vector,vector
-#% description: Name of streams ordered map created
-#% required : yes
-#%end
-#%option
-#% key: textoutput
-#% type: string
-#% gisprompt: new_file,file,output
-#% key_desc: name
-#% description: Name for output text file
-#% required: yes
-#%end
-#%option
-#% key: bkgrmap
-#% type: string
-#% key_desc: bkgrmap
-#% gisprompt: old,cell,raster
-#% description: Name of background map to plot
-#% 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
-
-dem=$GIS_OPT_DEM
-xcoor=$GIS_OPT_XCOOR
-ycoor=$GIS_OPT_YCOOR
-thr=$GIS_OPT_THR
-output=$GIS_OPT_OUTPUT
-textoutput=$GIS_OPT_TEXTOUTPUT
-bkgrmap=$GIS_OPT_BKGRMAP
-
-#check presence of x and y coordinates if you want to use the un-interactive mode
-
-if [ $GIS_FLAG_I -eq 0 ] && [ -z $xcoor ] && [ -z $ycoor ]; then
-echo "
-
-
-If you don't want to use the interactive-mode you must select the outlets' coodinates
-
-
-"
-exit 0
-
-fi
-
-#check presence of raster MASK, put it aside
-MASKFOUND=0
-eval `g.findfile element=cell file=MASK`
-if [ "$file" ] ; then
- g.message "Raster MASK found, temporarily disabled"
- g.rename MASK,"${TMPNAME}_origmask" --quiet
- MASKFOUND=1
-fi
-
-eval `g.gisenv`
-: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
-LOCATION=$GISDBASE/$LOCATION_NAME/$MAPSET
-
-#uncomment this line if you need to remove some maps remaining from previous runs
-#g.remove rast=rnetwork,rnetwork_thin,rnetwork_b,rnetwork_l,rnetwork_patch,rnetwork_thin2,drainage,MASK,basin vect=vnetwork,vnetwork_b,vnetwork_b_add,vnetwork2,vnetwork_dangle
-
-#get resolution of DEM
-res=`g.region -p | grep nsres | cut -f2 -d':' | tr -d ' '`
-
-#find raster streams and plot them
-r.watershed elevation=$dem threshold=$thr stream=rnetwork drainage=drainage --overwrite
-r.null map=rnetwork setnull=0
-r.thin input=rnetwork output=rnetwork_th iterations=200 --overwrite
-r.mapcalc "rnetwork_thin=rnetwork_th/rnetwork_th"
-d.mon stop=x0
-d.mon start=x0
-
-if [ "$bkgrmap" ]; then
-d.rast map=$bkgrmap
-fi
-
-
-d.rast -o map=rnetwork_thin
-
-#use the interactive mode
-if [ $GIS_FLAG_I -eq 1 ] ; then
-
- #ask the user for zooming and choosing the outlet cell
- #-------
- echo "
-
-
- Now, please zoom to the area where you want to put the outlet cross section
-
-
- "
-
- d.zoom
-
- echo "
-
-
- Now, please click on the cell representing the cross section
-
-
- "
-
- #prevent from choosing a wrong cell outside from the streams
- cat=2
- while [ "$cat" != "1" ]
- do
- result=`d.what.rast -t -1 rnetwork_thin`
- coor=`echo $result | cut -f1 -d ' '`
- x=`echo $coor | cut -f1 -d':'`
- y=`echo $coor | cut -f2 -d':'`
- cat=`echo $result | cut -f3 -d ' ' | tr -d :`
- done
-
- #come back to the previous zoom
- d.zoom -r
-
-#do not use the interactive mode
-else
- x=$xcoor
- y=$ycoor
-fi
-
-
-#find the basin and set the mask
-r.water.outlet drainage=drainage basin=basin easting=$x northing=$y
-d.rast basin
-r.null map=basin setnull=0
-r.mapcalc "MASK=basin"
-d.rast MASK
-
-#clean raster map and convert it to vector
-g.region rast=MASK
-r.to.vect input=rnetwork_thin output=vnetwork feature=line --overwrite
-d.vect map=vnetwork
-v.type input=vnetwork output=vnetwork_b type=line,boundary --overwrite
-v.centroids input=vnetwork_b output=vnetwork_b_add option=add layer=1 cat=1 step=1 --overwrite
-newres=`echo "$res/4" | bc -l`
-g.region res=$newres
-v.to.rast input=vnetwork_b_add output=rnetwork_b use=val type=area layer=1 value=1 rows=4096 --overwrite
-v.to.rast input=vnetwork output=rnetwork_l use=val type=line layer=1 value=1 rows=4096 --overwrite
-r.patch input=rnetwork_b,rnetwork_l output=rnetwork_patch --overwrite
-r.thin input=rnetwork_patch output=rnetwork_thin2 iterations=200 --overwrite
-r.to.vect input=rnetwork_thin2 output=vnetwork2 feature=line --overwrite
-
-#remove dangles and create polylines
-soglia=`echo "sqrt(2*($res^2))+1" | bc -l`
-v.clean input=vnetwork2 output=vnetwork_dangle type=line tool=rmdangle thresh=$soglia --overwrite
-v.build.polylines input=vnetwork_dangle output=vnetwork_poly cats=first --overwrite
-d.vect vnetwork_poly
-
-#enlarge region to ensure v.strahler get raster values at streams ends and use v.strahler
-g.region -a res=$res n=n+$soglia s=s-$soglia e=e+$soglia w=w-$soglia
-v.strahler input=vnetwork_poly output=$output dem=$dem sloppy=0 layer=1 txout=$textoutput
-
-#remove temporary files
-g.remove rast=rnetwork,rnetwork_thin,rnetwork_b,rnetwork_l,rnetwork_patch,rnetwork_thin2,drainage,MASK,basin vect=vnetwork,vnetwork_b,vnetwork_b_add,vnetwork2,vnetwork_dangle
-
Deleted: grass-addons/grass6/vector/v.strahler/r.strahler.sh.html
===================================================================
--- grass-addons/grass6/vector/v.strahler/r.strahler.sh.html 2012-06-10 23:21:13 UTC (rev 52029)
+++ grass-addons/grass6/vector/v.strahler/r.strahler.sh.html 2012-06-11 08:33:16 UTC (rev 52030)
@@ -1,39 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em>r.strahler.sh</em> Creates a vector map of Strahler ordered streams of a single basin starting from a DEM.
-Once selected the area where you want to work by visualizing it in the GRASS monitor, the programs extracts lines form the DEM by using <em><a href="r.watershed.html">r.watershed</a></em> , makes you able to select an outlet for the basin (or takes the outlet coords manually inserted at the beginning), individuates the basin by using <em><a href="r.water.outlet.html">r.water.outlet</a></em> , cleans the topology, executes <em><a href="v.strahler.html">v.strahler</a></em> for the extracted network. <br />
-The program returns an ouput text file (see <em><a href="v.strahler.html">v.strahler</a></em> ) and the basin network ordered by Horton-Strahler.
-The program is also able to plot a background map during the elaboration, so, if you have some aerial photos or the same DEM, it's easier to select the right outlet cross-section in the map.<br />
-If the <b>-i</b> interactive flag is given then you are working in interactive mode and the program makes you able to choose the outlet cross section directly from the monitor, instead the program works in un-interactive mode but it needs the outlet coords manually inserted.
-
-<h2>EXAMPLE</h2>
-
-The interactive-mode syntax:
-
-<div class="code"><pre>
- r.strahler.sh -i dem=dem_20 thr=1000
- output=ordered_net textoutput=text_net bkgrmap=dem_20
-</pre></div>
-
-The un-interactive-mode syntax:
-
-<div class="code"><pre>
- r.strahler.sh dem=dem_20 xcoor=2291350.34 ycoor=4765192.22 thr=1000
- output=ordered_net textoutput=text_net bkgrmap=dem_20
-</pre></div>
-
-<h2>BUGS</h2>
-The program has some stability problems for low area threshold; from sperimental analysis the program works well for thresholds greater than the ones in graphic (below) for given cells of the studied-area in the DEM.<br />
-<img src="images/thrVScell.jpg"> <br />
-
-<h2>SEE ALSO</h2>
-<em><a href="v.strahler.html">v.strahler</a></em><br>
-<em><a href="r.broscoe.sh.html">r.broscoe.sh</a></em><br>
-
-<h2>AUTHORS</h2>
-
-Ivan Marchesini and Annalisa Minelli, Univ. Perugia. <br>
-
-
-<p>
-<i>Last changed: $Date$</i>
Deleted: grass-addons/grass6/vector/v.strahler/v.mainchannel
===================================================================
--- grass-addons/grass6/vector/v.strahler/v.mainchannel 2012-06-10 23:21:13 UTC (rev 52029)
+++ grass-addons/grass6/vector/v.strahler/v.mainchannel 2012-06-11 08:33:16 UTC (rev 52030)
@@ -1,154 +0,0 @@
-#!/bin/sh
-#
-############################################################################
-#
-# MODULE: v.mainchannel
-# AUTHOR(S): Annalisa Minelli, Ivan Marchesini
-# PURPOSE: Find the main channel of a vector stream network (useful to compare with a network ordered by Horton-Strahler)
-# COPYRIGHT: (C) 2005 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: Find the main channel of a vector stream network (useful to compare with a network ordered by Horton-Strahler)
-#% keywords: strahler, main_channel, vector
-#%End
-#%option
-#% key: dem
-#% type: string
-#% key_desc: dem
-#% gisprompt: old,cell,raster
-#% description: Name of DEM raster map
-#% required : yes
-#%END
-#%option
-#% key: xcoor
-#% type: double
-#% description: x coord of outlet (read by dem in the pixel where is the vector file cross section)
-#% required : yes
-#%end
-#%option
-#% key: ycoor
-#% type: double
-#% description: y coord of outlet (read by dem in the pixel where is the vector file cross section)
-#% required : yes
-#%end
-#%option
-#% key: cost
-#% type: integer
-#% description: Distance between isolines (separated by <comma>)
-#% required : yes
-#%end
-#%option
-#% key: input
-#% type: string
-#% gisprompt: old,vector,vector
-#% description: Vector map of Strahler-ordered streams
-#% required : yes
-#%end
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,vector,vector
-#% description: Name of the main_channel vector map
-#% required : yes
-#%END
-#%flag
-#% key: d
-#% label: Display results
-#% description: Plot the result of elaboration in relation to the strahler ordered stream network
-#%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
-
-dem="$GIS_OPT_DEM"
-xcoor="$GIS_OPT_XCOOR"
-ycoor="$GIS_OPT_YCOOR"
-cost="$GIS_OPT_COST"
-input="$GIS_OPT_INPUT"
-output="$GIS_OPT_OUTPUT"
-
-### setup enviro vars ###
-eval `g.gisenv`
-: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
-
-g.remove vect=end_basin,out_1del,out_1add,out_basin,out_isodistance,out_2,out_2del,out_2points,out_2addp,out_3D,out_path
-
-echo "$xcoor|$ycoor|1" | v.in.ascii out=end_basin
-
-thr=`g.region -p | grep nsres | cut -f2 -d':' | tr -d ' '`
-thresh=`echo "$thr*1.5" | bc -l | cut -f1 -d'.'`
-
-
-#### setup temporary files
-TMP="`g.tempfile pid=$$`"
-if [ $? -ne 0 ] || [ -z "$TMP" ] ; then
- g.message -e "Unable to create temporary files"
- exit 1
-fi
-
-v.category input=$input output=out_1del type=line option=del cat=1 layer=1 step=1
-v.category input=out_1del output=out_1add type=line option=add cat=1 layer=1 step=1
-v.net input=out_1add points=end_basin output=out_basin operation=connect \
- alayer=1 nlayer=2 thresh="$thresh"
-v.net.iso input=out_basin output=out_isodistance type=line \
- alayer=1 nlayer=2 ccats=1-99999999 costs=$cost
-
-echo "v.net.iso executed"
-
-v.category input=out_isodistance opt=report -g > "$TMP.1"
-grep 'line' "$TMP.1" | cut -f5 -d ' ' > "$TMP.2"
-
-v.extract input=out_isodistance output=out_2 type=line layer=1 new=-1 file="$TMP.2"
-v.category input=out_2 output=out_2del option=del \
- type=point,line,boundary,centroid,area cat=1 layer=1 step=1
-v.net input=out_2del output=out_2points operation=nodes alayer=1 nlayer=2
-v.category input=out_2points output=out_2addp type=point option=add \
- cat=1 layer=2 step=1
-
-v.drape input=out_2addp type=point,line rast=$dem method=nearest output=out_3D
-v.db.addtable map=out_3D table=tab_coord layer=2 \
- 'columns=cat integer, x double precision, y double precision, z integer'
-v.to.db map=out_3D type=point layer=2 qlayer=1 option=coor units=meters 'column=x, y, z'
-v.univar -g map=out_3D type=point column=z layer=2 percentile=90 > "$TMP.3"
-
-max=`cat "$TMP.3" | grep max | cut -f2 -d'='`
-
-db_driver=`db.connect -p | grep driver | cut -f2 -d':'`
-db_database=`db.connect -p | grep database | cut -f2 -d':'`
-x_max=`db.select -c table=tab_coord database=$db_database driver=$db_driver 'sql=select x from tab_coord where z='$max''`
-y_max=`db.select -c table=tab_coord database=$db_databese driver=$db_driver 'sql=select y from tab_coord where z='$max''`
-
-echo "x_max e y_max valgono: $x_max e $y_max"
-
-echo "1 $x_max $y_max $xcoor $ycoor" | v.net.path input=out_isodistance \
- output=$output type=line alayer=1 nlayer=2
-
-
-if [ "$GIS_FLAG_D" -eq 1 ] ; then
- #plot the result of elaboration in relation to the strahler ordered stream network
- #d.mon stop=x0
- d.mon start=x0
- d.vect "$input" display=shape,cat
- d.vect map="$output" color=red lcolor=black fcolor=170:170:170 \
- display=shape type=line icon=basic/x size=5 width=2 layer=1 lsize=8
- d.barscale tcolor=black bcolor=white at=2,2
-fi
-
-
-#g.remove vect=end_basin,out_1del,out_1add,out_basin,out_2,out_2del,out_2points,out_2addp,out_3D
-
-\rm "$TMP" "$TMP".[1-3]
-
-#you can modify the "cost" distances in order to have main channel including the max_ordered channel
Copied: grass-addons/grass6/vector/v.strahler/v.mainchannel.sh (from rev 51972, grass-addons/grass6/vector/v.strahler/v.mainchannel)
===================================================================
--- grass-addons/grass6/vector/v.strahler/v.mainchannel.sh (rev 0)
+++ grass-addons/grass6/vector/v.strahler/v.mainchannel.sh 2012-06-11 08:33:16 UTC (rev 52030)
@@ -0,0 +1,154 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: v.mainchannel
+# AUTHOR(S): Annalisa Minelli, Ivan Marchesini
+# PURPOSE: Find the main channel of a vector stream network (useful to compare with a network ordered by Horton-Strahler)
+# COPYRIGHT: (C) 2005 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: Find the main channel of a vector stream network (useful to compare with a network ordered by Horton-Strahler)
+#% keywords: strahler, main_channel, vector
+#%End
+#%option
+#% key: dem
+#% type: string
+#% key_desc: dem
+#% gisprompt: old,cell,raster
+#% description: Name of DEM raster map
+#% required : yes
+#%END
+#%option
+#% key: xcoor
+#% type: double
+#% description: x coord of outlet (read by dem in the pixel where is the vector file cross section)
+#% required : yes
+#%end
+#%option
+#% key: ycoor
+#% type: double
+#% description: y coord of outlet (read by dem in the pixel where is the vector file cross section)
+#% required : yes
+#%end
+#%option
+#% key: cost
+#% type: integer
+#% description: Distance between isolines (separated by <comma>)
+#% required : yes
+#%end
+#%option
+#% key: input
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Vector map of Strahler-ordered streams
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Name of the main_channel vector map
+#% required : yes
+#%END
+#%flag
+#% key: d
+#% label: Display results
+#% description: Plot the result of elaboration in relation to the strahler ordered stream network
+#%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
+
+dem="$GIS_OPT_DEM"
+xcoor="$GIS_OPT_XCOOR"
+ycoor="$GIS_OPT_YCOOR"
+cost="$GIS_OPT_COST"
+input="$GIS_OPT_INPUT"
+output="$GIS_OPT_OUTPUT"
+
+### setup enviro vars ###
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+
+g.remove vect=end_basin,out_1del,out_1add,out_basin,out_isodistance,out_2,out_2del,out_2points,out_2addp,out_3D,out_path
+
+echo "$xcoor|$ycoor|1" | v.in.ascii out=end_basin
+
+thr=`g.region -p | grep nsres | cut -f2 -d':' | tr -d ' '`
+thresh=`echo "$thr*1.5" | bc -l | cut -f1 -d'.'`
+
+
+#### setup temporary files
+TMP="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMP" ] ; then
+ g.message -e "Unable to create temporary files"
+ exit 1
+fi
+
+v.category input=$input output=out_1del type=line option=del cat=1 layer=1 step=1
+v.category input=out_1del output=out_1add type=line option=add cat=1 layer=1 step=1
+v.net input=out_1add points=end_basin output=out_basin operation=connect \
+ alayer=1 nlayer=2 thresh="$thresh"
+v.net.iso input=out_basin output=out_isodistance type=line \
+ alayer=1 nlayer=2 ccats=1-99999999 costs=$cost
+
+echo "v.net.iso executed"
+
+v.category input=out_isodistance opt=report -g > "$TMP.1"
+grep 'line' "$TMP.1" | cut -f5 -d ' ' > "$TMP.2"
+
+v.extract input=out_isodistance output=out_2 type=line layer=1 new=-1 file="$TMP.2"
+v.category input=out_2 output=out_2del option=del \
+ type=point,line,boundary,centroid,area cat=1 layer=1 step=1
+v.net input=out_2del output=out_2points operation=nodes alayer=1 nlayer=2
+v.category input=out_2points output=out_2addp type=point option=add \
+ cat=1 layer=2 step=1
+
+v.drape input=out_2addp type=point,line rast=$dem method=nearest output=out_3D
+v.db.addtable map=out_3D table=tab_coord layer=2 \
+ 'columns=cat integer, x double precision, y double precision, z integer'
+v.to.db map=out_3D type=point layer=2 qlayer=1 option=coor units=meters 'column=x, y, z'
+v.univar -g map=out_3D type=point column=z layer=2 percentile=90 > "$TMP.3"
+
+max=`cat "$TMP.3" | grep max | cut -f2 -d'='`
+
+db_driver=`db.connect -p | grep driver | cut -f2 -d':'`
+db_database=`db.connect -p | grep database | cut -f2 -d':'`
+x_max=`db.select -c table=tab_coord database=$db_database driver=$db_driver 'sql=select x from tab_coord where z='$max''`
+y_max=`db.select -c table=tab_coord database=$db_databese driver=$db_driver 'sql=select y from tab_coord where z='$max''`
+
+echo "x_max e y_max valgono: $x_max e $y_max"
+
+echo "1 $x_max $y_max $xcoor $ycoor" | v.net.path input=out_isodistance \
+ output=$output type=line alayer=1 nlayer=2
+
+
+if [ "$GIS_FLAG_D" -eq 1 ] ; then
+ #plot the result of elaboration in relation to the strahler ordered stream network
+ #d.mon stop=x0
+ d.mon start=x0
+ d.vect "$input" display=shape,cat
+ d.vect map="$output" color=red lcolor=black fcolor=170:170:170 \
+ display=shape type=line icon=basic/x size=5 width=2 layer=1 lsize=8
+ d.barscale tcolor=black bcolor=white at=2,2
+fi
+
+
+#g.remove vect=end_basin,out_1del,out_1add,out_basin,out_2,out_2del,out_2points,out_2addp,out_3D
+
+\rm "$TMP" "$TMP".[1-3]
+
+#you can modify the "cost" distances in order to have main channel including the max_ordered channel
More information about the grass-commit
mailing list