#!/bin/sh
############################################################################
#
# MODULE:       i.landsat.trim
# AUTHOR(S):    Alexander Muriy
#               (Institute of Environmental Geoscience, Moscow, Russia)  
#               e-mail: amuriy AT gmail DOT com 
# VERSION: 0.2
#
#               Some code from page http://grass.osgeo.org/wiki/LANDSAT
#
# PURPOSE:      Trims the "fringe" from the borders of Landsat images, for each band
#               separately or with the MASK where coverage exists for all bands. 
#               Optionally saves vector footprints of trimmed rasters and MASK. 
#               Works with Landsat 5, Landsat 7 (SLC-on). 
#
# COPYRIGHT:    (C) 2012 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: Trims the "fringe" from the borders of Landsat images, for each band separately or with the MASK where coverage exists for all bands. Optionally saves vector footprints of trimmed rasters and MASK. Works with Landsat 5, Landsat 7 (SLC-on).
#% keywords: imagery, landsat, raster, vector
#%End
#%Option
#% key: input
#% type: string
#% required: no
#% multiple: no
#% label: Name of input raster band(s)
#% description: Example: L5170028_02820070521_B10
#% gisprompt: old,cell,raster
#%End
#%Option
#% key: input_base
#% type: string
#% required: no
#% multiple: no
#% label: Base name of input raster bands
#% description: Example: L5170028_02820070521
#%End
#%Option
#% key: input_prefix
#% type: string
#% required: no
#% multiple: no
#% label: Prefix name of input raster bands
#% description: Example: 'B.' for B.1, B.2, ...
#%End
#%Option
#% key: output_prefix
#% type: string
#% required: yes
#% multiple: no
#% label: Prefix for output raster maps 
#% description: Example: 'trim' generates B.1.trim, B.2.trim, ...
#%End
#%Option
#% key: rast_buffer
#% type: integer
#% required: no
#% multiple: no
#% description: Distance for raster buffering (in meters)
#% answer: 300
#%End
#%Option
#% key: gener_thresh
#% type: integer
#% required: no
#% multiple: no
#% description: Threshold for generalizing of vector footprints or coverage MASK (in meters)
#% answer: 3000
#%End
#%Flag
#% key: m
#% description: Trim raster(s) with the MASK where coverage exists for all bands   
#%End
#%Flag
#% key: g
#% description: Trim raster(s) with the generalized footprint from all bands   
#%End
#%Flag
#% key: a
#% description: Process all bands 
#%End
#%Flag
#% key: f
#% description: Save vector footprint(s) of trimmed raster bands or coverage MASK 
#%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

############################################################
cleanup()
{ 
    g.region region=${TMP}_OLD_REGION --q
    g.remove region=${TMP}_OLD_REGION --q
    g.mremove -f rast=MASK,"*I_LANDSAT_TRIM*" --q
    g.mremove -f vect="*I_LANDSAT_TRIM*" --q
}
############################################################
## what to do in case of user break:
exitprocedure()
{
    echo ""
    echo "User break!"
    cleanup
    exit 1
}

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

############################################################
## functions ##

TMP=I_LANDSAT_TRIM

band_mask() 
{
    rast_thresh=$(($GIS_OPT_RAST_BUFFER/3))
    
    if [ $GIS_FLAG_M -eq 1 ]; then
	g.region rast=${TMP}.cover res=$rast_thresh
	g.message message="Creating raster mask ..."
    	r.mapcalc "${TMP}.mask = if(${TMP}.cover > 0, 1, null())"
    else
	g.region rast=$MAP res=$rast_thresh
	if [ $GIS_FLAG_G -eq 0 ]; then	
	    g.message message="Creating raster mask ..."
	fi
	r.mapcalc "${TMP}.mask = if($MAP > 0, 1, null())"
    fi

    if [ $GIS_FLAG_G -eq 0 ]; then
	g.message message="Creating raster buffer with distance <"$GIS_OPT_RAST_BUFFER"> ..."
    fi
    r.buffer in=${TMP}.mask out=${TMP}.buffer distances="$GIS_OPT_RAST_BUFFER" --o --q
    r.mapcalc "${TMP}.buffer_2 = if(${TMP}.buffer == 2, 1, null())"
    g.region res="$GIS_OPT_RAST_BUFFER"
    if [ $GIS_FLAG_G -eq 0 ]; then
	g.message message="Growing raster and cut mask ..."
    fi
    r.grow in=${TMP}.buffer_2 out=${TMP}.grow radius=1 --o --q
    r.mapcalc "${TMP}.mask.cut = if(isnull(${TMP}.grow), ${TMP}.mask, null())" 
    if [ $GIS_FLAG_G -eq 0 ]; then
	g.message message="Vectorize raster mask ..."
    fi
    r.to.vect in=${TMP}.mask.cut out=${TMP}_mask feature=area --o --q > /dev/null 2>&1
    g.rename vect=${TMP}_mask,${TMP}_mask_nocentr --o --q > /dev/null 2>&1
    v.edit ${TMP}_mask_nocentr tool=delete type=centroid cats=0-999999 --o --q
    v.centroids in=${TMP}_mask_nocentr out=${TMP}_mask --o --q > /dev/null 2>&1
}

band_mask_gener() 
{
    main_poly=$(v.to.db -p ${TMP}_mask opt=area --q \
    	| sort -n -t'|' -k2 | tail -n1 | cut -d'|' -f1)
    if [ $GIS_FLAG_G -eq 0 ]; then 
	g.message message="Extract main area from mask ..."
    fi
    v.extract in=${TMP}_mask out=${TMP}_main list=$main_poly --o --q > /dev/null 2>&1
    
    v.type in=${TMP}_main out=${TMP}_lines type=boundary,line --o --q > /dev/null 2>&1
    v.category in=${TMP}_lines out=${TMP}_lines_nocats type=line --o --q > /dev/null 2>&1
    v.category in=${TMP}_lines_nocats out=${TMP}_lcats type=line --o --q > /dev/null 2>&1
    v.edit ${TMP}_lcats tool=merge type=line cats=0-99999 --q
    v.edit ${TMP}_lcats tool=delete type=centroid cats=0-99999 --q
    
    if [ $GIS_FLAG_G -eq 0 ]; then 
	g.message message="Generalizing vectors ..."
    fi
    thresh=$(($GIS_OPT_GENER_THRESH/20))
    v.generalize in=${TMP}_lcats out=${TMP}_lcats_gener method=douglas thresh=$thresh --o --q > /dev/null 2>&1
    
    thresh2=$(($GIS_OPT_GENER_THRESH/5))
    v.generalize in=${TMP}_lcats_gener out=${TMP}_lcats_gener2 method=douglas \
    	threshold=$thresh2 --o --q > /dev/null 2>&1
    
    thresh3=$(($GIS_OPT_GENER_THRESH/3))
    v.generalize in=${TMP}_lcats_gener2 out=${TMP}_lcats_gener3 method=douglas \
    	threshold=$thresh3 --o --q > /dev/null 2>&1
    
    v.type in=${TMP}_lcats_gener3 out=${TMP}_bounds type=line,boundary --o --q > /dev/null 2>&1
    v.centroids in=${TMP}_bounds out=${TMP}_mask_gener --o --q > /dev/null 2>&1
    
    if [ $GIS_FLAG_F -eq 1 ]; then
	if [ $GIS_FLAG_G -eq 0 ]; then
	    VECT_NAME=$(echo "${MAP}_${GIS_OPT_OUTPUT_PREFIX}" | tr '[[:punct:]]' '_') 
	    g.message message="Save vector mask to <"$VECT_NAME"> ..."
	    g.rename vect=${TMP}_mask_gener,${VECT_NAME} --o --q > /dev/null 2>&1
	    g.message message="------------------------------------------------------------"
	fi
    fi
    
    if [ $GIS_FLAG_G -eq 1 ]; then
	VECT_NAME="${MAP}_${TMP}_gener"
	g.rename vect=${TMP}_mask_gener,${VECT_NAME} --o --q > /dev/null 2>&1
    fi
}

band_mask_cut() 
{
    g.region rast="$MAP"
    g.message message="Rasterizing vector mask ..."
    v.to.rast in=${TMP}_mask_gener out=${TMP}_mask_gener use=cat --o --q > /dev/null 2>&1
    g.message message="Write trimmed raster to <"${MAP}.${GIS_OPT_OUTPUT_PREFIX}"> ..."
    r.mask in=${TMP}_mask_gener --q 
    r.mapcalc ""${MAP}.${GIS_OPT_OUTPUT_PREFIX}" = $MAP"
    r.mask -r --q
    g.message message="------------------------------------------------------------"
}

coverage_mask() 
{
    g.region rast=$(g.mlist rast pat="${GIS_OPT_INPUT_BASE}*${GIS_OPT_INPUT_PREFIX}*" sep=',') 
    g.message message="Creating coverage raster mask ..."        
    r.series in=$(g.mlist rast pat="${GIS_OPT_INPUT_BASE}*${GIS_OPT_INPUT_PREFIX}*" \
	sep=',') -n out=${TMP}.thresh method=threshold thresh=1 --o 
    r.mapcalc "${TMP}.cover = if(isnull(${TMP}.thresh))"
}

coverage_mask_gener()
{
    band_mask
    
    v.centroids in=${TMP}_mask out=${TMP}_fill --o --q > /dev/null 2>&1
    v.edit ${TMP}_fill tool=delete type=centroid cats=0-999999 --o --q
    v.category in=${TMP}_fill out=${TMP}_acats option=add cat=1 step=0 type=area  --o --q > /dev/null 2>&1
    v.extract in=${TMP}_acats out=${TMP}_extract type=area new=1 list=1-999999  --o --q > /dev/null 2>&1
    v.dissolve in=${TMP}_extract out=${TMP}_dissolve --o --q > /dev/null 2>&1
    v.category in=${TMP}_dissolve out=${TMP}_dissolve_nocats opt=del type=centroid --o --q > /dev/null 2>&1
    v.category in=${TMP}_dissolve_nocats out=${TMP}_dissolve_newcats type=centroid --o --q > /dev/null 2>&1
    v.db.addtable ${TMP}_dissolve_newcats col="area double" --o --q > /dev/null 2>&1
    v.to.db ${TMP}_dissolve_newcats opt=area col=area --o --q > /dev/null 2>&1
    
    main_poly=$(v.to.db -p ${TMP}_dissolve_newcats opt=area --q \
	| sort -n -t'|' -k2 | tail -n1 | cut -d'|' -f1)
    g.message message="Extract main area from mask ..."
    v.extract in=${TMP}_dissolve_newcats out=${TMP}_cut list=$main_poly --o --q > /dev/null 2>&1

    g.message message="Generalizing vectors ..."
    v.clean in=${TMP}_cut out=${TMP}_gener tool=prune thresh="$GIS_OPT_GENER_THRESH" --o --q > /dev/null 2>&1
    v.type in=${TMP}_gener out=${TMP}_lines type=boundary,line --o --q > /dev/null 2>&1
    v.category in=${TMP}_lines out=${TMP}_lcats --o --q > /dev/null 2>&1
    v.edit ${TMP}_lcats tool=merge type=line cats=0-999999  --q
    v.edit ${TMP}_lcats tool=delete type=centroid cats=0-999999 --q 

    thresh2=$(($GIS_OPT_GENER_THRESH/3))
    v.generalize in=${TMP}_lcats out=${TMP}_lcats_gener method=douglas threshold=$thresh2 --o --q > /dev/null 2>&1
    v.type in=${TMP}_lcats_gener out=${TMP}_bounds type=line,boundary --o --q > /dev/null 2>&1
    v.centroids in=${TMP}_bounds out=${TMP}_area --o --q > /dev/null 2>&1
    
    if [ $GIS_FLAG_F -eq 1 ]; then
	g.message message="Save coverage mask to <"${GIS_OPT_INPUT_BASE}_cover_mask"> ..."
	g.rename vect=${TMP}_area,${GIS_OPT_INPUT_BASE}_cover_mask --o --q
	g.message message="------------------------------------------------------------"
    fi
}

coverage_mask_cut()
{
    g.region rast=$MAP
    
    eval $(g.gisenv)
    eval $(g.findfile element=cell mapset=$MAPSET file="${TMP}_area")
    if [ ! "$file" ]; then
	g.message message="Rasterizing coverage mask ..."
	v.to.rast in=${TMP}_area out=${TMP}_area use=cat --q > /dev/null 2>&1
    fi
    
    r.mask in=${TMP}_area --q
    g.message message="Write trimmed raster to <"${MAP}.${GIS_OPT_OUTPUT_PREFIX}"> ..."
    r.mapcalc "${MAP}.${GIS_OPT_OUTPUT_PREFIX} = ${MAP}"
    r.mask -r --q
    g.message message="------------------------------------------------------------"
}

gener_footprint()
{
    g.message message="=== Create generalized footprint for all bands  ==="
    g.region res="$GIS_OPT_RAST_BUFFER"

    g.mlist vect pat="${GIS_OPT_INPUT_BASE}*_${TMP}_gener" | while read TRIM; do 
	start=$count
	v.type in="$TRIM" out="${TMP}_line${count}" type=boundary,line --o --q 
	v.edit "${TMP}_line${count}" tool=delete type=centroid cats=0-999 --q 
	v.to.rast in="${TMP}_line${count}" out="${TMP}_rast${count}" use=cat --o --q
	r.grow in="${TMP}_rast${count}" out="${TMP}_grow${count}" radius=3 --o --q
	count=$(($start+1))
    done
    
    for patch in $(g.mlist rast pat="${TMP}_grow*" sep=','); do
    	r.patch in="$patch" out="${TMP}_patch" --o --q 
    done
    
    r.mapcalc "${TMP}_patch_null = if(isnull(${TMP}_patch),1,null())" 
    r.to.vect in="${TMP}_patch_null" out="${TMP}_patch_null" feature=area --o --q
    
    main_poly=$(v.to.db -p "${TMP}_patch_null"  opt=area --q \
    	| sort -n -t'|' -k2 | tail -n1 | cut -d'|' -f1)
    
    v.extract in="${TMP}_patch_null" out="${TMP}_patch_null_main" list="$main_poly" --o --q
    thresh=$(($GIS_OPT_GENER_THRESH/5))
    v.generalize in="${TMP}_patch_null_main" out="${TMP}_patch_null_gener" method=douglas threshold=$thresh --o --q > /dev/null 2>&1
    dist_buffer=$(($GIS_OPT_RAST_BUFFER*3))
    v.buffer -s in="${TMP}_patch_null_gener" out="${TMP}_patch_buffer" distance="$dist_buffer" --o --q 
    
    if [ $GIS_FLAG_F -eq 1 ]; then
    	g.message message="Save generalized footprint to <"${GIS_OPT_INPUT_BASE}_gener_mask"> ..."
    	g.rename vect=${TMP}_patch_buffer,${GIS_OPT_INPUT_BASE}_gener_mask --o --q 
    	g.message message="------------------------------------------------------------"
    fi
}

gener_footprint_cut()
{
    g.region rast="$MAP"
    g.message message="Rasterizing generalized vector footprint ..." 
    v.to.rast in=${TMP}_patch_buffer out=${TMP}_patch_buffer use=cat --o --q > /dev/null 2>&1 
    g.message message="Write trimmed raster to <"${MAP}.${GIS_OPT_OUTPUT_PREFIX}"> ..."
    r.mask in=${TMP}_patch_buffer --q  
    r.mapcalc ""${MAP}.${GIS_OPT_OUTPUT_PREFIX}" = $MAP"
    r.mask -r --q
    g.message message="------------------------------------------------------------"
}


############################################################
## main ##

g.region save=${TMP}_OLD_REGION --q

oIFS="$IFS"

if [ $GIS_FLAG_M -eq 1 -o $GIS_FLAG_G -eq 1 ] && [ -z "$GIS_OPT_INPUT_BASE" ]; then
    g.message -e message="Use flags \"-g\" and \"-m\" with option \"input_base\" (and optionally with \"input_prefix\")"
    cleanup
    exit 1
fi

if [ $GIS_FLAG_M -eq 1 -a $GIS_FLAG_G -eq 1 ]; then
    g.message -e message="Don't use flags \"-g\" and \"-m\" together"
    cleanup
    exit 1
fi

if [ $GIS_FLAG_M -eq 1 -a $GIS_FLAG_F -eq 0 ] && [ -z "$GIS_OPT_INPUT" -a $GIS_FLAG_A -eq 0 ]; then
    g.message -e message="Specify raster(s) with option \"input\" or use flag \"-a\" to choose all rasters"
    cleanup
    exit 1
fi

if [ $GIS_FLAG_G -eq 1 ] && [ -z "$GIS_OPT_INPUT" -a $GIS_FLAG_A -eq 0 ]; then
    g.message -e message="Specify raster(s) with option \"input\" or use flag \"-a\" to choose all rasters"
    cleanup
    exit 1
fi

if [ $GIS_FLAG_A -eq 1 ]; then
    if [ $GIS_FLAG_G -eq 1 -a $GIS_FLAG_F -eq 1 ]; then
	g.message -e message="Don't use flag \"-a\" with flags \"-f\" and  \"-g\" "
	cleanup
	exit 1
    fi 
    if [ $GIS_FLAG_M -eq 1 -a $GIS_FLAG_F -eq 1 ]; then
	g.message -e message="Don't use flag \"-a\" with flags \"-f\" and  \"-m\" "
	cleanup
	exit 1
    fi 
fi  

if [ $GIS_FLAG_A -eq 1 ]; then
    if [ -z "$GIS_OPT_INPUT" -a -n "$GIS_OPT_INPUT_BASE" ]; then
	GIS_OPT_INPUT=$(g.mlist rast pat="${GIS_OPT_INPUT_BASE}*${GIS_OPT_INPUT_PREFIX}*" sep=',')
    elif [ -z "$GIS_OPT_INPUT" -a -z "$GIS_OPT_INPUT_BASE" ]; then
	g.message -e message="Choose input raster(s) with \"input\" option or with \"input_base\" and \"input_prefix\" options"
	cleanup
	exit 1
    fi
fi

if [ $GIS_FLAG_M -eq 1 ]; then
    coverage_mask
    coverage_mask_gener
    
    IFS=','

    if [ $GIS_FLAG_F -eq 0 ]; then
	for MAP in $GIS_OPT_INPUT ; do
	    g.message message="=== Processing raster: <"$MAP"> ===" 
	    coverage_mask_cut
	done
    fi

    IFS="$oIFS"

elif [ $GIS_FLAG_G -eq 0 ]; then    
    
    IFS=','
    
    for MAP in $GIS_OPT_INPUT ; do
	g.message message="=== Processing raster: <"$MAP"> ==="
	band_mask
	band_mask_gener	
	
	if [ $GIS_FLAG_F -eq 0 ]; then
	    band_mask_cut
	fi
    done
    
    IFS="$oIFS"
fi 

if [ $GIS_FLAG_G -eq 1 ]; then
    if [ -n "$GIS_OPT_INPUT_BASE" ]; then
	GENER_INPUT=$(g.mlist rast pat="${GIS_OPT_INPUT_BASE}*${GIS_OPT_INPUT_PREFIX}*" sep=',')
    fi
    
    IFS=','
    
    g.message message="=== Make rasters masks and generalizing them ==="
    for MAP in $GENER_INPUT; do
	g.message message="Processing raster <"$MAP">"
	band_mask
	band_mask_gener
    done
    
    IFS="$oIFS"

    gener_footprint
    
    if [ "$GIS_FLAG_F" -eq 0 ]; then
	if [ -z "$GIS_OPT_INPUT" ]; then
	    GENER_F_INPUT=$(g.mlist rast pat="${GIS_OPT_INPUT_BASE}*${GIS_OPT_INPUT_PREFIX}*" sep=',')
	else
	    GENER_F_INPUT="$GIS_OPT_INPUT"
	fi
	
	IFS=','

	for MAP in $GENER_F_INPUT; do
	    gener_footprint_cut
 	done
	IFS="$oIFS"
    fi
fi	

IFS="$oIFS"
	
# cleanup
cleanup

exit 0