/* VModis.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 <cassert>

#include "debug.h"
#include "common.h"
#include "VModis.h"
#include "VModis_latlon_resolve.h"
#include "VModis_interpol.h"
#include "hdffiledata.h"

const double
VModis::wavelengths[] =  { 0.00, 3.78, 3.97, 3.98, 4.06, 4.45, 4.53, 6.78, 7.34, 8.55, 9.72, 11.02, 12.04, 13.36, 13.68, 13.92, 14.21 };

// interface to J. Descloitres' interpolation code for Modis geolocation data
void
VModis::interpol(int sample_nrows, int sample_ncols, const coord_type *sample_lat, const coord_type *sample_lon,
		 int full_nrows, int full_ncols, coord_type *full_lat, coord_type *full_lon) {

  PRODUCT product;
  int georow1, georow2;
  double t;
  double *lat3 = new double[sample_ncols];
  double *lon3 = new double[sample_ncols];

  product.Nl_geo = sample_nrows;
  product.Np_geo = sample_ncols;
  product.sds.Nl = full_nrows;
  product.sds.Np = full_ncols;
  product.projparam.type = INTERP_SWATH;

  const int rowsperscan = (int) (10*product.sds.Np/1354. + 0.5);
  const int georowsperscan = (int) (10 * product.Np_geo / 1354. + 0.5);
  const int georowsamplingstep = 10/georowsperscan;
  
  get_step_and_offset(product.sds.Nl, product.Nl_geo, &product.rowoffset, &product.rowstep, rowsperscan, georowsamplingstep, FALSE /* crosstrack */);
  get_step_and_offset(product.sds.Np, product.Np_geo, &product.coloffset, &product.colstep, rowsperscan, georowsamplingstep, TRUE /* crosstrack */);

  for (int irow = 0 ; irow < full_nrows ; ++irow) {
    const int iscan = irow/rowsperscan;
    get_georows(irow, &product, iscan, &georow1, &georow2, &t);
    
    assert(0 <= georow1 && georow1 < sample_nrows);
    assert(0 <= georow2 && georow2 < sample_nrows);

    const coord_type *lon1 = &sample_lon[georow1*sample_ncols];
    const coord_type *lon2 = &sample_lon[georow2*sample_ncols];
    const coord_type *lat1 = &sample_lat[georow1*sample_ncols];
    const coord_type *lat2 = &sample_lat[georow2*sample_ncols];
    get_virtual_row(lon1, lat1, lon2, lat2, t, product.Np_geo, lon3, lat3);
    resample_row_geolocation(INTERP_CENTERS, lon3, lat3, &product, &lon_[irow*full_ncols], &lat_[irow*full_ncols]);

  }

  delete[] lat3;
  delete[] lon3;
  
}

double
VModis::ichannel_to_wavelength(int ichannel) const {

  assert(ichannel >= 0 && ichannel < sizeof(wavelengths)/sizeof(*wavelengths));

  return wavelengths[ichannel];

}

VModis::VModis(const char *filename, const char *dataset, int ichannel) 
 try : VHdf(filename, dataset, ichannel)
{

  PDEBUG;

  if (strstr(dataset, "EV_1KM_Emissive") != NULL) {
    content_ = GRID_CONTENT_RADIANCES;
  }
  else {
    content_ = GRID_CONTENT_DEFAULT;
  }

  wavelength_ = ichannel_to_wavelength(ichannel_);
  read_lat_lon_time();

}
catch (const std::exception &e) {
  cerr << APPNAME ": " << __FILE__ << ": " << __LINE__ << ": " << e.what() << endl;
}

void VModis::get_calibration(double &slope, double &offset) const {
  PDEBUG;

  float32 slope32;
  float32 offset32;
  float64 slope64;
  float64 offset64;

  slope = 1.;
  offset = 0.;

  try {
    HDFFileData f_data(data_filename_, "r");
    Hdf_sds sds = f_data.get_hdf_file()->get_sds(sds_data_);

    if (strstr(sds_data_, "EV_1KM_Emissive")  != NULL || strstr(sds_data_, "EV_Band26") != NULL) {
      Hdf_attr scales_attr = sds.get_attribute("radiance_scales");
      Hdf_attr offset_attr = sds.get_attribute("radiance_offsets");
      assert(ichannel_ > 0);
      scales_attr.get_value(&slope32, ichannel_ - 1);
      offset_attr.get_value(&offset32, ichannel_ - 1);
      slope = slope32;
      offset = offset32;
    }
    else {
      Hdf_attr scale_factor = sds.get_attribute("scale_factor");
      scale_factor.get_value(&slope64);
      offset64 = 0.;

      slope = slope64;
      offset = offset64;
    }
  }
  catch (bad_file &e) {
    Debug(cerr << APPNAME << ": " << __PRETTY_FUNCTION__ << ": " << e.what() << endl;);
    throw;
  }
  catch (bad_sds_name &e) {
    Debug(cerr << APPNAME << ": " << __PRETTY_FUNCTION__ << ": " << e.what() << endl;);
    throw;
  }
  catch (const char *err_msg) {
    Debug(cerr << APPNAME << ": " << __PRETTY_FUNCTION__ << ": " << err_msg << endl;);
    throw;
  }
  catch (const string& err_msg) {
    Debug(cerr << APPNAME << ": " << __PRETTY_FUNCTION__ << ": " << err_msg << endl;);
    throw;
  }
  catch (...) {
    cerr << APPNAME << ": " __FILE__ << ": " << __LINE__ << ": unexpected exception" << endl;
    throw;
  }
  
}

void VModis::read_lat_lon_time() {
  const char *vdata_name = "Level 1B Swath Metadata";
  const char *field_name = "EV Sector Start Time";

  assert(sds_lat_ && *sds_lat_);
  assert(sds_lon_ && *sds_lon_);
  assert(lat_ == NULL);
  assert(lon_ == NULL);
  assert(time_ == NULL);

  try {
    HDFFileData f_data(data_filename_, "r");
    vector<long> sample_dimensions = f_data.get_sds_dimension(sds_lat_);
    vector<long> data_dimensions = f_data.get_sds_dimension(sds_data_);
    int sample_nrows = sample_dimensions[0];
    int sample_ncols = sample_dimensions[1];
    int data_nrows;
    int data_ncols;
     
    if (data_dimensions.size() == 2) {
      data_nrows = data_dimensions[0];
      data_ncols = data_dimensions[1];
    }
    else if (data_dimensions.size() == 3) {
      data_nrows = data_dimensions[1];
      data_ncols = data_dimensions[2];
    }
    else {
      std::cerr << APPNAME << ": " << __PRETTY_FUNCTION__ << ": unexpected value for data_dimensions.size(): " << data_dimensions.size() << std::endl;
      exit (EXIT_FAILURE);
    }

    vector<long> vdata_dimensions = f_data.get_vdata_dimension(vdata_name);
    
    assert(vdata_dimensions.size() == 2);
    int nscan = vdata_dimensions[0];
    int nrows_per_scan = data_nrows/nscan;
    
    // if (data_nrows % nrows_per_scan != 0) {
    if (nrows_per_scan != 10) {
      std::cerr << APPNAME << ": " << data_filename_ << ": inconsistency between number of rows in data (" << data_nrows << ") and number of scans (" << nscan << ")" << std::endl;
      exit (EXIT_FAILURE);
    }
    
    float64 *time_buffer = NULL;
    time_buffer = static_cast<float64 *>(f_data.read_vdata(time_buffer, vdata_name, field_name));

    time_ = new time_type[data_nrows*data_ncols];
    
    for (int irow = 0 ; irow < data_nrows ; ++irow) {
      for (int icol = 0 ; icol < data_ncols ; ++icol) {
	const int iscan = irow/nrows_per_scan;
	time_[irow*data_ncols + icol] = time_buffer[iscan];
      }
    }

    free(time_buffer); time_buffer = NULL;

    float32 *sample_lat = NULL;
    float32 *sample_lon = NULL;

    sample_lat = f_data.read_data(sample_lat, "Latitude");
    sample_lon = f_data.read_data(sample_lon, "Longitude");

    lat_ = new coord_type[data_nrows*data_ncols];
    lon_ = new coord_type[data_nrows*data_ncols];
    assert(lat_);
    assert(lon_);

    interpol(sample_nrows, sample_ncols, sample_lat, sample_lon, data_nrows, data_ncols, lat_, lon_);
   
    delete[] sample_lat; sample_lat = NULL;
    delete[] sample_lon; sample_lon = NULL;

  }
  catch (bad_file &e) {
    Debug(cerr << APPNAME << ": " << __PRETTY_FUNCTION__ << ": " << e.what() << endl;);
    throw;
  }
  catch (bad_sds_name &e) {
    Debug(cerr << APPNAME << ": " << __PRETTY_FUNCTION__ << ": " << e.what() << endl;);
    throw;
  }
  catch (...) {
    cerr << APPNAME << ": " __FILE__ << ": " << __LINE__ << ": unexpected exception" << endl;
    throw;
  }
}

VModis::~VModis() {

  PDEBUG;

}