/***************************************************************************
 *   Copyright (C) 2005 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 "iirfiledata.h"

const string IIRFileData::latitude_sds_name="Latitude";
const string IIRFileData::longitude_sds_name="Longitude";
// string IIRFileData::time_sds_name="Lidar_Shot_Time";
const int IIRFileData::IIR_swath_size=69;
const float IIRFileData::colocation_tolerance=0.05;
const int IIRFileData::nb_record_max=25000;

IIRFileData::IIRFileData(const string &name,const string & mode):
    FileData(name,mode),SatelliteFileData(name,mode),HDFFileData(name,mode) {
    init();
}
void IIRFileData::init() {
    lat_data=NULL;
    lon_data=NULL;
    time_data=NULL;

    // extract the acquisition date from the file name
    string short_filename = get_tail(get_name());
    parse_filename(short_filename);
    // set the time SDS name depending of the product
    if ( product_type == L1 )
        time_sds_name = "Lidar_Shot_Time";
    else
        time_sds_name = "LIDAR_Shot_Time";

    // read some attibutes from the metadata
    load_hdf_file();
    set_lat_lon_index_max();
    set_time_coverage();
    free_hdf_file();
}
IIRFileData::~IIRFileData() {
    free_geolocation_data();
    free_hdf_file();
    free_hdf_metadata();
}
void IIRFileData::set_lat_lon_index_max() {
    Hdf_file* hdf = get_hdf_file();
    if (hdf!=NULL) {
        lat_lon_index_max[0] = hdf->get_sds(latitude_sds_name.c_str()).get_dimension(0)-1;
        lat_lon_index_max[1] = hdf->get_sds(latitude_sds_name.c_str()).get_dimension(1)-1;
    }
}
void IIRFileData::check_filename( const string& short_filename ) const {
    // check the extension
    string extension=short_filename.substr(short_filename.size()-4,4);
    if (extension!=".hdf") {
      invalid_filename e(__FILE__,__LINE__,short_filename,"IIR");
      throw e;
    }
    // check the product radix
    int first_minus_pos = short_filename.find_first_of("-");
    string radix=short_filename.substr(0,first_minus_pos);
    if (radix!="CAL_IIR_L1" &&
        radix!="CAL_IIR_L2_Track" && radix!="CAL_IIR_L2_Swath") {
      invalid_filename e(__FILE__,__LINE__,short_filename,"IIR");
      throw e;
    }
    // check the day/night flag
    char day_mode=short_filename[short_filename.size()-5];
    if (day_mode!='D' && day_mode!='N') {
      invalid_filename e(__FILE__,__LINE__,short_filename,"IIR");
      throw e;
    }
}
void IIRFileData::parse_filename( const string& short_filename ) {
    check_filename(short_filename); // throw a parse_IIR_filename_error if the filename doesn't matches the pattern
    // set the product type
    int first_minus_pos = short_filename.find_first_of("-");
    string radix=short_filename.substr(0,first_minus_pos);
    if (radix=="CAL_IIR_L1")
      product_type=L1;
    else if (radix=="CAL_IIR_L2_Track")
      product_type=L2_TRACK;
    else if (radix=="CAL_IIR_L2_Track")
      product_type=L2_SWATH;
    else
      product_type=UNDEFINED;
    // set the day/night mode
    char s_day_mode=short_filename[short_filename.size()-5];
    if ( s_day_mode == 'D' )
        day_mode=true;
    else if ( s_day_mode == 'N')
        day_mode=false;
    // set the acquisition date
    string s_date = short_filename.substr(short_filename.size()-25,21);
    date->set_date_str(s_date,"%Y-%m-%eT%H-%M-%S");
}
const int IIRFileData::get_track_size() const {
    return lat_lon_index_max[0]+1;
}
const int IIRFileData::get_swath_size( ) const {
    if (product_type==L1 || product_type==L2_SWATH)
      // 2D images
      return IIR_swath_size;
    return 1;
}
void IIRFileData::get_latitude_data( float * data ) {
    int start[] = {0,1};
    int edges[] = {get_track_size(),get_swath_size()};
    read_data(static_cast<void*>(data),latitude_sds_name.c_str(),start,NULL,edges);
}
void IIRFileData::get_longitude_data( float * data ) {
    int start[] = {0,1};
    int edges[] = {get_track_size(),get_swath_size()};
    read_data(static_cast<void*>(data),longitude_sds_name.c_str(),start,NULL,edges);
}
void IIRFileData::get_time_data( double * data ) {
    int start[] = {0,1};
    int edges[] = {get_track_size(),get_swath_size()};
    read_data(static_cast<void*>(data),time_sds_name.c_str(),start,NULL,edges);
}
const bool IIRFileData::is_day() const {
    return day_mode;
}
const int IIRFileData::get_nearest_point_index( const float & lat, const float & lon, const double & lidar_shot_time) {
    int nearest_index=-1;

    bool data_already_loaded=is_geolocation_data_loaded();
    if (!data_already_loaded)
        load_geolocation_data();

    int nb_pts = get_track_size();
    float* p_lat = lat_data;
    float* p_lon = lon_data;

    // search for the nearest point
    double delta;
    if (nb_pts!=0) {
        nearest_index=0;
        delta=sqrt( pow(*(p_lat+0)-lat,2) + pow(*(p_lon+0)-lon,2) );
        double current_delta;
        for (int i = 1 ; i<nb_pts ; ++i) {
            current_delta=sqrt(  pow(*(p_lat+i)-lat,2) + pow(*(p_lon+i)-lon,2) );
            if (current_delta<delta) {
                delta=current_delta;
                nearest_index=i;
            }
        }
    }
    if (!data_already_loaded)
        free_geolocation_data();
    return nearest_index;
}

const bool IIRFileData::get_index(const float & lat, const float & lon, const double & time, int * nearest_pix_idx, const float coloc_tolerance) {
    nearest_pix_idx[0]=-1; nearest_pix_idx[1]=-1; // default output

    bool data_already_loaded=is_geolocation_data_loaded();
    if (!data_already_loaded)
        load_geolocation_data();

    double delta_lat,delta_lon,delta,best_delta;
    delta_lat=delta_lon=delta=best_delta=9999.; // 9999 is bigger than the biggest delta that can be
    int _nearest_pix_idx[]={-1,-1};

    /*** use a generic search on sorted datas ***/
    // bounds of the colocation frame
    float lat_min=lat-coloc_tolerance; // lower latitude acceptable
    float lon_min=lon-coloc_tolerance; // lower longitude acceptable
    PixelType pix_min(&lat_min,&lon_min,vector<int>(0));
    float lat_max=lat+coloc_tolerance; // biggest latitude acceptable
    float lon_max=lon+coloc_tolerance; // biggest longitude acceptable
    PixelType pix_max(&lat_max,&lon_max,vector<int>(0));

    // set the range of colocated points candidates
    I_Pixel i_pix_min=lower_bound(v_pixel.begin(),v_pixel.end(),pix_min);
    I_Pixel i_pix_max=upper_bound(v_pixel.begin(),v_pixel.end(),pix_max);
//     cout<<"Pix Min ["<<i_pix_min->get_val()[0]<<","<<i_pix_min->get_val()[1]<<"]"<<endl;
//     cout<<"Pix Max ["<<i_pix_max->get_val()[0]<<","<<i_pix_max->get_val()[1]<<"]"<<endl;
//     cout<<"Distance "<<distance(i_pix_min,i_pix_max)<<endl;
    // cross the candidates to match the nearest one
    for (I_Pixel i_pix=i_pix_min;i_pix!=i_pix_max;++i_pix) {
      delta_lat=abs(i_pix->get_lat()-lat);
      delta_lon=abs(i_pix->get_lon()-lon);
      if (delta_lat<coloc_tolerance && delta_lon<coloc_tolerance) {
        // the pixel is in the colocalisation frame
//         cout<<"Possible point ("<<i_pix->get_lat()<<","<<i_pix->get_lon()<<") ["<<i_pix->get_val()[0]<<","<<i_pix->get_val()[1]<<"]"<<endl;
        delta=sqrt(pow(delta_lat,2)+pow(delta_lon,2));
        if (delta<best_delta) {
            best_delta=delta;
            _nearest_pix_idx[0]=i_pix->get_val()[0];
            _nearest_pix_idx[1]=i_pix->get_val()[1];
        }
      }
    }

    if (!data_already_loaded)
        free_geolocation_data();

    nearest_pix_idx[0]=_nearest_pix_idx[0];
    nearest_pix_idx[1]=_nearest_pix_idx[1];
    return (nearest_pix_idx[0]!=-1 && nearest_pix_idx[1]!=-1);
}
void IIRFileData::set_time_coverage() {
    int start[]={get_track_size()-1,0};
    int edges[]={1,1};
    int rank = 2;

    double time_coverage;
    // read the time of the last track measure
    read_data(&time_coverage,time_sds_name.c_str(),start,NULL,edges,rank);
    // compute the duration using the time of the first track measure
    time_coverage-=get_date().get_TAI93_time();
}
const bool IIRFileData::contains_data_at( const float & lat, const float & lon, const double & time ) const
{
    double start_time = get_date().get_TAI93_time();
    double end_time = start_time+time_coverage;
    bool in_time = (time>=start_time && time<=end_time);
    bool found = in_time;
    return found;
}
const bool IIRFileData::is_geolocation_data_loaded() const {
    return (lat_data!=NULL && lon_data!=NULL && time_data!=NULL);
}
void IIRFileData::free_geolocation_data() {
    if (lat_data!=NULL) {
        delete[] lat_data;
        lat_data=NULL;
    }
    if (lon_data!=NULL) {
        delete[] lon_data;
        lon_data=NULL;
    }
    if (time_data!=NULL) {
        delete[] time_data;
        time_data=NULL;
    }
}
void IIRFileData::load_geolocation_data() {
    if (!is_geolocation_data_loaded()) {
        // hdf_file is needed
        bool hdf_file_already_loaded = is_hdf_file_loaded();
        if (!hdf_file_already_loaded)
            load_hdf_file();

        // read lat, lon data
        int start[] = {0,0};
        int edges[] = {get_track_size(),get_swath_size()};
        if (lat_data==NULL)
            lat_data=static_cast<float*>(read_data(static_cast<void*>(lat_data),latitude_sds_name.c_str(),start,NULL,edges));
        if (lon_data==NULL)
            lon_data=static_cast<float*>(read_data(static_cast<void*>(lon_data),longitude_sds_name.c_str(),start,NULL,edges));
        if (time_data==NULL)
            time_data=static_cast<double*>(read_data(static_cast<void*>(time_data),time_sds_name.c_str(),start,NULL,edges));

        // build the sorted pixels list
        load_v_pixel();

        if (!hdf_file_already_loaded)
            free_hdf_file();
    }
}

void IIRFileData::load_v_pixel()
{
    vector<int> sds_index(2); // indexes [y,x] of a pixel in a MODIS SDS
    int y_max=lat_lon_index_max[0]+1,x_max=lat_lon_index_max[1]+1;
    int buffer_idx=0; // linear index to access 2D lat/lon data arrays
    for (int y_idx=0;y_idx<y_max;++y_idx) {
        for (int x_idx=0;x_idx<x_max;++x_idx) {
        buffer_idx=y_idx*x_max+x_idx;
        sds_index[0]=y_idx;
        sds_index[1]=x_idx;
        v_pixel.push_back(PixelType(lat_data+buffer_idx,lon_data+buffer_idx,sds_index));
        }
    }
    sort(v_pixel.begin(),v_pixel.end());
}

const bool IIRFileData::get_index(const float &lat, const float& lon, int &nearest_pix_idx, const float coloc_tolerance)
{
    int nearest_pix_idx2[]={-1,-1};

    get_index(lat, lon, -1.0, nearest_pix_idx2, coloc_tolerance);
    if(nearest_pix_idx2[0]==-1 || nearest_pix_idx2[1]==-1) return false;

    int x_max=lat_lon_index_max[1]+1;
    nearest_pix_idx=nearest_pix_idx2[0]*x_max+nearest_pix_idx2[1];

    return true;
}

const float IIRFileData::get_nearest_point_distance( const float & lat, const float & lon, const float coloc_tolerance ) {
    bool data_already_loaded=is_geolocation_data_loaded();
    if (!data_already_loaded)
        load_geolocation_data();

    int nearest_point_idx;
    double nearest_point_distance=-1.;
    if (get_index(lat,lon,nearest_point_idx,coloc_tolerance))
        // a point in the coincidence frame has been found
        nearest_point_distance=sqrt(pow(lat_data[nearest_point_idx]-lat,2) + pow(lon_data[nearest_point_idx]-lon,2));

    if (!data_already_loaded)
        free_geolocation_data();

    return nearest_point_distance;
}

const bool IIRFileData::contain_location(const float &lat,const float &lon, const double &tolerance){
  int nearest_idx=-1;
  return get_index(lat,lon,nearest_idx,colocation_tolerance);
}
const bool IIRFileData::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);
}

#if IIRFILEDATA_TEST != 0
int main(int argc, char *argv[]) {
    string test_filename("/icartoon4/CALIPSO/CALIPSO_MODIS-LITE/CAL_IIR_L2_Swath-BETA-V01.2003-12-31T22-41-30ZN.hdf");
    IIRFileData iirfd(test_filename);
    return 0;
}
#endif