/***************************************************************************
 *   Copyright (C) 2010 by Nicolas PASCAL   *
 *   nicolas.pascal@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.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#include "pmfiledata.h"

const float PMFileData::colocation_tolerance = 0.01;
const int16 PMFileData::sz_grid_y = 3240;
const int16 PMFileData::sz_grid_x = 1296;
const string PMFileData::latitude_name  = "Latitude";
const string PMFileData::longitude_name = "Longitude";

PMFileData::PMFileData ( const string &name , const string & mode ) :
    FileData ( name, mode ), SatelliteFileData ( name, mode ), HDFFileData ( name, mode ) {
    init();
}
void PMFileData::init() {
    lat_data  = NULL;
    lon_data  = NULL;
    time_data = NULL;
    irow_shift_data = NULL;
    product   = "";
    version   = "";
    time_coverage = -1.;
    orbit_nb = -1;

    // parse the file name to know the product, the acquisition start time...
    parse_filename ( get_tail ( get_name() ) );

    load_hdf_file();

    // set the attibutes
    set_lat_lon_index_max();
    set_time_coverage();
    read_node_lon ();

    free_hdf_file();
}

void PMFileData::read_node_lon() {
    Hdf_attr file_attr = hdf_file->get_attr ( "CenterLongitudeDeg" );
    file_attr.get_values ( &node_lon );
}

PMFileData::~PMFileData() {
    free_geolocation_data();
    free_hdf_metadata();
    free_hdf_file();
}
void PMFileData::set_time_coverage( ) {
    //-- Read the end time in the hdf file attributes
    // date has a format like "2009-01-01:05:00:16.460Z" "YYYYMMDDHHMMSS" -> 24 chars
    char schar_end_time[24+1];
    fill_n ( schar_end_time, 0, 25 );
    hdf_file->get_attr("EndingAcquisitionDate").get_values ( schar_end_time );
    string s_endtime(schar_end_time);
    //-- Build a date with this end time
    Date end_time;
    end_time.set_date_str ( s_endtime.substr ( 0, 19 ) ,"%Y-%m-%d:%H:%M:%S" );
    //-- Substrack end and start time to compute the time coverage
    time_coverage=end_time.get_epoch_time()-date->get_epoch_time();
}
void PMFileData::free_geolocation_data() {
    delete[] lat_data;
    lat_data=NULL;
    delete[] lon_data;
    lon_data=NULL;
    delete[] time_data;
    time_data=NULL;
    delete[] irow_shift_data;
    irow_shift_data = NULL;
}

void PMFileData::load_geolocation_data() {
    if (!is_geolocation_data_loaded()) {
        // hdf_file is needed, but the method assume to conserve the caller's state
        bool hdf_file_already_loaded = is_hdf_file_loaded();
        if (!hdf_file_already_loaded)
            load_hdf_file();
        //-- Read lat,lon,time data
        if (lat_data==NULL)
            lat_data=(float*)read_data(lat_data,latitude_name.c_str());
        if (lon_data==NULL)
            lon_data=(float*)read_data(lon_data,longitude_name.c_str());
        // time data unused for this product
        if ( irow_shift_data == NULL )
            irow_shift_data = ( int16 * ) read_data ( irow_shift_data, "ISinus_Grid_Row_Shift_Full_Resolution" );
        if (!hdf_file_already_loaded)
            free_hdf_file();
    }
}

void PMFileData::set_lat_lon_index_max() {
    lat_lon_index_max [0] = PMFileData::sz_grid_y;
    lat_lon_index_max [1] = PMFileData::sz_grid_y;
}

void PMFileData::parse_filename( const string& short_filename ) {
    check_filename(short_filename); // throw a parse_CLOUDSAT_filename_error if the filename doesn't matches the pattern

    using namespace MyTools; // for has_type function

    int sz = short_filename.size();
    size_t i = short_filename.find_first_of ( ".", 0 );

    // Product name
    product = short_filename.substr ( 0, i );
    // timestamp : the first 14 characters are digits and represents the orbit start date (in UTC time)
    ++i;
    string s_date = short_filename.substr ( i, 14 );
    date->set_date_str ( s_date,"%Y%m%d%H%M%S");
    // orbit number
    i += 17;
    orbit_nb = to_num < int > ( short_filename.substr ( i, 6 ) );
    // version
    i += 7;
    version = short_filename.substr( i, sz - i - 4 );
}

bool PMFileData::check_filename( const string & short_filename )
{
    int sz = short_filename.size();
    using namespace MyTools; // for has_type function

    size_t i = short_filename.find_first_of ( ".", 0 );

    string s_product = short_filename.substr ( 0, i );
    if ( s_product != string ( "PM01_L2" ) &&
         s_product != string ( "PM02_L2" ) ) {
      invalid_filename e(__FILE__,__LINE__,short_filename,"POLDER-MODIS");
      throw e;
    }

    // timestamp : the first 14 characters are digits and represents the orbit start date (in UTC time)
    ++i;
    string s_date = short_filename.substr ( i, 16 );
    if ( ! has_type < unsigned long long > ( s_date )) {
      invalid_filename e(__FILE__,__LINE__,short_filename,"POLDER-MODIS");
      throw e;
    }
    // orbit number
    i += 17;
    string s_orbit_nb=short_filename.substr( i, 6 );
    if (!has_type<unsigned long>(s_orbit_nb)) {
      invalid_filename e(__FILE__,__LINE__,short_filename,"POLDER-MODIS");
      throw e;
    }

    // extension
    string s_extension = short_filename.substr ( sz - 4, 4 );
    if ( s_extension != string ( ".hdf" ) ) {
      invalid_filename e(__FILE__,__LINE__,short_filename,"POLDER-MODIS");
      throw e;
    }

    return true;
}
const bool PMFileData::contain_location(const float &lat,const float &lon, const double &tolerance){
  int nearest_idx[]={-1,-1};
  return get_index(lat,lon,nearest_idx,colocation_tolerance);
}
const bool PMFileData::contain_data( const float & lat, const float & lon, const double & time, const double &colocation_tolerance ) {
  return contain_location(lat,lon,colocation_tolerance)&&contain_time(time);
}

const bool PMFileData::get_index( const float & lat, const float & lon, int igrid[], const float colocation_tolerance ) {
    assert( irow_shift_data != NULL );
    bool already_loaded_data=is_geolocation_data_loaded();
    if (!already_loaded_data)
        load_geolocation_data();

    fill_n ( igrid, 2, -1 );

    float32 newlon = lon - node_lon;
    if ( newlon > 180.0  ) newlon -= 360. ;
    if ( newlon < -180.0 ) newlon += 360. ;

    geolocation_to_grid ( lat, newlon, igrid );

//     printf ( "\tigrid [%ld,%ld]\n", igrid[0], igrid[1] );

    if ( ( igrid [0] != -1 ) && ( igrid [1] != -1 ) ) {
        /* change Fortran indexing ( from 1 ) to C way ( from 0 ) */
        igrid [0] -= 1, igrid [1] -= 1;
        /* Translation in column to get the column number in the required subset */
        igrid [1] -= 2 * PMFileData::sz_grid_x;
        if ( igrid [1] < 0 || igrid [1] > PMFileData::sz_grid_x )
            // invalidate
            fill_n ( igrid, 2, -1 );
    }

//     printf ( "\tigrid [%ld,%ld]\n", igrid[0], igrid[1] );

    if (!already_loaded_data)
        free_geolocation_data();

    return ( igrid [0] != -1 && igrid [1] != -1 );
}

const bool PMFileData::get_index( const float & lat, const float & lon, int &nearest_pix_idx, const float colocation_tolerance ) {
    nearest_pix_idx = -1;

    int idata[] = { -1, -1 };
    get_index ( lat, lon, idata );

    if ( idata [0] != -1 && idata [1] != -1 )
        nearest_pix_idx = idata [0] * sz_grid_x +  idata [1];

    return ( nearest_pix_idx != -1 );
}

void PMFileData::geolocation_to_grid( const float * lat_lon, int * grid_idx )
{
    assert ( grid_idx != NULL );
    assert ( lat_lon  != NULL );
    int nlines = PMFileData::sz_grid_y;
    float lat = lat_lon[0];
    float lon = lat_lon[1];
    if (lat_lon[0]>=-90. && lat_lon[0]<=90. && lat_lon[1]>=-180. && lat_lon[1]<=180.) {
        int iline = lround ( 18. * ( 90. - lat ) + 0.5 );
        int ni    = lround ( nlines * sin ( M_PI * ( iline - 0.5 ) / nlines ) );
        int icol  = lround ( nlines + 0.5 + ni * lon / 180. );
        grid_idx [0] = iline, grid_idx [1] = icol;
    } else {
        fill_n ( grid_idx, 2, -1 );
    }
}

void PMFileData::grid_to_geolocation( const int * grid_idx, float * lat_minmax, float * lon_minmax )
{
    int grid_factor = 1;
    // grid middle pixel lat,lon
    double lat = -9999.;
    double lon = -9999.;
    double cell_side=static_cast<double>(grid_factor)/18.;
    double delta = 0.5*cell_side; // 2*delta = size of a grid cell side
    double nb_grid_lat=3240./static_cast<double>(grid_factor);
    double nb_grid_lon=6480./static_cast<double>(grid_factor);
    int ni=-1;
    int line=grid_idx[0];
    int col= grid_idx[1];
    if (grid_idx!=NULL) {
        if (lat_minmax!=NULL) {
            if (line>0 && line<=(nb_grid_lat) ) { // lat range can be computed
                lat=90.-(static_cast<double>(line)-.5)*cell_side;
                ni=lround((nb_grid_lat)*cos(radians(lat)));
                lat_minmax[0]=lat-delta;
                lat_minmax[1]=lat+delta;
            } else {
                lat_minmax[0]=lat_minmax[1]=-9999.;
            }
            if (lon_minmax!=NULL) {
                if (ni!=-1 && col>0 && col<=nb_grid_lon) { // lon range can be computed
                    lon=(180./ni)*(static_cast<double>(col)-nb_grid_lon/2.-.5);
                    //lon=(180./ni)*(static_cast<double>(col)+1./ni);
                    double delta_lon=180./ni;
                    lon_minmax[0]=lon-delta_lon;
                    lon_minmax[1]=lon+delta_lon;
                } else {
                    lon_minmax[0]=lon_minmax[1]=-9999.;
                }
            }
        }
    }
}