[GRASS-SVN] r54721 - in grass-addons/grass6/vector: . v.breach
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Jan 20 10:17:18 PST 2013
Author: msieczka
Date: 2013-01-20 10:17:17 -0800 (Sun, 20 Jan 2013)
New Revision: 54721
Added:
grass-addons/grass6/vector/v.breach/
grass-addons/grass6/vector/v.breach/Makefile
grass-addons/grass6/vector/v.breach/description.html
grass-addons/grass6/vector/v.breach/v.breach
Log:
Submitting v.breach addon shell script.
Added: grass-addons/grass6/vector/v.breach/Makefile
===================================================================
--- grass-addons/grass6/vector/v.breach/Makefile (rev 0)
+++ grass-addons/grass6/vector/v.breach/Makefile 2013-01-20 18:17:17 UTC (rev 54721)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.breach
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass6/vector/v.breach/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass6/vector/v.breach/description.html
===================================================================
--- grass-addons/grass6/vector/v.breach/description.html (rev 0)
+++ grass-addons/grass6/vector/v.breach/description.html 2013-01-20 18:17:17 UTC (rev 54721)
@@ -0,0 +1,65 @@
+<h2>DESCRIPTION</h2>
+
+The <em>out_ln</em> vector map contains all input <em>vect</em> lines segmented so that each line goes per one input <em>rast</em> cell per one input <em>vect</em> line. Line directions are preserved. There are 2 layers in <em>out_ln</em>.
+<p>
+Categories in layer 1 of the <em>out_ln</em> lines are unchanged compared to <em>vect</em>. In layer 2, each line has a unique category. These categories increase at a step of 1 along each input line, down the input line direction. Following attributes are stored for each line in layer 2:
+<p>
+<b>lcat</b> - category of the input <em>vect</em> line <br>
+<b>z</b> - original elevation of the <em>rast</em> cell under the given <em>out_ln</em> line <br>
+<b>z_breach</b> - elevation of the <em>rast</em> cell under the given <em>out_ln</em> line <b>breached</b>, so that the elevation gradient down the watercourse is assured; <b>minus</b> additional <em>depth</em> specified <br>
+<p>
+The <b>z_breach</b> is calculated in a following way: One line from the <em>vect</em> input is processed at a time. It is divided into segments, one per each input <em>rast</em> cell it flows through. Segments are being read down the line direction. If the elevation of the cell under the current segment >= elevation of the previous cell, the elevation of the line segment = the elevation of the previous cell - 0.000001 m. Otherwise, the elevation of the current segment = the elevation of the current cell. 0.000001 is the lowest value possible to use, because GRASS database modules round down everything beyond 6 decimal places (at least <em>v.db.select</em> and <em>db.execute</em> used in the script do).
+<p>
+There is one point in the <em>out_pt</em> per one <em>out_ln</em> line. Categories in layer 1 and 2 of <em>out_pt</em> are identical to the cats of the closest <em>out_ln</em> line. Points are located in the middles of <em>out_ln</em> lines, except the first and last point of the given cat in layer 1, which are located, respectively, at the beginning and at the end of the given input <em>vect</em> line. For this reason, where 2 or more input <em>vect</em> lines are connected, the same number of points lie on each other. Attributes stored are the same as in the layer 2 of <em>out_ln</em>, plus:
+<p>
+<b>along</b> - the distance from the starting node of the given input <em>vect</em> line <br>
+<b>x</b> - X coordinate <br>
+<b>y</b> - Y coordinate
+<p>
+The <em>out_pt</em>'s <b>z_breach</b> attribute can be used as an aid in creating a hydrologically sound DEM. However, mind the multiple, identically located points at the connections between the input <em>vect</em> lines, and the multiple points per one input <em>rast</em> cell in case when 2 or more <em>vect</em> lines flow through a given cell.
+<p>
+In the 1st case, as the <b>z_breach</b> attribute of the <b>duplicates</b> is <b>identical</b> anyway, you might want to have only one of the points (eg. if your interpolation program can't handle duplicate input points). To remove duplicates use <em>v.clean tool=rmdupl</em>. To assess the latter case, you will need to modify your input <em>vect</em> so that max one stream flows through any <em>rast</em> cell, or increase the resolution of your input DEM so that it meets the level of detail in your input drainage network.
+<p>
+
+<h2>NOTES</h2>
+
+<ul>
+ 1. <em>v.breach</em> is a Bash and Awk script. Requires GRASS 6.x.<br>
+ <p>
+ 2. The script was developed for a <b>metric</b> coordinate system, but I suppose it should work in a feet (or other unit) based system as well, only it *might* require a few modifications in parts where it directly refers to meters. It won't work in lat/lon locations though for sure, without serious modifications. Summarising, candidate troublemakers in the script might be:
+ <p>
+ r.buffer distances= <br>
+ v.parallel distance= <br>
+ v.distance <br>
+ v.to.db option=length <br>
+ v.to.db option=coor <br>
+ v.segment <br>
+ sort -n <br>
+ all the Awk maths <br>
+ <p>
+ 3. Each input <em>vect</em> line must have a <b>unique</b> category in <b>layer 1</b>. Use <em>v.category</em> to add categories, if needed.
+ <p>
+ 4. An input line must not cross itself. Otherwise the script will go into an infinite loop.
+ <p>
+ 5. For the same reason, any input lines must not constitute a loop.
+ <p>
+ 6. Any of the input <em>vect</em> lines must not stand out of the input <em>rast</em> DEM.
+ <p>
+ 7. The input DEM must have rectangular cells.
+<p>
+ 8. Input lines direction is the key. Use <em>v.flip</em> from GRASS AddOns to flip any if needed.
+</ul>
+
+<h2>SEE ALSO</h2>
+
+<em>
+ <a href="r.carve.html">r.carve</a>,
+ <a href="v.category.html">v.category</a>,
+ <a href="v.clean.html">v.clean</a>
+</em>
+
+<h2>AUTHOR</h2>
+Maciej Sieczka
+
+<p>
+<i>Last changed: $Date$</i>
Property changes on: grass-addons/grass6/vector/v.breach/description.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass6/vector/v.breach/v.breach
===================================================================
--- grass-addons/grass6/vector/v.breach/v.breach (rev 0)
+++ grass-addons/grass6/vector/v.breach/v.breach 2013-01-20 18:17:17 UTC (rev 54721)
@@ -0,0 +1,550 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: v.breach
+#
+# AUTHOR(S): Maciej Sieczka
+#
+# PURPOSE: Create vector maps of lines and points of continously lowering
+# elevation down the input watercourses.
+#
+# VERSION: 5.9.7, developed over GRASS 6.3 CVS 2007.01.24, SVN 2008.01.22
+#
+# COPYRIGHT: (c) 2007, 2008 Maciej Sieczka
+#
+# NOTES: This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+# CHANGELOG:
+#
+# 5.9.6: first public release
+# 5.9.7: workaround for a "g.region res= vect= -a" given all together might not
+# work as expectected, depending on the uinitial region settings
+
+#%Module
+#% description: Create vector maps of lines and points of continously lowering elevation down the input watercourses.
+#%END
+
+#%option
+#% key: vect
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Input watercourses vector
+#% required : yes
+#%END
+
+#%option
+#% key: rast
+#% gisprompt: old,cell,raster
+#% type: string
+#% description: Input DEM raster
+#% required : yes
+#%END
+
+#%option
+#% key: out_ln
+#% type: string
+#% gisprompt: new,dig,vector
+#% description: Output watercourses vector
+#% required : yes
+#%END
+
+#%option
+#% key: out_pt
+#% type: string
+#% gisprompt: new,dig,vector
+#% description: Output points vector
+#% required : yes
+#%END
+
+#%option
+#% key: depth
+#% type: double
+#% answer: 0
+#% description: Additional stream depth (meters)
+#% required : yes
+#%END
+
+##%option
+##% key: output_raster
+##% gisprompt: new,cell,raster
+##% type: string
+##% description: Output breached DEM raster
+##% required : no
+##%END
+
+echo "This script needs to be updated to
+exit
+
+# called from Grass?
+if test "$GISBASE" = ""; then
+ echo "ERROR: 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
+ echo "ERROR: awk required, please install awk or gawk first." 1>&2
+ exit 1
+fi
+
+# set environment so that awk works properly in all languages
+unset LC_ALL
+LC_NUMERIC=C
+export LC_NUMERIC
+
+INL="$GIS_OPT_VECT"
+INLNOAT=`echo "$GIS_OPT_VECT" | sed 's/@.*//'`
+INR="$GIS_OPT_RAST"
+OUTP="$GIS_OPT_OUT_PT"
+OUTL="$GIS_OPT_OUT_LN"
+DPTH="$GIS_OPT_DEPTH"
+#OUTR="$GIS_OPT_OUTPUT_RASTER"
+
+# first of all check DEM cells proportions
+eval `r.info -s map="${INR}"`
+if [ "$ewres" -ne "$nsres" ] ; then
+ echo "ERROR: The input DEM must have rectangular cells." 1>&2
+ exit 1
+fi
+
+eval `g.gisenv`
+# check if output vectors exist
+for i in "${OUTP}" "${OUTL}"; do
+ g.findfile elem=vector file="${i}" mapset="${MAPSET}" > /dev/null
+ if [ $? -eq 0 ] ; then
+ echo "ERROR: The output vector <"${i}"> already exists in current mapset." 1>&2
+ exit 1
+ fi
+done
+
+# all above OK - proceed: set up temporary files
+TMP="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMP" ] ; then
+ echo "ERROR: Unable to create temporary files." 1>&2
+ exit 1
+fi
+
+PROG=`basename $0 | sed 's/\./_/g'`
+UNQ=`basename $TMP | sed 's/\./_/g'`"_${PROG}"
+
+# check if intermediate vectors exist
+for i in "${INLNOAT}_tmp_${UNQ}"\
+ "${OUTP}_grid_${UNQ}"\
+ "${OUTP}_grid_addcat_${UNQ}"\
+ "${OUTP}_grid_addcat_bnd_patch_break_ln_2l_tmp_${UNQ}"\
+ "${OUTP}_grid_addcat_bnd_patch_break_ln_tmp_${UNQ}"\
+ "${OUTP}_grid_addcat_bnd_patch_break_tmp_${UNQ}"\
+ "${OUTP}_grid_addcat_bnd_patch_tmp_${UNQ}"\
+ "${OUTP}_grid_addcat_bnd_tmp_${UNQ}"\
+ "${OUTP}_segm_2l_l_center_tmp_${UNQ}"\
+ "${OUTP}_segm_2l_l_tmp_${UNQ}"\
+ "${OUTP}_segm_2l_r_center_tmp_${UNQ}"\
+ "${OUTP}_segm_2l_r_tmp_${UNQ}"\
+ "${OUTP}_segm_2l_shorties_tmp_${UNQ}"\
+ "${OUTP}_shorties_break_onlycat_prune_poly_prune_${UNQ}"\
+ "${OUTP}_shorties_break_onlycat_prune_poly_prune_addcat_tmp_${UNQ}"\
+ "${OUTP}_shorties_break_onlycat_prune_poly_prune_tmp_${UNQ}"\
+ "${OUTP}_shorties_break_onlycat_prune_poly_tmp_${UNQ}"\
+ "${OUTP}_shorties_break_onlycat_prune_tmp_${UNQ}"\
+ "${OUTP}_shorties_break_onlycat_tmp_${UNQ}"\
+ "${OUTP}_shorties_break_tmp_${UNQ}"\
+ "${OUTP}_shorties_tmp_${UNQ}"\
+ "${OUTP}_tmp_${UNQ}"; do
+ g.findfile elem=vector file="${i}" mapset="${MAPSET}" > /dev/null
+ if [ $? -eq 0 ] ; then
+ echo "ERROR: The intermediate vector <"${i}"> already exists in current mapset." 1>&2
+ exit 1
+ fi
+done
+
+# check if intermediate rasters exist
+for i in "${OUTP}_rast_${UNQ}"\
+ "${OUTP}_rast_buf_${UNQ}"\
+ "${OUTP}_rast_buf_uniq_${UNQ}"; do
+ g.findfile elem=cell file="${i}" mapset="${MAPSET}" > /dev/null
+ if [ $? -eq 0 ] ; then
+ echo "ERROR: The intermediate raster <"${i}"> already exists in current mapset." 1>&2
+ exit 1
+ fi
+done
+
+# check if intermediate region def exists
+i="region_${UNQ}"
+g.findfile elem=windows file="${i}" mapset="${MAPSET}" > /dev/null
+ if [ $? -eq 0 ] ; then
+ echo "ERROR: The intermediate region definition <"${i}"> already exists in current mapset." 1>&2
+ exit 1
+ fi
+
+
+# define the cleanup procedure
+cleanup()
+{
+ \rm -f $TMP
+# \rm -f $TMP.${PROG}
+ \rm -f $TMP.${PROG}.lcats
+ \rm -f $TMP.${PROG}.offsets
+ \rm -f $TMP.${PROG}.pruned
+ \rm -f $TMP.${PROG}.segm
+ \rm -f $TMP.${PROG}.raw
+ \rm -f $TMP.${PROG}.input
+ \rm -f $TMP.${PROG}.lower
+ \rm -f $TMP.${PROG}.lower.sortbynodes
+ \rm -f $TMP.${PROG}.common
+ \rm -f $TMP.${PROG}.sql
+ \rm -f $TMP.${PROG}.reclassit
+ \rm -f $TMP.${PROG}.outl
+ \rm -f $TMP.${PROG}.outl_raw
+ \rm -f $TMP.${PROG}.outl_raw_cor
+
+g.remove rast="${OUTP}_rast_${UNQ}","${OUTP}_rast_buf_${UNQ}","${OUTP}_rast_buf_uniq_${UNQ}" > /dev/null
+g.remove region="region_${UNQ}" > /dev/null
+g.remove vect="${OUTP}_grid_${UNQ}",\
+"${OUTP}_grid_addcat_${UNQ}",\
+"${OUTP}_grid_addcat_bnd_patch_break_ln_2l_tmp_${UNQ}",\
+"${OUTP}_grid_addcat_bnd_patch_break_ln_tmp_${UNQ}",\
+"${OUTP}_grid_addcat_bnd_patch_break_tmp_${UNQ}",\
+"${OUTP}_grid_addcat_bnd_patch_tmp_${UNQ}",\
+"${OUTP}_grid_addcat_bnd_tmp_${UNQ}",\
+"${OUTP}_segm_2l_l_center_tmp_${UNQ}",\
+"${OUTP}_segm_2l_l_tmp_${UNQ}",\
+"${OUTP}_segm_2l_r_center_tmp_${UNQ}",\
+"${OUTP}_segm_2l_r_tmp_${UNQ}",\
+"${OUTP}_segm_2l_shorties_tmp_${UNQ}",\
+"${OUTP}_shorties_break_onlycat_prune_poly_prune_${UNQ}",\
+"${OUTP}_shorties_break_onlycat_prune_poly_prune_addcat_tmp_${UNQ}",\
+"${OUTP}_shorties_break_onlycat_prune_poly_prune_tmp_${UNQ}",\
+"${OUTP}_shorties_break_onlycat_prune_poly_tmp_${UNQ}",\
+"${OUTP}_shorties_break_onlycat_prune_tmp_${UNQ}",\
+"${OUTP}_shorties_break_onlycat_tmp_${UNQ}",\
+"${OUTP}_shorties_break_tmp_${UNQ}",\
+"${OUTP}_shorties_tmp_${UNQ}",\
+"${OUTP}_tmp_${UNQ}",\
+"${INLNOAT}_tmp_${UNQ}" > /dev/null
+}
+
+# what to do in case of user break:
+exitprocedure()
+{
+ echo "User break!" 1>&2
+ cleanup
+ exit 1
+}
+# shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 15
+
+
+
+### DO IT ###
+
+# set the region to match whole input vector but respect current resolution
+res=$nsres
+g.region vect="${INL}" res=$res save="region_${UNQ}" -a -u > /dev/null
+# we need to repeat this for there is a bug/feature that setting bot res=, -a
+# and vect= might not work as expected
+g.region region="region_${UNQ}" vect="${INL}" res=$res save="region_${UNQ}" -a -u --o > /dev/null
+
+WIND_OVERRIDE="region_${UNQ}"
+export WIND_OVERRIDE
+
+eval `g.region -g` > /dev/null
+
+# make a raster that will encompass all the input lines; each cell must have an
+# unique value
+
+buf=`echo $res | awk '{printf "%.8f\n",$0*2}'`
+
+v.to.rast input="${INL}" output="${OUTP}_rast_${UNQ}" use=val layer=1 value=1 > /dev/null
+r.buffer input="${OUTP}_rast_${UNQ}" output="${OUTP}_rast_buf_${UNQ}" distances=$buf units=meters > /dev/null
+r.mapcalc ""${OUTP}_rast_buf_uniq_${UNQ}"=if(not(isnull("${OUTP}_rast_buf_${UNQ}")),row()*col(),null())" > /dev/null
+
+# transform it into vector areas - we will use their boundaries
+r.to.vect input="${OUTP}_rast_buf_uniq_${UNQ}" output="${OUTP}_grid_${UNQ}" feature=area > /dev/null
+
+# add categories to bnds so we can select them
+v.category input="${OUTP}_grid_${UNQ}" output="${OUTP}_grid_addcat_${UNQ}" type=boundary option=add cat=1 layer=1 step=1 > /dev/null
+
+# get the input lines categories
+#lcats=`v.category input="${INL}" option=print`
+v.category input="${INL}" option=print > $TMP.${PROG}.lcats
+
+
+ # THE LOOP PER EACH LINE CAT #
+
+#for i in $lcats; do
+cat $TMP.${PROG}.lcats | while read i; do
+
+ v.extract -t input="${INL}" output="${INLNOAT}_tmp_${UNQ}" layer=1 list=$i type=line new=-1 --o > /dev/null
+
+ # select only boundaries overlaping with the line of the current CAT
+ v.select -t ainput="${OUTP}_grid_addcat_${UNQ}" atype=boundary alayer=1 binput="${INLNOAT}_tmp_${UNQ}" btype=line blayer=1 output="${OUTP}_grid_addcat_bnd_tmp_${UNQ}" operator=overlap --o > /dev/null
+
+ # patch input with selected grid-boundaries
+ v.patch input="${INLNOAT}_tmp_${UNQ}","${OUTP}_grid_addcat_bnd_tmp_${UNQ}" output="${OUTP}_grid_addcat_bnd_patch_tmp_${UNQ}" --o > /dev/null
+
+ # divide lines into segments using that - 1 line segment per 1 cell
+ v.clean input="${OUTP}_grid_addcat_bnd_patch_tmp_${UNQ}" output="${OUTP}_grid_addcat_bnd_patch_break_tmp_${UNQ}" type=line,boundary tool=break --o > /dev/null
+
+ # extract lines
+ v.extract -t input="${OUTP}_grid_addcat_bnd_patch_break_tmp_${UNQ}" output="${OUTP}_grid_addcat_bnd_patch_break_ln_tmp_${UNQ}" type=line new=-1 --o > /dev/null
+
+ # add cats to layer 2; save each line into a separate file as we'll
+ # them later (in the loop at end of the script)
+ v.category input="${OUTP}_grid_addcat_bnd_patch_break_ln_tmp_${UNQ}" output="${OUTP}_grid_addcat_bnd_patch_break_ln_2l_tmp_${UNQ}" type=line option=add layer=2 --o > /dev/null
+
+ # prepare the input for v.segment that will extract center points of
+ # each line segment; excluding the first and the last offset
+ v.to.db -p map="${OUTP}_grid_addcat_bnd_patch_break_ln_2l_tmp_${UNQ}" layer=2 type=line option=length units=me column=dummy_vtodb | awk -F "|" 'NR>1 {printf "%s","P "$1" "$1" "; printf "%.16f\n",$2/2}' > $TMP.${PROG}.offsets
+
+ # create parallel lines on the sides of the original one
+ v.parallel input="${OUTP}_grid_addcat_bnd_patch_break_ln_2l_tmp_${UNQ}" output="${OUTP}_segm_2l_r_tmp_${UNQ}" distance=0.00001 --o > /dev/null
+
+ v.parallel input="${OUTP}_grid_addcat_bnd_patch_break_ln_2l_tmp_${UNQ}" output="${OUTP}_segm_2l_l_tmp_${UNQ}" distance=-0.00001 --o > /dev/null
+
+ # Pipe offsets into v.segment. The output are points located exactly in the
+ # middle of the segments.
+ cat $TMP.${PROG}.offsets | v.segment llayer=2 input="${OUTP}_segm_2l_r_tmp_${UNQ}" output="${OUTP}_segm_2l_r_center_tmp_${UNQ}" --o > /dev/null
+
+ cat $TMP.${PROG}.offsets | v.segment llayer=2 input="${OUTP}_segm_2l_l_tmp_${UNQ}" output="${OUTP}_segm_2l_l_center_tmp_${UNQ}" --o > /dev/null
+
+ # connect the closest points - creates mega-short lines going across the
+ # input lines
+ v.distance -p from="${OUTP}_segm_2l_r_center_tmp_${UNQ}" to="${OUTP}_segm_2l_l_center_tmp_${UNQ}" from_type=point to_type=point from_layer=1 to_layer=1 output="${OUTP}_segm_2l_shorties_tmp_${UNQ}" dmax=-1 upload=cat column=dummy_vdistance --o > /dev/null
+
+ # patch the input and shorties...
+ v.patch input="${OUTP}_segm_2l_shorties_tmp_${UNQ}","${INLNOAT}_tmp_${UNQ}" output="${OUTP}_shorties_tmp_${UNQ}" --o > /dev/null
+
+ # ...break that to create the nodes...
+ v.clean input="${OUTP}_shorties_tmp_${UNQ}" output="${OUTP}_shorties_break_tmp_${UNQ}" type=line tool=break --o > /dev/null
+
+ # ...and extract only lines with cats - shorties will be excluded, nodes
+ # will remain
+ v.extract -t input="${OUTP}_shorties_break_tmp_${UNQ}" output="${OUTP}_shorties_break_onlycat_tmp_${UNQ}" type=line layer=1 new=0 --o > /dev/null
+
+# It would be cool to extract the nodes in order-down-the-line-direction with
+# v.to.points now, but it fails (the order is not preserved). Workaround by
+# removing all the vertices, turning remaining nodes into vertices and
+# v.to.points these - the resulting points categories will have cats growing
+# down the flowpath as required for further processing.
+
+# remove all the vertices - nodes will remain...
+
+###v.clean input="${OUTP}_shorties_break_onlycat" output="${OUTP}_shorties_break_onlycat_prune" type=line tool=prune thresh=???!!!
+
+# yes, it would be nice to be able just to prune but the thresh is way too funky
+# for me to comprehend it; v.out.ascii | grep+awk and back instead:
+
+ v.out.ascii "${OUTP}_shorties_break_onlycat_tmp_${UNQ}" format=standard | awk 'NR>10' | grep -C1 "L \| 1 " | awk '/L / {$2=" 2"} /--/ {next} {print}' > $TMP.${PROG}.pruned
+
+ v.in.ascii -n format=standard input=$TMP.${PROG}.pruned output="${OUTP}_shorties_break_onlycat_prune_tmp_${UNQ}" --o
+
+ # ...so we can turn them into vertices now...
+ v.build.polylines input="${OUTP}_shorties_break_onlycat_prune_tmp_${UNQ}" output="${OUTP}_shorties_break_onlycat_prune_poly_tmp_${UNQ}" --o > /dev/null
+
+ # but, there is a bug in v.build.polylines (BT#4247); workaround - prune
+ # *doubled* vertices it creates (v.clean tool=prune works OK for thresh=0
+ # at least)
+
+ # UPDATE, 12.03.2007: this bug is supposedly fixed in 6.3 CVS few days
+ # ago. I haven't had time to check it yet. Anyway, the fix was not
+ # backported to 6.2, thus the workaround is still needed; it won't do
+ # any harm in 6.3 CVS, only adds to processing time.
+
+ v.clean input="${OUTP}_shorties_break_onlycat_prune_poly_tmp_${UNQ}" output="${OUTP}_shorties_break_onlycat_prune_poly_prune_tmp_${UNQ}" type=line tool=prune thresh=0 --o > /dev/null
+
+ # add a category back (to keep the original one - v.build.polylines
+ # removes any cats)
+ v.category input="${OUTP}_shorties_break_onlycat_prune_poly_prune_tmp_${UNQ}" output="${OUTP}_shorties_break_onlycat_prune_poly_prune_addcat_tmp_${UNQ}" type=line option=add cat=${i} --o > /dev/null
+
+ # add the result into ASCII vector (excluding cats for layer 2 - they
+ # will have to be created from scratch, unique per each segment -
+ # currently there doubled cats due to cats for layer 2 have been added
+ # to each single line; all the lines are stored there we'll import it
+ # when all done
+ v.out.ascii input="${OUTP}_shorties_break_onlycat_prune_poly_prune_addcat_tmp_${UNQ}" format=standard | awk 'NR>10' | awk 'NR==1 {$2=$2-2} NR==3 {next} {print}' | tac | awk 'NR==3 {next} {print}' | tac >> $TMP.${PROG}.segm
+
+done
+
+# phew, OK, now:
+
+# Create the lines vector for extracting points the middles of the segments
+# created above. Those points will be used for breaching DEM down the flowpath.
+v.in.ascii -n format=standard input=$TMP.${PROG}.segm output="${OUTP}_shorties_break_onlycat_prune_poly_prune_${UNQ}" > /dev/null
+
+# extract the line vertices and start/end nodes
+v.to.points -v input="${OUTP}_shorties_break_onlycat_prune_poly_prune_${UNQ}" type=line output="${OUTP}" llayer=1 > /dev/null
+
+# add columns to store the XY and the Z sampled from DEM, and later, the breached Z
+v.db.addcol map="${OUTP}" layer=2 'columns=x double, y double, z double, z_breach double' > /dev/null
+
+# upload XY
+v.to.db map="${OUTP}" type=point layer=2 option=coor column=x,y,z > /dev/null
+
+# sample Z from input DEM
+v.what.rast vector="${OUTP}" raster="${INR}" layer=2 column=z > /dev/null
+
+# print to text file
+v.db.select -c map="${OUTP}" layer=2 column=cat,lcat,x,y,z fs=" " | awk '{printf "%s %s %.6f %.6f %.6f\n",$1,$2,$3,$4,$5}' > $TMP.${PROG}.raw
+
+
+
+### DO THE MAGIC ###
+
+lower()
+{
+# set each following "Z" to previous, if it is higher
+awk '
+{
+if (NR==1)
+
+ {along=$2; first=$5; printf $1" "$2" "$3" "$4" "; printf "%.6f\n",$5}
+
+else
+ if ($5>=first && $2==along)
+
+ {($5=first-0.000001); first=$5; printf "%s %s %.6f %.6f %.6f\n",$1,$2,$3,$4,$5}
+
+ else
+
+ {along=$2; first=$5; printf "%s %s %.6f %.6f %.6f\n",$1,$2,$3,$4,$5}
+}' $TMP.${PROG}.input > $TMP.${PROG}.lower
+}
+
+common()
+{
+# sort the lowered points by XY, then these by Z - in order to group the
+# identical points (common nodes of more lines)...
+sort -n -k 3,4 -k 5 $TMP.${PROG}.lower > $TMP.${PROG}.lower.sortbynodes
+# ...and set the "z" to the lowest value of all identical nodes
+awk '
+{
+if (NR==1)
+
+ {x=$3; y=$4; z=$5; printf "%s %s %.6f %.6f %.6f\n",$1,$2,$3,$4,$5}
+
+else
+ if ($3==x && $4==y)
+
+ {$5=z; printf "%s %s %.6f %.6f %.6f\n",$1,$2,$3,$4,$5}
+
+ else
+
+ {x=$3; y=$4; z=$5; printf "%s %s %.6f %.6f %.6f\n",$1,$2,$3,$4,$5}
+}' $TMP.${PROG}.lower.sortbynodes | sort -n -k 1,1 > $TMP.${PROG}.common
+}
+
+cp $TMP.${PROG}.raw $TMP.${PROG}.input
+lower
+common
+
+# compare; if identical we are done, otherwise repeat
+cmp $TMP.${PROG}.common $TMP.${PROG}.input > /dev/null
+
+while [ "$?" -eq "1" ]; do
+
+cp $TMP.${PROG}.common $TMP.${PROG}.input
+lower
+common
+cmp $TMP.${PROG}.common $TMP.${PROG}.input > /dev/null
+
+done
+
+# upload breached Z into the table
+
+tbl=`v.db.connect -g "${OUTP}" | awk '{print $2}'`
+dbs=`v.db.connect -g "${OUTP}" | awk '{print $4}'`
+drv=`v.db.connect -g "${OUTP}" | awk '{print $5}'`
+
+# upload SQL commands into txt file first = ~20 times faster db.execute
+for z in `awk -v DPTH="$DPTH" '{printf "%f ",$5-DPTH}' $TMP.${PROG}.input`; do
+ row=`expr $row + 1`
+ echo "UPDATE $tbl SET z_breach=$z WHERE cat=$row;" >> $TMP.${PROG}.sql
+done
+
+db.execute input=$TMP.${PROG}.sql database=$dbs driver=$drv > /dev/null
+
+# Transfer 'z' and 'z_breach' from points to $OUTL's lines' segments. v.distance
+# does only anytype->point transfers, so a workaround is needed. Additionaly it
+# is complicated that at each lines' connection there as many points as lines
+# connected there. Thus each line has to be processed separately.
+
+ # ANOTHER LOOP PER EACH LINE CAT #
+
+cat $TMP.${PROG}.lcats | while read i; do
+
+### LINES
+
+ v.extract -t input="${INL}" output="${INLNOAT}_tmp_${UNQ}" layer=1 list=$i type=line new=-1 --o > /dev/null
+
+ # select only boundaries overlaping with the line of the current CAT
+ v.select -t ainput="${OUTP}_grid_addcat_${UNQ}" atype=boundary alayer=1 binput="${INLNOAT}_tmp_${UNQ}" btype=line blayer=1 output="${OUTP}_grid_addcat_bnd_tmp_${UNQ}" operator=overlap --o > /dev/null
+
+ # patch input with selected grid-boundaries
+ v.patch input="${INLNOAT}_tmp_${UNQ}","${OUTP}_grid_addcat_bnd_tmp_${UNQ}" output="${OUTP}_grid_addcat_bnd_patch_tmp_${UNQ}" --o > /dev/null
+
+ # divide lines into segments using that - 1 line segment per 1 cell
+ v.clean input="${OUTP}_grid_addcat_bnd_patch_tmp_${UNQ}" output="${OUTP}_grid_addcat_bnd_patch_break_tmp_${UNQ}" type=line,boundary tool=break --o > /dev/null
+
+ # extract lines
+ v.extract -t input="${OUTP}_grid_addcat_bnd_patch_break_tmp_${UNQ}" output="${OUTP}_grid_addcat_bnd_patch_break_ln_tmp_${UNQ}" type=line new=-1 --o > /dev/null
+
+ # add cats to layer 2; save each line into a separate file as we'll
+ # them later (in the loop at end of the script)
+ v.category input="${OUTP}_grid_addcat_bnd_patch_break_ln_tmp_${UNQ}" output="${OUTP}_grid_addcat_bnd_patch_break_ln_2l_tmp_${UNQ}" type=line option=add layer=2 --o > /dev/null
+
+
+### POINTS
+
+ # extract points of given lcat to a separate file
+ v.extract -t input="${OUTP}" output="${OUTP}_tmp_${UNQ}" type=point layer=2 new=-1 where=lcat="${i}" --o > /dev/null
+
+ # find points and lines' segments closest to each other
+ v.distance -p from="${OUTP}_tmp_${UNQ}" to="${OUTP}_grid_addcat_bnd_patch_break_ln_2l_tmp_${UNQ}" from_type=point to_type=line from_layer=2 to_layer=2 dmax=-1 upload=cat column=dummy | awk 'NR>1' | sed 's/|/ /' > $TMP.${PROG}.reclassit
+
+ # add to ASCII vector file...
+ v.out.ascii format=standard in="${OUTP}_grid_addcat_bnd_patch_break_ln_2l_tmp_${UNQ}" | awk 'NR>10' > $TMP.${PROG}.outl_raw
+
+ cat $TMP.${PROG}.reclassit | while read j; do
+
+ is=`echo $j | awk '{print $2}'`
+ to=`echo $j | awk '{print $1}'`
+
+ # ...reclassing categories as needed (with sed), adding "skrzat"
+ # markers to avoid multiple substitutions
+ sed "s/^ 2 $is /skrzat 2 $to/" $TMP.${PROG}.outl_raw > $TMP.${PROG}.outl_raw_cor
+ cp $TMP.${PROG}.outl_raw_cor $TMP.${PROG}.outl_raw
+ done
+
+ # remove "skrzat" markers
+ sed "s/skrzat//" $TMP.${PROG}.outl_raw > $TMP.${PROG}.outl_raw_cor
+
+ # patch lines together into an ASCII vect - this will be the $OUTL module's output
+ cat $TMP.${PROG}.outl_raw_cor >> $TMP.${PROG}.outl
+
+done
+
+# import segments
+v.in.ascii -n format=standard input=$TMP.${PROG}.outl output="${OUTL}" > /dev/null
+
+# now that categories in points and lines outputs are the same, we can copy the
+# table from points for lines
+db.copy from_driver=$drv from_database=$dbs to_driver=$drv to_database=$dbs to_table="${OUTL}_2" select="SELECT cat,lcat,z,z_breach FROM $tbl" > /dev/null
+
+v.db.connect map="${OUTL}" layer=2 driver=$drv database=$dbs table="${OUTL}_2" key=cat > /dev/null
+
+### ALL DONE ###
+
+cleanup
+
+echo 1>&2
+echo "Done." 1>&2
+echo 1>&2
Property changes on: grass-addons/grass6/vector/v.breach/v.breach
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/x-sh
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list