#!/bin/sh

############################################################################
#
# MODULE:	i.fusion.brovey
# AUTHOR(S):	Markus Neteler. <neteler itc it>
# PURPOSE:	Brovey transform to merge
#                 - LANDSAT-7 MS (2, 4, 5) and pan (high res)
#                 - SPOT MS and pan (high res)
#                 - QuickBird MS and pan (high res)
#
# COPYRIGHT:	(C) 2002,2004 by the GRASS Development Team
#
#		This program is free software under the GNU General Public
#		License (>=v2). Read the file COPYING that comes with GRASS
#		for details.
#
# REFERENCES:
#             (?) Roller, N.E.G. and Cox, S., 1980. Comparison of Landsat MSS
#                and merged MSS/RBV data for analysis of natural vegetation.
#                Proc. of the 14th International Symposium on Remote Sensing
#                of Environment, San Jose, Costa Rica, 23-30 April, pp. 1001-1007.
#
#               for LANDSAT 5: see Pohl, C 1996 and others
#
# TODO:         add overwrite test at beginning of the script
#############################################################################

#%Module
#%  description: Brovey transform to merge multispectral and high-res panchromatic channels
#%  keywords: raster, imagery, fusion
#%End
#%Flag
#%  key: l
#%  description: sensor: LANDSAT
#%END
#%Flag
#%  key: q
#%  description: sensor: QuickBird
#%END
#%Flag
#%  key: s
#%  description: sensor: SPOT
#%END
#%option
#% key: ms1
#% type: string
#% gisprompt: old,cell,raster
#% description: raster input map (green: tm2 | qbird_green | spot1)
#% required : yes
#%end
#%option
#% key: ms2
#% type: string
#% gisprompt: old,cell,raster
#% description: raster input map (NIR: tm4 | qbird_nir | spot2)
#% required : yes
#%end
#%option
#% key: ms3
#% type: string
#% gisprompt: old,cell,raster
#% description: raster input map (MIR; tm5 | qbird_red | spot3)
#% required : yes
#%end
#%option
#% key: pan
#% type: string
#% gisprompt: old,cell,raster
#% description: raster input map (etmpan | qbird_pan | spotpan)
#% required : yes
#%end
#%option
#% key: outputprefix
#% type: string
#% gisprompt: new,cell,raster
#% description: raster output map prefix (e.g. 'brov')
#% required : yes
#%end

# 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

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

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

PROG=`basename $0`

#### check if we have awk
if [ ! -x "`which awk`" ] ; then
    echo "$PROG: awk required, please install awk or gawk first" 2>&1
    exit 1
fi

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

cleanup()
{
 #restore settings:
 g.region region=i.fusion.brovey.$TMP
 g.remove region=i.fusion.brovey.$TMP > /dev/null
}

TMP=$$

if [ $GIS_FLAG_S -eq 0 -a $GIS_FLAG_L -eq 0 -a $GIS_FLAG_Q -eq 0 ] ; then
   echo "Please select a flag to specify the satellite sensor"
   exit 1
fi

#get PAN resolution:
PANRES=`r.info $GIS_OPT_PAN | grep Res | tr -s ' ' ' ' | cut -d':' -f4 | cut -d' ' -f2 | tr '\n' ' '  | awk '{printf "%f\n", ($1+$2) / 2.}'`
if [ $? -ne 0 ] ; then
 echo "An error occurred."
 exit 1
fi

#save current settings:
g.region save=i.fusion.brovey.$TMP > /dev/null

echo "Temporarily setting raster resolution to PAN resolution: $PANRES"
g.region res=$PANRES -ap > /dev/null

echo "Performing Brovey transformation..."

# The formula was originally developed for LANDSAT-TM5 and SPOT, 
# but it also works well with LANDSAT-TM7
# LANDSAT formula:
#  r.mapcalc "brov.red=1. *  tm.5 / (tm.2 + tm.4 + tm.5) * etmpan"
#  r.mapcalc "brov.green=1. * tm.4 /(tm.2 + tm.4 + tm.5) * etmpan"
#  r.mapcalc "brov.blue=1. * tm.2 / (tm.2 + tm.4 + tm.5) * etmpan"
#
# SPOT formula:
# r.mapcalc "brov.red= 1.  * spot.ms.3 / (spot.ms.1 + spot.ms.2 + spot.ms.3) * spot.p"
# r.mapcalc "brov.green=1. * spot.ms.2 / (spot.ms.1 + spot.ms.2 + spot.ms.3) * spot.p"
# r.mapcalc "brov.blue= 1. * spot.ms.1 / (spot.ms.1 + spot.ms.2 + spot.ms.3) * spot.p"
# note: for RGB composite then revert brov.red and brov.green!

echo -n "  Calculating $GIS_OPT_OUTPUTPREFIX.red: ...(step 1/4)"
r.mapcalc   $GIS_OPT_OUTPUTPREFIX.red="1. * $GIS_OPT_MS3 / ($GIS_OPT_MS1 + $GIS_OPT_MS2 + $GIS_OPT_MS3) * $GIS_OPT_PAN"
if [ $? -ne 0 ] ; then
 echo "An error occurred."
 exit 1
fi

echo -n "  Calculating $GIS_OPT_OUTPUTPREFIX.green: ...(step 2/4)"
r.mapcalc   $GIS_OPT_OUTPUTPREFIX.green="1. * $GIS_OPT_MS2 / ($GIS_OPT_MS1 + $GIS_OPT_MS2 + $GIS_OPT_MS3) * $GIS_OPT_PAN"
if [ $? -ne 0 ] ; then
 echo "An error occurred."
 exit 1
fi

echo -n "  Calculating $GIS_OPT_OUTPUTPREFIX.blue: ...(step 3/4)"
r.mapcalc   $GIS_OPT_OUTPUTPREFIX.blue="1. * $GIS_OPT_MS1 / ($GIS_OPT_MS1 + $GIS_OPT_MS2 + $GIS_OPT_MS3) * $GIS_OPT_PAN"
if [ $? -ne 0 ] ; then
 echo "An error occurred."
 exit 1
fi

# Maybe?
#r.colors   $GIS_OPT_OUTPUTPREFIX.red col=grey
#r.colors   $GIS_OPT_OUTPUTPREFIX.green col=grey
#r.colors   $GIS_OPT_OUTPUTPREFIX.blue col=grey
#to blue-ish, therefore we modify
#r.colors $GIS_OPT_OUTPUTPREFIX.blue col=rules << EOF
#5 0 0 0
#20 200 200 200
#40 230 230 230
#67 255 255 255
#EOF

if [ $GIS_FLAG_S -eq 1 ] ; then
 #apect table is nice for SPOT:
 echo "Assigning color tables for SPOT ...(step 4/4)"
 r.colors $GIS_OPT_OUTPUTPREFIX.blue col=aspect
 r.colors $GIS_OPT_OUTPUTPREFIX.green col=aspect
 r.colors $GIS_OPT_OUTPUTPREFIX.red col=aspect
 echo "Fixing output names..."
 g.rename $GIS_OPT_OUTPUTPREFIX.green,$GIS_OPT_OUTPUTPREFIX.$TMP > /dev/null
 g.rename $GIS_OPT_OUTPUTPREFIX.red,$GIS_OPT_OUTPUTPREFIX.green > /dev/null
 g.rename $GIS_OPT_OUTPUTPREFIX.$TMP,$GIS_OPT_OUTPUTPREFIX.red > /dev/null
else
 #aspect table is nice for LANDSAT and QuickBird:
 echo "Assigning color tables for LANDSAT or QuickBird ...(step 4/4)"
 r.colors $GIS_OPT_OUTPUTPREFIX.blue col=aspect
 r.colors $GIS_OPT_OUTPUTPREFIX.green col=aspect
 r.colors $GIS_OPT_OUTPUTPREFIX.red col=aspect
fi

echo "Restoring region settings."
cleanup

echo "Following pan-sharpened output maps have been generated:"
echo "  $GIS_OPT_OUTPUTPREFIX.red"
echo "  $GIS_OPT_OUTPUTPREFIX.green"
echo "  $GIS_OPT_OUTPUTPREFIX.blue"
echo ""
echo "To visualize output, run:"
echo "  g.region -p rast=$GIS_OPT_OUTPUTPREFIX.red ; d.erase"
echo "  d.rgb r=$GIS_OPT_OUTPUTPREFIX.red g=$GIS_OPT_OUTPUTPREFIX.green b=$GIS_OPT_OUTPUTPREFIX.blue"
echo ""
echo "If desired, combine channels with 'r.composite' to a single map"