#!/bin/sh
############################################################################
#
# MODULE:       i.landsat.trim
# AUTHOR(S):    Alexander Muriy
#               (Institute of Environmental Geoscience, Moscow, Russia)  
#               e-mail: amuriy AT gmail DOT com 
#
#               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: 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 "User break!"
    cleanup
    exit 1
}

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

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

TMP=I_LANDSAT_TRIM

if_exist()
{
    eval "$(g.gisenv)"
    eval "$(g.findfile mapset="$MAPSET" element=cell file="$MAP")"
    if [ ! "$file" ]; then
	g.message -e "Raster map <"$MAP"> not found"
	cleanup
	exit 1
    fi
}

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
	g.message message="Creating raster mask ..."
	r.mapcalc "${TMP}.mask = if($MAP > 0, 1, null())"
    fi
    
    g.message message="Creating raster buffer with distance <"$GIS_OPT_RAST_BUFFER"> ..."
    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"
    g.message message="Growing raster and cut mask ..."
    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())" 
    g.message message="Vectorize raster mask ..."
    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)
    g.message message="Extract main area from mask ..."
    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
    
    g.message message="Generalizing vectors ..."
    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
	VECT_NAME=$(echo "${MAP}_${GIS_OPT_OUTPUT_PREFIX}" | tr '[[:punct:]]' '_') 
	g.message message="Save vector mask to <"$VECT_NAME"> ..."
	g.copy vect=${TMP}_mask_gener,${VECT_NAME} --o --q > /dev/null 2>&1
	g.message message="------------------------------------------------------------"
    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.copy vect=${TMP}_area,${GIS_OPT_INPUT_BASE}_cover_mask --o --q
	g.message message="------------------------------------------------------------"
    fi
}

coverage_mask_cut()
{
    g.region rast=$MAP
    g.message message="Rasterizing coverage mask ..."
    v.to.rast in=${TMP}_area out=${TMP}_area use=cat --o --q > /dev/null 2>&1
    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="------------------------------------------------------------"
}

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

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

if [ $GIS_FLAG_A -eq 1 ]; then
    if [ -n "$GIS_OPT_INPUT" ]; then
	g.message -e "Use flag \"-a\" with options \"input_base\" and \"input_prefix\", but without \"input\""
	cleanup
	exit 1
    elif [ -n "$GIS_OPT_INPUT_BASE" ]; then
	GIS_OPT_INPUT=$(g.mlist rast pat="${GIS_OPT_INPUT_BASE}*${GIS_OPT_INPUT_PREFIX}*" sep=',')
    fi
    if [ $GIS_FLAG_F -eq 1 -a $GIS_FLAG_M -eq 1 ]; then
	g.message -e "Use flag \"-f\" without flag \"-a\""
	cleanup
	exit 1
    fi
fi

if [ -n "$GIS_OPT_INPUT" ]; then 
    if [ $GIS_FLAG_M -eq 1 ]; then
	if [ -n "$GIS_OPT_INPUT_BASE" ]; then
	    coverage_mask
	    
	    ORIG_IFS="$IFS"
	    IFS=','
	    
	    for MAP in $GIS_OPT_INPUT ; do
	    	if_exist
		g.message message="=== Processing raster: <"$MAP"> ==="
	    	coverage_mask_gener
		
		if [ $GIS_FLAG_F -eq 0 ]; then
		    coverage_mask_cut
		fi
	    done
	    
	    IFS="$ORIG_IFS"

	else
	    g.message -e "Use flag \"-m\" with options \"input_base\" and \"input_prefix\""
	    cleanup
	    exit 1
	fi

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

elif [ $GIS_FLAG_M -eq 1 ] && [ $GIS_FLAG_F -eq 1 ]; then
    if [ -n "$GIS_OPT_INPUT_BASE" ]; then
	coverage_mask
	coverage_mask_gener
    fi
fi    




## cleanup
cleanup

exit 0