#!/bin/sh
export GIS_LOCK=$$
PATH=$PATH:/grass/scripts
GISDBASE="C:/GRASSDATA"
OUTPATH="C:/GRASSDATA/Scripts/LRS/"
MAPSET="LRS"
LOCATION_NAME="latlon"
OUTPUT_PATH=$OUTPATH"output/"$MAPSET"/"
g.mapset mapset=$MAPSET location=$LOCATION_NAME gisdbase=$GISDBASE
db.connect driver=sqlite database=$GISDBASE/$LOCATION_NAME/$MAPSET/sqlite.db
################# PREPARE OUTPUT FOLDER #################################
echo "Checking for top-level output folder"
  # if output directory doesn't exist, make it
if [ ! -d "$OUTPATH/output" ]; then
  echo "Creating top-level output folder"
  mkdir $OUTPATH/output
fi

echo "Checking for output/"$MAPSET" folder"
  # if output/$MAPSET directory exists, delete it
if [ -d "$OUTPATH/output/"$MAPSET ]; then
  echo "Removing output/"$MAPSET" folder"
  rm -r $OUTPATH/output/$MAPSET
fi

  # create output/$MAPSET directory
echo "Creating output/"$MAPSET" folder"
mkdir $OUTPATH/output/$MAPSET
#=========================================  CREATE STREAM NETWORK SECTION ============================================
g.region -a region=watershed
BDY="12462500"
r.mask -o input=tmp_BDY_rast maskcats=$BDY

r.watershed -m --overwrite elevation=DEM_mask accumulation=tmp_accum drainage=tmp_flowdir basin=basin stream=stream threshold=2000 memory=300
  # vector directions from r.to.vect are all over the place. r.stream.extract creates a vector file
  # which is SOOOO much better; vector direction is all downstream. Note the IDs are different from 
  # r.watershed. You cannot use these stream maps in comparison to maps created by the other scripts.
r.stream.extract --overwrite elevation=DEM_mask threshold=1000 stream_vect=stream_ext
  # note that running v.clean, v.generalize, and absolutley v.build.polylines messes up 
  # the vector directions and the network becomes useless really. Upstream and downstream
  # run into each other. Also the v.clean and v.generalize mess up some of the junctions.
  #
  # nothing can be easy, can it?  Ok - the arcs added to the network linking the gauges to the stream
  # have their direction backward ... you'll never get a gauge-to-gauge path as one gauge will
  # always have a downstream arc and the other will have an upstream one thus getting the route
  # shut off.  Need to snap the points onto the vector line BEFORE creating the network.
  # First need to make a "clean" version of the stream network so other intersections do not come into play.
v.clean --overwrite input=stream_ext output=stream_clean tool=break #,rmdupl,rmdangle,rmline
  # use v.distance to create the vectors connecting from the gauge locations perpendicular to the stream vector
v.distance --overwrite from=gauges from_type=point to=stream_clean to_type=line output=connectors upload=dist,cat column="dist,to_cat"
  # patch the stream vector and the gauge connector vector to get intersections that are not "cleaned"
v.patch input=stream_clean,connectors out=lines_merged --overwrite 
  # run clean with the err option to generate the intersection points
v.clean input=lines_merged output=snap_points error=crosspoints_err tool=break --overwrite
  # there should be 15 intersection points, but for some reason there are only 11 in the error output file...
  # does anyone now what is going on? Is there a threshold for creating the error? Some missing points
  # are less than 2 meters from the stream and some are 260 meters ... very strange.