/* main.cpp */

/*
remap
Copyright (C) 2006 Fabrice Ducos, fabrice.ducos@icare.univ-lille1.fr

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.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#include <iostream>
#include <cstdlib>

#include <unistd.h> /* getopt */
#include <cmath>
#include <exception>
#include "common.h"
#include "parse_argument.h"
#include "filetypes.h"
#include "grid.h"
#include "reproject.h"
#include "hdf_utils.h"
#include "radiances_to_temperatures.h"
#include "VEquiRect.h"

#define DEFAULT_OUTPUT_FILENAME "remap.out.hdf"

using namespace std;

bool g_verbose = false;
bool g_extended_output = false;
bool g_brillance_temperatures = false;
char output_filename[STRING_MAXLEN + 1];
distance_type g_distance = 1.; // distance of colocalization, in kilometers
time_type g_max_dtime = INFINITY; // time of colocalization, in seconds

bool user_defined_input = false;
coord_type g_equirect_latitude_orig = DEFAULT_COORD_FILL_VALUE;
coord_type g_equirect_longitude_orig = DEFAULT_COORD_FILL_VALUE;
distance_type g_equirect_resolution = 0.; // resolution of user defined equirectangular grids, in kilometers
int g_equirect_nrows = 0;
int g_equirect_ncols = 0;

distance_type g_distance_fill_value = DEFAULT_DISTANCE_FILL_VALUE;
time_type g_time_fill_value = DEFAULT_TIME_FILL_VALUE;
data_type g_data_fill_value = DEFAULT_DATA_FILL_VALUE;
coord_type g_coord_fill_value = DEFAULT_COORD_FILL_VALUE;

float64 g_output_slope = 0.01;
float64 g_output_offset = 0.;

void usage() {

  cerr << "usage: " APPNAME " [OPTIONS] dataset_ref[channel_ref]@file_ref dataset1[channel1]@file1 dataset2[channel2]@file2..." << endl << endl;
  cerr << "OPTIONS:" << endl;
  cerr << "  -c slope,offset                  calibration factors (default: slope = " << g_output_slope << " offset = " << g_output_offset << ")" << endl;
  cerr << "  -d <distance>                    distance of colocalization, in km (default is " << g_distance << " km)" << endl;
  cerr << "  -h                               help (displays supported types of files and implemented user's defined projections)" << endl;
  cerr << "  -p <user defined grid>           specifies a user defined projection instead of a file defined projection (use -h for details)" << endl;
  cerr << "  -m <spec. for missing value>     specifies values for missing coordinates, time or data (use -h for details)" << endl;
  cerr << "  -o <filename>                    output filename (default is " DEFAULT_OUTPUT_FILENAME ")" << endl;
  cerr << "  -t <max_dtime>                   max value of dtimes to reference, in minutes" << endl;
  cerr << "                                   (may be inf, for infinity, or a strictly positive real number, default is " << g_max_dtime/60. << ")" << endl;
  cerr << "  -v                               verbose mode" << endl;
  cerr << "  -x                               extended output (record delta in time and distance)" << endl;
  /* option -T still unstable: not to be documented */
  //cerr << "  -T                               converts radiances into brillance temperatures" << endl;
  cerr << endl;
  cerr << "Supported types of files are documented in the online help (option -h)" << endl;
  cerr << "version of " APPNAME ": " << VERSION << endl;
  exit (EXIT_FAILURE);
}

/* options */
void parse_options(int *argc, char **argv[]) {
  int option;
  const char *options_list = "c:d:hp:m:o:r:t:vx";
  
  while ((option = getopt(*argc, *argv, options_list)) != -1 ) {
    switch(option) {
    case 'c' :
      {
	if (sscanf(optarg, "%lf,%lf", &g_output_slope, &g_output_offset) == 2) {
	  /* do nothing */
	}
	else {
	  cerr << APPNAME ": invalid argument for -c option: " << optarg << endl;
	  exit (EXIT_FAILURE);
	}
      }
      break;
    case 'd' : 
      {
	/* distance of colocalization */
	g_distance = atof(optarg);
	if (g_distance <= 0.) {
	  cerr << APPNAME ": invalid value for distance of colocalization (should be positive)" << endl;
	  exit (EXIT_FAILURE);
	}
      }
      break;
    case 'h' :
      {
	 print_supported_filetypes();
	 cout << endl;
	 cout << "Option -m (may be used several times, for each type of missing value):" << endl;
	 cout << "  coord=<value for missing coordinates>      simple precision floating point number (default is " << g_coord_fill_value << ", nan = not a number, inf = infinity)" << endl;
	 cout << "  distance=<value for missing distances>     simple precision floating point number (default is " << g_distance_fill_value << ", nan = not a number, inf = infinity)" << endl;
	 cout << "  time=<value for missing times>             double precision floating point number (default is " << g_time_fill_value  << ", nan = not a number, inf = infinity)" << endl;
	 cout << "  data=<value for missing data>              16-bit signed integer, between " << DATA_TYPE_MIN << " and " << DATA_TYPE_MAX << " (default is " << g_data_fill_value << ")" << endl;
	 cout << endl;
	 cout << "Option -p:" << endl;
	 cout << "  bbox=lat1,lon1,lat2,lon2,resolution         bounding box for an equirectangular projection: upper left corner (1), lower right corner (2)" << endl;
	 cout << "  tbox=lat,lon,weight,height,resolution       target box for an equirectangular projection (lat,lon are the center's coordinates)" << endl;
	 cout << "  lcbox=lat,lon,weight,height,resolution      equirectangular projection specified by the upper left corner's coordinates, weight, height" << endl;
	 cout << endl;
	 cout << "(latitudes and longitudes are to be given in degrees, weights and heights in pixels, resolution in km/pixel)" << endl;
	 exit (EXIT_SUCCESS);
      }
      break;
    case 'm' :
      {
	if (sscanf(optarg, "data=%hud", &g_data_fill_value) == 1) {
	  /* do nothing */
	}
	else if (sscanf(optarg, "coord=%f", &g_coord_fill_value) == 1) {
	  /* do nothing */
	}
	else if (sscanf(optarg, "distance=%f", &g_distance_fill_value) == 1) {
	  /* do nothing */
	}
	else if (sscanf(optarg, "time=%lf", &g_time_fill_value) == 1) {
	  /* do nothing */
	}
	else {
	  cerr << APPNAME ": invalid specification for -m option: " << optarg << " (will be ignored)" << endl;
	}
      }
      break;
    case 'o' :
      {
	/* specifies output filename */
	strncpy(output_filename, optarg, STRING_MAXLEN);
      }
      break;
    case 'p' :
      {
	double lat, lon;
	double lat1, lon1;
	double lat2, lon2;
	int weight, height;
	double resolution;

	if (sscanf(optarg, "bbox=%lf,%lf,%lf,%lf,%lf", &lat1, &lon1, &lat2, &lon2, &resolution) == 5) {
	  /* bounding box */

	  const coord_type stride = resolution*VEquiRect::DEGREES_PER_KM; // in degrees per pixel
	  const coord_type lat_min = min(lat1, lat2);
	  const coord_type lat_max = max(lat1, lat2);
	  const coord_type lon_min = min(lon1, lon2);
	  const coord_type lon_max = max(lon1, lon2);

	  lat = lat_max;
	  lon = lon_min;
	  height = lround((lat_max - lat_min)/stride);
	  weight = lround((lon_max - lon_min)/stride);

	}
	else if (sscanf(optarg, "tbox=%lf,%lf,%d,%d,%lf", &lat, &lon, &weight, &height, &resolution) == 5) {
	  /* target box */
	  cerr << APPNAME << ": target box not yet implemented, sorry" << endl;
	  exit (EXIT_FAILURE);
	}
	else if (sscanf(optarg, "lcbox=%lf,%lf,%d,%d,%lf", &lat, &lon, &weight, &height, &resolution) == 5) {
	  /* left corner box */
	}
	else {
	  cerr << APPNAME ": invalid specification for option -p" << endl;
	  exit (EXIT_FAILURE);
	}
	if ((MIN_LAT <= lat && lat <= MAX_LAT) == false) {
	  cerr << APPNAME ": equirectangular projection: left corner latitude (in degrees) out of bounds: " << lat << endl;
	  exit (EXIT_FAILURE);
	}
	
	if ((MIN_LON <= lat && lat <= MAX_LON) == false) {
	  cerr << APPNAME ": equirectangular projection: left corner longitude (in degrees) out of bounds: " << lat << endl;
	  exit (EXIT_FAILURE);
	}

	if ((resolution > 0) == false) {
	  cerr << APPNAME ": equirectangular projection: resolution (in km/pixel) not positive: " << resolution << endl;
	  exit (EXIT_FAILURE);
	}

	if ((weight > 0) == false) {
	  cerr << APPNAME ": weight (in pixels) not positive: " << weight << endl;
	  exit (EXIT_FAILURE);
	}

	if ((height > 0) == false) {
	  cerr << APPNAME ": height (in pixels) not positive: " << height << endl;
	  exit (EXIT_FAILURE);
	}
	
	user_defined_input = true;
	g_equirect_latitude_orig = lat;
	g_equirect_longitude_orig = lon;
	g_equirect_resolution = resolution;
	g_equirect_nrows = height;
	g_equirect_ncols = weight;
	g_max_dtime = INFINITY; // disable time filtering
      }
      break;
    case 'T' :
      {
	g_brillance_temperatures = true;
      }
      break;
    case 't' :
      {
	if (strncasecmp(optarg, "inf", STRING_MAXLEN) == 0) {
	  g_max_dtime = INFINITY;
	}
	else {
	  g_max_dtime = atof(optarg)*60.;
	  if (g_max_dtime <= 0.) {
	    cerr << APPNAME ": invalid value for max_dtime: " << optarg << " (shoud be positive, in minutes)" << endl;
	    exit (EXIT_FAILURE);
	  }
	}
      }
      break;
    case 'v' :
      {
	/* enables verbose mode */
	g_verbose = true;
      }
      break;
    case 'x' :
      { 
	g_extended_output = true;
      }
      break;
    case '?' :
      {
	fprintf(stderr, APPNAME ": option %c : unknown\n\n",optopt);
	usage();
      }
      break;
    default :
      {
	fprintf(stderr, APPNAME ": parse options : fatal error\n");
	exit(8);
      }
      break;
    } /* switch */
  } /* while */
  *argc -= optind - 1;
  *argv += optind - 1;
}

int main(int argc, char *argv[]) {

  int err;
  grid_type target_grid;
  grid_type previous_src_grid;
  int first_iarg_to_reproject = 0;

  parse_options(&argc, &argv);

  if (argc < 2) usage();

  if (! *output_filename) {
    strncpy(output_filename, DEFAULT_OUTPUT_FILENAME, STRING_MAXLEN);
  }
  
  memset(&target_grid, 0, sizeof(target_grid));

  if (user_defined_input) {
    target_grid.file[0] = '\0'; // redundant but harmless
    first_iarg_to_reproject = 1;
  }
  else {
    err = parse_argument(1, argv, &target_grid);
    if (err != 0) {
      Debug(cerr << APPNAME << ": parse_argument failed" << endl;);
      exit (EXIT_FAILURE);
    }
    first_iarg_to_reproject = 2;
  }

  target_grid.is_target = true;
  err = load_grid(&target_grid);
  if (err != 0) {
    cerr << APPNAME << ": failed to load reference " << argv[1] << ", abort" << endl;
    exit (EXIT_FAILURE);
  }
 
  strncpy(target_grid.file, output_filename, STRING_MAXLEN);
  hdf_create_empty_file(output_filename);

  err = save_latlontime(&target_grid);
  if (err != 0) {
    cerr << APPNAME << ": failed to save lat,lon,time in " << output_filename << ", abort" << endl;
    exit (EXIT_FAILURE);
  }

  if (first_iarg_to_reproject == 2) {
    assert(target_grid.data);
    if (g_verbose == true) {
      cerr << APPNAME << ": reference: " << argv[1] << endl;
    }

    if (g_brillance_temperatures == true && target_grid.content == GRID_CONTENT_RADIANCES && target_grid.wavelength > 0.) {
      convert_radiances_to_temperatures(&target_grid);
    }
    
    err = save_grid(&target_grid, false /* don't save delta-time and delta-distance which are irrelevant for the reference */);
    if (err != 0) {
      cerr << APPNAME << ": failed to save entry " << argv[1] << ", abort" << endl;
      exit (EXIT_FAILURE);
    }
  }
  
  memset(&previous_src_grid, 0, sizeof(previous_src_grid));
  copy_grid_footprint(&previous_src_grid, &target_grid);

  assert(first_iarg_to_reproject >= 1);

  for (int iarg = first_iarg_to_reproject ; iarg < argc ; iarg++) {
    try {
      grid_type src_grid;
      
      memset(&src_grid, 0, sizeof(src_grid));
      err = parse_argument(iarg, argv, &src_grid);
      if (err != 0) {
	Debug(cerr << APPNAME << ": parse_argument failed" << endl;);
	exit (EXIT_FAILURE);
      }
      
      if (g_verbose == true) {
	cout << APPNAME << ": processing of " << argv[iarg] << endl;
      }
      
      err = load_grid(&src_grid);
      
      if (err != 0) {
	cerr << APPNAME << ": failed to load entry " << argv[iarg] << ", ignore it" << endl;
	continue;
      }
      
      if (have_same_grid_coordinates(&previous_src_grid, &src_grid) == false) {
	reproject_coordinates(&src_grid, &target_grid);
      }

      reproject_measures(&src_grid, &target_grid);

      if (g_brillance_temperatures == true && target_grid.content == GRID_CONTENT_RADIANCES && target_grid.wavelength > 0.) {
	convert_radiances_to_temperatures(&target_grid);
      }

      //err = save_grid(&target_grid, g_extended_output && ! have_same_grid_coordinates(&previous_src_grid, &src_grid));
      err = save_grid(&target_grid, g_extended_output);
      if (err != 0) {
	cerr << APPNAME << ": failed to save entry " << argv[iarg] << ", ignore it" << endl;
	continue;
      }
      copy_grid_footprint(&previous_src_grid, &src_grid);
      destroy_grid(&src_grid);
    } // try
    catch (const char *err_msg) {
      cerr << APPNAME ": " << err_msg << endl;
      continue;
    }
    catch (const exception &e) {
      cerr << APPNAME ": " << e.what() << endl;
      continue;
    }
    catch (...) {
      cerr << APPNAME ": unexpected exception" << endl;
      exit (EXIT_FAILURE);
    }
  } // for (int iarg)
 
  destroy_grid(&target_grid); 
  return EXIT_SUCCESS;

}