#!/bin/sh

# s.to.poly.sh - create coarse-grained area vectors around site data
# 
#   tom poindexter, feb 2000
#
#   usage: s.to.poly.sh input=<sitefile> output=<vector> 
#     options:
#           filter=<filtername> | filter=none (default = 11x11 isolation) 
#           maxpoly=<maximum number of resulting area polygons> (default=99999)
#           maxiter=<maximum iterations> (default=5)

#### Check that user is running GRASS
if test "$GISRC" = ""
then
        echo "$0: You must be running GRASS to execute $0" >&2
        exit 1
fi

#### Check GRASS variables are set
eval `g.gisenv`

# our variables
site=''
poly=''
filt=''
maxpoly=99998
maxiter=5

# process args
args="$*"
for arg in $args
do
    SAVEIFS="$IFS"
    IFS='='
    set $arg
    IFS="$SAVEIFS"
    case "$1" in
        input)
	    site=$2
	    ;;
        output)
	    poly=$2
	    ;;
        maxpoly)
	    maxpoly=$2
	    ;;
        maxiter)
	    maxiter=$2
	    ;;
        filter)
	    filt=$2
            if [ "$filt" != "none" ] ; then
                if [ ! -f $filt ] ; then
                    echo "filter $filt not found"
                    exit
                fi
            fi
	    ;;
    esac
done

if [ -z "$site" -o -z "$poly" ] ; then
    echo "usage:  $0  input=<sitefile>  output=<vector>"
    echo "options:"
    echo "     filter=<filtername> | filter=none (default = 11x11 isolation)" 
    echo "     maxpoly=<maximum number of resulting polygons> (default=99999)"
    echo "     maxiter=<maximum iterations> (default=1)"
    exit
fi

# temp name for intermediate rasters
tmp=tmp__$$

# convert site file to raster
echo "converting site file to raster"
s.to.rast -s -q input=$site output=$tmp.r
input=$tmp.r

if [ -z "$filt" ] ; then 
    filt=/tmp/$tmp.filt
    # create a smoothing filter to throw out isolated points
    cat >/tmp/$tmp.filt <<EOF
MATRIX 11
1 1 1 1 1 1 1 1 1 1 1
1 1 1 2 2 2 2 2 1 1 1
1 1 2 2 3 3 3 2 2 1 1
1 2 2 3 3 3 3 3 2 2 1
1 2 3 3 0 0 0 3 3 2 1
1 2 3 3 0 0 0 3 3 2 1
1 2 3 3 0 0 0 3 3 2 1
1 2 2 3 3 3 3 3 2 2 1
1 1 2 2 3 3 3 2 2 1 1
1 1 1 2 2 2 2 2 1 1 1
1 1 1 1 1 1 1 1 1 1 1
DIVISOR 50
TYPE P
EOF
fi

# filter raster file to clean up outliers
if [ "$filt" != "none" ] ; then
    echo "filtering raster"
    r.mfilter input=$tmp.r output=$tmp.f filter=$filt >/dev/null 2>&1
    input=$tmp.f
    g.remove rast=$tmp.r >/dev/null
fi
rm -f /tmp/$tmp.filt

# reduce raster to 0/1 values
r.mapcalc "$tmp.i0 = if($input)" >/dev/null 2>&1
g.remove rast=$tmp.f >/dev/null

curr_rast=$tmp.i0
curr_cnt=99999
iter_cnt=0

# iterate conversion to polygons, growing each time to reduce polygon count

while [ $curr_cnt -gt $maxpoly -a $maxiter -gt $iter_cnt ] ; do

    iter_cnt=`expr $iter_cnt + 1`
    echo "iteration $iter_cnt"

    # convert current raster to smooth polygons (-l)
    r.poly input=$curr_rast output=$poly -l >/dev/null

    # see how many polygon where created
    area_line=`v.support option=build map=$poly 2>/dev/null | \
						grep 'Number of areas:'`
    set -- $area_line
    curr_cnt=$4
    
    echo "    polygon count = $curr_cnt"

    if [ $curr_cnt -gt $maxpoly -a $maxiter -gt $iter_cnt ] ; then
        echo "    growing raster"

        # cycle raster through r.grow, binary 0/1 (-b) output
        next=$tmp.i$iter_cnt
        r.grow input=$curr_rast output=$next -b >/dev/null 2>&1
        g.remove rast=$curr_rast >/dev/null
        curr_rast=$next
    fi

done

g.remove rast=$curr_rast >/dev/null

echo "done."