#!/bin/sh
############################################################################
#
# MODULE:       v.ldm
# AUTHOR(S):    Alexander Muriy
#               (Institute of Environmental Geoscience, Moscow, Russia)  
#               e-mail: amuriy AT gmail DOT com 
#
# PURPOSE:      Compute "Linear Directional Mean" of vector lines, display arrow 
#               on the graphic monitor, save to vector line and update attribute table
#               with LDM parameters.
#
# COPYRIGHT:    (C) 2011 Alexander Muriy / GRASS Development Team
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
############################################################################
#%Module 
#%  description: Compute "Linear Directional Mean" of vector lines, displays arrow on the graphic monitor, save to vector line and update attribute table with LDM parameters. 
#%  keywords: display, graphics, vector, symbology
#%End
#%Option
#%  key: map
#%  type: string
#%  required: yes
#%  multiple: no
#%  key_desc: name
#%  description: Name of input vector map
#%  gisprompt: old,vector,vector
#%End
#%Option
#%  key: layer
#%  type: integer
#%  required: no
#%  description: Layer number 
#%  answer: 1
#%End
#%Option
#%  key: output
#%  type: string
#%  required: no
#%  multiple: no
#%  key_desc: name
#%  description: Name of output vector map
#%  gisprompt: new,vector,vector
#%End
#%Option
#%  key: width
#%  type: integer
#%  required: no
#%  description: Width of LDM line with arrow
#%  answer: 5
#%End
#%Option
#%  key: color
#%  type: string
#%  required: no
#%  description: Color for LDM line with arrow (standard color name or R:G:B)
#%  answer: black
#%  gisprompt: old_color,color,color
#%End
#%Option
#%  key: graph
#%  type: string
#%  required: no
#%  multiple: no
#%  key_desc: name
#%  description: Output file to save Linear Directional Mean graphics
#%  gisprompt: new_file,file,output
#%End
#%Flag
#%  key: n
#%  description: Don't show LDM graphics on the screen
#%End
#%Flag
#%  key: s 
#%  description: Save LDM graphics to file (for use with d.graph)
#%End
#%Flag
#%  key: l
#%  description: Save LDM line as vector map and create attribute table with LDM parameters
#%End
#%Flag
#%  key: g
#%  description: Print in shell script style
#%End


if [ -z "$GISBASE" ] ; then
    echo "You must be in GRASS GIS to run this program." 1>&2
    exit 1
fi

if [ "$1" != "@ARGS_PARSED@" ] ; then
    exec g.parser "$0" "$@"
fi


## set environment so that awk works properly in all languages ##
unset LC_ALL
export LC_NUMERIC=C


############################################################
cleanup()
{ 
    g.mremove -f vect=V_LDM_HULL --q
    rm -f $TMP1*
    
}
############################################################
## what to do in case of user break:
exitprocedure()
{
    echo "User break!"
    cleanup
    exit 1
}

## shell check for user break (signal list: trap -l)
trap "exitprocedure" 2 3 15

############################################################
## check for awk
if [ ! -x "$(which awk)" ] ; then
    g.message -e "<awk> required, please install <awk> or <gawk> first"
    exit 1
fi

## check for cs2cs
if [ ! -x "$(which cs2cs)" ] ; then
    g.message -e "cs2cs required, please install <cs2cs> package first"
    exit 1
fi

## does map exist?
eval "$(g.gisenv)"
eval "$(g.findfile mapset="$MAPSET" element=vector file="$GIS_OPT_MAP")"
if [ ! "$file" ]; then
    g.message -e "Vector map <"$GIS_OPT_MAP"> not found"
    exit 1
fi

## check for lines in input map
eval "$(v.info -t "$GIS_OPT_MAP")"

if [ "$lines" -eq 0 ]; then
    g.message -e "Input vector map not contain lines"
    exit 1
fi

# for vars in "$(v.info -t "$GIS_OPT_MAP" | cut -d"=" -f1)"; do
#     unset -v "$vars"
# done

############################################################
## setup temporary files ##

# create temporary file
TMP1="$(g.tempfile pid=$$)"
if [ "$?" -ne 0 ] || [ -z "$TMP1" ] ; then
    g.message -e "ERROR: Unable to create temporary files."
    exit 1
fi

############################################################
## DO IT

## compute LDM characteristics and print them

## make convex hull around lines
v.hull -a in=$GIS_OPT_MAP out=V_LDM_HULL --q > /dev/null 2>&1

CENTER_COORDS=$(v.to.db -p V_LDM_HULL opt=coor type=centroid --q | cut -d'|' -f2,3 | tr '|' ' ')
CENTER_X=$(echo $CENTER_COORDS | cut -d' ' -f1 | awk '{printf "%0.6f\n",$1}')
CENTER_Y=$(echo $CENTER_COORDS | cut -d' ' -f2 | awk '{printf "%0.6f\n",$1}')


# CENTER_COORDS=$(g.region -u vect=$GIS_OPT_MAP -cg | cut -d"=" -f2 | tr '\n' ' ')
# CENTER_X=$(echo ${CENTER_COORDS} | cut -d" " -f1)
# CENTER_Y=$(echo ${CENTER_COORDS} | cut -d" " -f2)


if [ "$GIS_FLAG_G" -eq 1 ]; then
    echo "mean_center=$(echo $CENTER_COORDS | tr ' ' ',')"
else 
    echo "Mean Center: $(echo $CENTER_COORDS | tr ' ' ',')"
fi


END_AZIM=$(v.to.db -p $GIS_OPT_MAP layer=$GIS_OPT_LAYER opt=azimuth units=radians --q | awk -F"|" -v count=${VECT_COUNT} '{
s1=sin($2); c1=cos($2); 
print s1,c1}' | awk '{
s2 += $1; c2 += $2} 
END {a1=atan2(s2,c2); 
a2=a1*57.2958; 
print a2}')

COMPASS_ANGLE=$(v.to.db -p $GIS_OPT_MAP layer=$GIS_OPT_LAYER opt=azimuth units=radians --q | awk -F"|" -v count=${VECT_COUNT} '{
s1=sin($2); c1=cos($2); 
print s1,c1}' | awk '{
s2 += $1; c2 += $2} 
END {a1=atan2(s2,c2); 
a2=a1*57.2958; 
if (a2 < 0) a3=-a2;
if (a2 > 0) a3=a2;
if (s2>0 && c2>0) ca=a3;
if (s2>0 && c2<0) ca=a3;
if (s2<0 && c2>0) ca=360-a3;
if (s2<0 && c2<0) ca=360-a3;
print ca}' )

if [ "$GIS_FLAG_G" -eq 1 ]; then
    echo "compass_angle=$(echo $COMPASS_ANGLE)"
else 
    echo "Compass Angle: $(echo $COMPASS_ANGLE)"
fi


LDM=$(v.to.db -p $GIS_OPT_MAP layer=$GIS_OPT_LAYER opt=azimuth units=radians --q | awk -F"|" -v count=${VECT_COUNT} '{
s1=sin($2); c1=cos($2); 
print s1,c1}' | awk '{
s2 += $1; c2 += $2} 
END {a1=atan2(s2,c2); 
a2=a1*57.2958; 
if (a2 < 0) a3=-a2;
if (a2 > 0) a3=a2;
if (s2>0 && c2>0) ldm=90-a3;
if (s2>0 && c2<0) ldm=450-a3;
if (s2<0 && c2>0) ldm=90+a3;
if (s2<0 && c2<0) ldm=90+a3;
print ldm}' )

if [ "$GIS_FLAG_G" -eq 1 ]; then
    echo "directional_mean=$(echo $LDM)"
else 
    echo "Directional Mean: $(echo $LDM)"
fi


VECT_COUNT=$(v.to.db -pc $GIS_OPT_MAP layer=$GIS_OPT_LAYER opt=count --q | cut -d":" -f2 | tr -d [[:blank:]])

CIRC_VAR=$(v.to.db -p $GIS_OPT_MAP layer=$GIS_OPT_LAYER opt=azimuth units=radians --q | awk -F"|" '{
s1=sin($2); c1=cos($2); 
print s1,c1}' | awk -v count="$VECT_COUNT" '{
s2 += $1; c2 += $2} 
END {s3=s2^2; c3=c2^2;
cv=1-((sqrt(s3+c3))/count); 
print cv}')

if [ $VECT_COUNT -gt "1" ]; then
    if [ "$GIS_FLAG_G" -eq 1 ]; then
	echo "circular_variance=$(echo $CIRC_VAR)"
    else    
	echo "Circular Variance: $(echo $CIRC_VAR)"
    fi
fi


### find start/end points of "mean" line 
END_AZIM_DMS=$(echo "${END_AZIM}" | cs2cs +proj=latlong +to +proj=latlong | sed -e 's/d/:/g' -e "s/'/:/g"  -e 's/"//g' | tr -d [[:alpha:]] | tr -d [[:blank:]] | sed 's/:0.000//g')

END_AZIM_NEG=$(echo $END_AZIM | cut -c1)

if [ $END_AZIM_NEG = "-" ]; then
    END_AZIM_DMS=$(echo "-"${END_AZIM_DMS})
fi

START_AZIM=$(echo $END_AZIM | awk '{print 180-$1}')

START_AZIM_DMS=$(echo "${START_AZIM}" | cs2cs +proj=latlong +to +proj=latlong | sed -e 's/d/:/g' -e "s/'/:/g"  -e 's/"//g' | tr -d [[:alpha:]] | tr -d [[:blank:]] | sed 's/:0.000//g')

MEAN_LENGTH=$(v.to.db -p $GIS_OPT_MAP layer=$GIS_OPT_LAYER option=length --q | awk -F"|" ' {s1 = NR} {s2 += $2} {x = s2/s1} END {print x}')

if [ "$GIS_FLAG_G" -eq 1 ]; then
    echo "mean_length=$(echo $MEAN_LENGTH)"
else    
    echo "Mean Length: $(echo $MEAN_LENGTH)" 
fi


HALF_MEAN_LENGTH=$(echo $MEAN_LENGTH | awk '{print $1/2}')

echo "N ${END_AZIM_DMS} E ${HALF_MEAN_LENGTH}" > $TMP1.end.cogo
END_COORDS=$(m.cogo in=$TMP1.end.cogo coord=$(echo ${CENTER_COORDS} | tr ' ' ','))

echo "N ${START_AZIM_DMS} W ${HALF_MEAN_LENGTH}" > $TMP1.start.cogo
START_COORDS=$(m.cogo in=$TMP1.start.cogo coord=$(echo ${CENTER_COORDS} | tr ' ' ','))


if [ "$GIS_FLAG_N" -eq 0 ]; then
    
## make "arrowhead" symbol (if not exist)
    echo "VERSION 1.0
BOX -0.5 -0.5 0.5 0.5
POLYGON
  RING
  FCOLOR NONE
    LINE
      0 0
      0.3 -1
    END
  END
POLYGON
  RING
  FCOLOR NONE
    LINE
      0 0
      -0.3 -1
    END
  END
END" > $TMP1.arrowhead_1

eval $(g.gisenv)

symbol_path="$GISDBASE"/"$LOCATION_NAME"/"$MAPSET"/symbol/arrows
if [ ! -d "$symbol_path" ]; then
    mkdir -p "$symbol_path"
    cp -f $TMP1.arrowhead_1 "$symbol_path"/arrowhead_1
fi


## write LDM graph file and display line of LDM with the arrow
echo "move $(echo ${START_COORDS})" > $TMP1.ldm.graph
echo "width $(echo ${GIS_OPT_WIDTH})" >> $TMP1.ldm.graph 

ARROW_SIZE=$(echo ${GIS_OPT_WIDTH} | awk '{print $1*7}')
ARROW_AZIM=$(echo ${END_AZIM} | awk '{print 360-$1}')

echo "color $(echo ${GIS_OPT_COLOR})" >> $TMP1.ldm.graph
echo "draw $(echo ${END_COORDS})" >> $TMP1.ldm.graph
echo "" >> $TMP1.ldm.graph
echo "rotation $(echo ${ARROW_AZIM})" >> $TMP1.ldm.graph
echo "width $(echo ${GIS_OPT_WIDTH})" >> $TMP1.ldm.graph
echo "symbol "arrows/arrowhead_1" $(echo ${ARROW_SIZE}) $(echo ${END_COORDS}) $(echo ${GIS_OPT_COLOR})" >> $TMP1.ldm.graph

if [ "$GIS_FLAG_S" -eq 0 ] && [ -n "$GIS_OPT_GRAPH" ]; then
    g.message -e "Please use flag <"-s"> to save LDM graphics in file <"$GIS_OPT_GRAPH">" 
    cleanup
    exit 1
fi
 

if [ "$GIS_FLAG_S" -eq 1 ]; then
    if [ -n "$GIS_OPT_GRAPH" ]; then
	cp -f $TMP1.ldm.graph "${GIS_OPT_GRAPH}"
	d.graph -m input="${GIS_OPT_GRAPH}" > /dev/null 2>&1
    else
	eval $(g.gisenv)
	GIS_OPT_GRAPH=""${GISDBASE}"/"${LOCATION_NAME}"/"${GIS_OPT_MAP}"_ldm.graph"
	cat $TMP1.ldm.graph > ${GIS_OPT_GRAPH}
	d.graph -m input="${GIS_OPT_GRAPH}" color=none > /dev/null 2>&1
    fi
else
    cat $TMP1.ldm.graph | d.graph -m > /dev/null 2>&1
fi

fi


## save LDM line to vector if flag "-l" set
if [ "$GIS_FLAG_L" -eq 1 ]; then
    if [ -n "$GIS_OPT_OUTPUT" ]; then
	echo ${START_COORDS} > $TMP1.line.coords
	echo ${END_COORDS} >> $TMP1.line.coords
	v.in.lines input="$TMP1.line.coords" output=${GIS_OPT_OUTPUT} fs=" " --q --o
	v.category in=${GIS_OPT_OUTPUT} out=${GIS_OPT_OUTPUT}_cats --q
	g.rename vect=${GIS_OPT_OUTPUT}_cats,${GIS_OPT_OUTPUT} --o --q > /dev/null 2>&1
	
	if [ -z "$CIRC_VAR" ]; then
	    v.db.addtable map=${GIS_OPT_OUTPUT} col="CompassA double,DirMean double,AveX double, \
AveY double,AveLen double" --q
	else
	    v.db.addtable map=${GIS_OPT_OUTPUT} col="CompassA double,DirMean double,CirVar double, \
AveX double,AveY double,AveLen double" --q
	fi

	echo "UPDATE ${GIS_OPT_OUTPUT} SET CompassA = ${COMPASS_ANGLE};" > $TMP1.ldm.sql 
	echo "UPDATE ${GIS_OPT_OUTPUT} SET DirMean = ${LDM};" >> $TMP1.ldm.sql 
	if [ -n "$CIRC_VAR" ]; then
	    echo "UPDATE ${GIS_OPT_OUTPUT} SET CirVar = ${CIRC_VAR};" >> $TMP1.ldm.sql 
	fi
	echo "UPDATE ${GIS_OPT_OUTPUT} SET AveX = ${CENTER_X};" >> $TMP1.ldm.sql 
	echo "UPDATE ${GIS_OPT_OUTPUT} SET AveY = ${CENTER_Y};" >> $TMP1.ldm.sql 
	echo "UPDATE ${GIS_OPT_OUTPUT} SET AveLen = ${MEAN_LENGTH};" >> $TMP1.ldm.sql 

	cat $TMP1.ldm.sql | db.execute
	
    else
	g.message -e "Please specify <output> file to save LDM vector line"
	cleanup
	exit 1
    fi
    
elif [ "$GIS_FLAG_L" -eq 0 ] && [ -n "$GIS_OPT_OUTPUT" ]; then
    g.message -e "Please use <"-l"> flag to save LDM vector line to <output> file"
    cleanup
    exit 1
fi

cleanup

exit 0