#!/bin/bash

# ==================================================================================================
# Purpose: Calculate the grid size of the given corner coordinates.
# 
# Input parameters:
#     1. The latitude binary file. The following filename convention is expected of the given file:
#        latitude_<nr columns>_<nr rows>.dat. It is assumed that the binary file values are 32-bit 
#        values. Thus the size of the file is (nr columns) * (nr rows) * 4. 
#     2. The longitude binary file. The same filename convention is expected for the longitude 
#        binary file as with the latitude binary file.
#     3. The projection type. This will be either equal to Normal, Polar-North or Polar-South. 
#     4. The resolution of the grid. This must be in meters. 
#     5. The output file in which the grid information is archived in. 
# 
# The output file will have the following fields:
# UL_LON=value # The Upper-Left longitude in degrees
# UL_LAT=value # The Upper-Left latitude in degrees
# LR_LON=value # The Lower-Right longitude in degrees
# LR_LAT=value # The Lower-Right latitude in degrees
#
# NOTES:
# This script use the "cs2cs" utility of the package proj.4. proj.4 can be downloaded from 
# http://trac.osgeo.org/proj
# 
# Also note that script depends on two IDL scripts namely bin_geo_ul_lr.pro and calc_grid_size.pro.
# 
# Author: Willem.Marais@ssec.wisc.edu
# Date: 09/14/2012
# ==================================================================================================

# Check if the IDL interpreter is available
if [ -z $IDL_CMD ]; then
    idl_path_str=$(which idl)
    if [ $? -ne 0 ]; then
        echo "(ERROR) Set the environmental variable IDL_CMD equal to the IDL interpreter path."
        exit 1
    fi
    
    IDL_CMD=$idl_path_str
fi

usage_str=$(basename $0)" <Latitude binary> <Longitude binary> <Projection type: CYLINDRICALEQUIDISTANT, POLARSTEREOGRAPHIC-North or POLARSTEREOGRAPHIC-South> <Resolution> <Output file>"
if [ $# -ne 5 ]; then
    echo "(ERROR) Usage: "$usage_str
    exit 1
fi

# --------------------------------------------------------------------------------------------------
# Set variables used throughout this script
# --------------------------------------------------------------------------------------------------
lat_bin_file=$1
lon_bin_file=$2
proj_name_str=$3
resolution=$4
grid_info_txt=$5

# Define the proj4 projection strings
proj4_from_str="+proj=latlong +datum=WGS84"
case $proj_name_str in 
    "CYLINDRICALEQUIDISTANT" )
        # Projection script when used for ms2gt0.5
        #proj4_to_str="+proj=eqc +datum=WGS84"
        # Projection script when used for ms2gt0.22
        proj4_to_str="+proj=eqc +a=6371007.181 +b=6371007.181"
        xy_offset_coord="-180.0 -90.0"
    ;;
    
    "POLARSTEREOGRAPHIC-North" )
        proj4_to_str="+proj=stere +lat_0=90  +a=6371007.181 +b=6371007.181"
        xy_offset_coord="225.0 +0.0"
    ;;

    "POLARSTEREOGRAPHIC-South" )
        proj4_to_str="+proj=stere +lat_0=-90 +a=6371007.181 +b=6371007.181"
        xy_offset_coord="-225.0 +0.0"
    ;;
    
    * )
        echo "(ERROR) Unknow projection type is given. It should be either 'Normal' or 'Polar-South, Polar-North'"
        exit 1
    ;;
esac

# --------------------------------------------------------------------------------------------------
# Get the Upper-Left and Lower-Right coordinates from the latitude and longitude binary files
# --------------------------------------------------------------------------------------------------
nr_cols_int=$(basename $lat_bin_file ".bin" | cut -d"_" -f 2)
nr_rows_int=$(basename $lat_bin_file ".bin" | cut -d"_" -f 3)
ul_lr_file_str="ul_lr_info.txt"

cat > runidl << EOF
bin_geo_ul_lr, '$lat_bin_file', '$lon_bin_file', $nr_cols_int, $nr_rows_int, '$ul_lr_file_str'
exit
EOF

if [ "$SCRIPT_DEBUG" == "YES" ]; then
    cat runidl
    $IDL_CMD runidl &> idl_output.txt
    idl_rtn=$?
    cat idl_output.txt
else
    $IDL_CMD runidl &> idl_output.txt
    idl_rtn=$?
fi

ul_lon=$(grep "UL_LON" $ul_lr_file_str | cut -d"=" -f 2)
ul_lat=$(grep "UL_LAT" $ul_lr_file_str | cut -d"=" -f 2)
lr_lon=$(grep "LR_LON" $ul_lr_file_str | cut -d"=" -f 2)
lr_lat=$(grep "LR_LAT" $ul_lr_file_str | cut -d"=" -f 2)

# --------------------------------------------------------------------------------------------------
# Project latitude and longitude coordinates
# --------------------------------------------------------------------------------------------------
cs2cs_exe_str="$MIPS_HOME/bin/cs2cs $proj4_from_str +to $proj4_to_str"
xy_oul_str=$(echo "$ul_lon $ul_lat" | $cs2cs_exe_str | awk '{print $1 ", " $2}')
xy_olr_str=$(echo "$lr_lon $lr_lat" | $cs2cs_exe_str | awk '{print $1 ", " $2}')
xy_offset_str=$(echo "$xy_offset_coord" | $cs2cs_exe_str | awk '{print $1 ", " $2}')

# --------------------------------------------------------------------------------------------------
# Calulate the grid size
# --------------------------------------------------------------------------------------------------
calc_grid_size_outf="calc_grid_size_output.txt"
cat > runidl << EOF
calc_grid_size, $xy_oul_str, $xy_olr_str, $xy_offset_str, $resolution, '$calc_grid_size_outf'
exit
EOF

if [ "$SCRIPT_DEBUG" == "YES" ]; then
    cat runidl
    $IDL_CMD runidl &> idl_output.txt
    idl_rtn=$?
    cat idl_output.txt
else
    $IDL_CMD runidl &> idl_output.txt
    idl_rtn=$?
fi

grep --ignore-case "Error" idl_output.txt &> /dev/null
if [ $? -eq 0 ] || [ $idl_rtn -ne 0 ]; then
    echo "(ERROR) IDL script calc_grid_size.pro failed"
    cat idl_output.txt
    exit 1
fi

# Collect the projected adjusted coordinates and the the calculate grid size
adj_ul_x=$(grep ADJ_UL_X $calc_grid_size_outf | cut -d"=" -f 2)
adj_lr_x=$(grep ADJ_LR_X $calc_grid_size_outf | cut -d"=" -f 2)
adj_ul_y=$(grep ADJ_UL_Y $calc_grid_size_outf | cut -d"=" -f 2)
adj_lr_y=$(grep ADJ_LR_Y $calc_grid_size_outf | cut -d"=" -f 2)
columns=$(grep COLS $calc_grid_size_outf | cut -d"=" -f 2)
rows=$(grep ROWS $calc_grid_size_outf | cut -d"=" -f 2)

# --------------------------------------------------------------------------------------------------
# Inverse project the adjusted coordinates to get the longitude and latitude values
# --------------------------------------------------------------------------------------------------
cs2cs_exe_str="$MIPS_HOME/bin/cs2cs $proj4_to_str +to $proj4_from_str -f %.6f"
adj_lon_lat_oul_str=$(echo "$adj_ul_x $adj_ul_y" | $cs2cs_exe_str | awk '{print $1" "$2}')
adj_lon_lat_olr_str=$(echo "$adj_lr_x $adj_lr_y" | $cs2cs_exe_str | awk '{print $1" "$2}')

echo "COLUMNS=$columns" > $grid_info_txt
echo "ROWS=$rows" >> $grid_info_txt
echo "UL_LON="$(echo $adj_lon_lat_oul_str | cut -d" " -f 1) >> $grid_info_txt
echo "UL_LAT="$(echo $adj_lon_lat_oul_str | cut -d" " -f 2) >> $grid_info_txt
echo "LR_LON="$(echo $adj_lon_lat_olr_str | cut -d" " -f 1) >> $grid_info_txt
echo "LR_LAT="$(echo $adj_lon_lat_olr_str | cut -d" " -f 2) >> $grid_info_txt

# --------------------------------------------------------------------------------------------------
# Remove files
# --------------------------------------------------------------------------------------------------
# rm -f idl_output.txt runidl $ul_lr_file_str $calc_grid_size_outf flsave.pro