#!/bin/bash

set -e
set -o pipefail
pushd `dirname "$0"` > /dev/null
DIR=`pwd`
popd > /dev/null
export PATH="$DIR:$PATH"

if [ "$#" -lt 5 ]; then

    if [ "$#" -eq 1 ]; then
	if [ "$1" == "-shortDescription" ]; then
	    echo "Generate a mapola from a regional lidar data set."
	    exit 0
	fi	
    fi

    echo "This program takes lidar data within a small region of the asteroid (in binary format as expected" 
    echo "by GMT blockmedian) and runs it through GMT's blockmedian and surface to generate a mapola. This"
    echo "program also converts the output of GMT (which is NetCDF format) into FITS format."
    echo
    echo "Usage: GMTSurface <input-file> <pixel-spacing> <bounds> <output-grd> <output-fit> [...]"
    echo
    echo "Where:"
    echo "  <input-file>    = binary file storing lidar data"
    echo "  <pixel-spacing> = pixels spacing in units of lidar data (GSD)"
    echo "  <bounds>        = bounds in units of lidar data to be passed to GMT's -R option: E.g. 0.5/2.5/1.25/3.25"
    echo "  <output-grd>    = path to GMT output (NetCDF format)"
    echo "  <output-fit>    = path to FIT output"
    echo "  [...]           = remaining arguments are passed unmodified to gmt surface"
    echo "                    '-fg' ensures geographic information is added to the NETCDF header for Long/Lat/Rad data"
    echo "                    'smooth' = blockmedian is done over 5 times the spaceing, followed by surface."
    echo "                               Otherwise blockmedian and surface have the same spacing."
    echo "                    'near'   = nearestneighbor algorithm is run over a radius 10 times spaceing"
    echo "                               where atleast one return is found in 50% of the radius (otherwise NaN)."
    echo
    exit 1
fi
echo "Starting GMTSurface"
gmtInputFile=$1
spacing=$2
bounds=$3
grd=$4
xyz=$4.xyz
fits=$5

#For lat-long cases
arr=($@)
for i in ${arr[@]}
do 
    if [[ $i == *"-f"* ]]
    then
	fgmt=$i
    fi
done


DIR=$(dirname $fits)
if [ ${#DIR} -gt 0 ] 
then  
    filename=${grd##*/}
    xgrd=$DIR/num_$filename
    ygrd=$DIR/std_$filename
else
    xgrd=num_$4
    ygrd=std_$4
fi

shift
shift
shift
shift
shift


cmd="blockmedian"
spc=$spacing
echo $cmd
echo $spc
# Mean case
arr=($@)
echo $arr
echo ${#arr[@]}
if [[ ${#arr[@]} -gt 0 ]] 
then 
    echo Hi $@
    for i in ${arr[@]}
    do 
	if [[ $i == *"smooth"* ]]
	then
	    spc=$(echo $spacing*3. | bc)
	    shift
	fi
    done
fi
echo $cmd
echo $spc

flag=1
#nearestneighbor
echo $@
for i in ${arr[@]}
do 
    if [[ $i == *"near"* ]]
    then
	# Run nearest neighbor algorithm
	search_radius=$(echo $spacing*10. | bc)
	echo $search_radius
	echo "gmt nearneighbor -N4  -bi4d $gmtInputFile -I$spacing -R$bounds -V -G$grd -S$search_radius"
	gmt nearneighbor -N4 -bi4d $gmtInputFile -I$spacing -R$bounds -V -G$grd -S$search_radius
	flag=0
    fi
done

if [[ $flag == 1 ]] 
then
    # Run GMT's surface to generate a mapola
    echo "gmt $cmd -bi4d $gmtInputFile -I$spc -R$bounds -V -bo | gmt surface -bi3 -I$spacing -R$bounds -V -G$grd $@"
    gmt $cmd -bi4d $gmtInputFile -I$spc -R$bounds -V -bo | gmt surface -bi3 -I$spacing -R$bounds -V -G$grd $@
fi


gmtVersion=$(gmt --version)
gmtMajor=$(echo $gmtVersion | cut -d'.' -f1)
gmtMinor=$(echo $gmtVersion | cut -d'.' -f2)
gmtPoint=$(echo $gmtVersion | cut -d'.' -f3)

#Anything before but not including GMT 5.4 is considered "legacy"
legacyGMT=true
if [ $gmtMajor -ge 6 ]; then
    legacyGMT=false
elif [ $gmtMajor -eq 5 ] && [ $gmtMinor -ge 4 ]; then
    legacyGMT=false
elif [ $gmtMajor -eq 5 ] && [ $gmtMinor -eq 3 ] && [ $gmtPoint -ge 2 ]; then
    legacyGMT=false
fi

if [ "$legacyGMT" = true ]; then
    echo "gmt blockmean -bi4d $gmtInputFile -I$spacing -R$bounds -V -Sn -Es -Wo -bo >$xyz"
    gmt blockmean -bi4d $gmtInputFile -I$spacing -R$bounds -V -Sn -Es -Wo -bo >$xyz
    nCols="7d"
else
    echo "gmt blockmean -bi4d $gmtInputFile -I$spacing -R$bounds -V -Sn -Es -bo >$xyz"
    gmt blockmean -bi4d $gmtInputFile -I$spacing -R$bounds -V -Sn -Es -bo >$xyz
    nCols="6d"
fi


if [ "$fgmt" ]; then 
    echo "gmt xyz2grd $xyz -G$xgrd -Dx/y/N -bi$nCols -I$spacing -R$bounds $fgmt"
    gmt xyz2grd $xyz -G$xgrd -Dx/y/N -bi$nCols -I$spacing -R$bounds $fgmt
    echo "gmt xyz2grd $xyz -G$ygrd -Dx/y/std -bi$nCols -i0,1,3 -I$spacing -R$bounds $fgmt"
    gmt xyz2grd $xyz -G$ygrd -Dx/y/std -bi$nCols -i0,1,3 -I$spacing -R$bounds $fgmt 
else
    echo "gmt xyz2grd $xyz -G$xgrd -Dx/y/N -bi$nCols -I$spacing -R$bounds "
    gmt xyz2grd $xyz -G$xgrd -Dx/y/N -bi$nCols -I$spacing -R$bounds 
    echo "gmt xyz2grd $xyz -G$ygrd -Dx/y/std -bi$nCols -i0,1,3 -I$spacing -R$bounds" 
    gmt xyz2grd $xyz -G$ygrd -Dx/y/std -bi$nCols -i0,1,3 -I$spacing -R$bounds  
fi



# convert GMT's output (netcdf format) to FITS
#NetCDF2FITS --ncFile $grd --fitFile $fits