/***************************************************************************
*   Copyright (C) 2005 by Nicolas PASCAL                                  *
*   pascal@icare-pc15                                                     *
*                                                                         *
*   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 "modisfiledata.h"

const int MODISFileData::NB_MAX_X_10KM  = 135;
const int MODISFileData::NB_MAX_Y_10KM  = 204;
const int MODISFileData::NB_MAX_X_5KM   = 271;
const int MODISFileData::NB_MAX_Y_5KM   = 410;
const int MODISFileData::NB_MAX_X_1KM   = 1354;
const int MODISFileData::NB_MAX_Y_1KM   = 2054;
const int MODISFileData::NB_MAX_X_250M  = 5416;
const int MODISFileData::NB_MAX_Y_250M  = 8216;

const string MODISFileData::latitude_sds_name="Latitude";
const string MODISFileData::longitude_sds_name="Longitude";

MODISFileData::MODISFileData(const string &name /*= ""*/, const string &mode/*= "r"*/) :
    FileData(name,mode),SatelliteFileData(name,mode),HDFFileData(name,mode) {
    init();
}
void MODISFileData::init() {
    // extract the acquisition date from the file name
    string short_filename = get_tail(get_name());
    lat_data=NULL;
    lon_data=NULL;
    lat_lon_index_max.resize(2); // lat,lon sds have a fixed size for MODIS
    gring.reserve(2); // gring is defined by 2 polygons maximum
    lat_gran_corner.resize(4);
    lon_gran_corner.resize(4);
    fill_n ( lat_gran_corner.begin(), 4, -9999. );
    fill_n ( lon_gran_corner.begin(), 4, -9999. );

    parse_filename(short_filename);

    init_colocation_tolerance ();

    load_hdf_file();
    load_hdf_metadata();

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

    free_hdf_metadata();
    free_hdf_file();
}

void MODISFileData::print_gring() const {
    cout<<"----- GRing -----"<<endl;
    for (uint i=0;i<gring.size();++i) {
        cout<<"-> Polygon "<<i<<"<- size "<<gring[i].size()<<endl;
        for (uint j=0;j<gring[i].size();++j)
            cout<<gring[i][j][0]<<","<<gring[i][j][1]<<endl;
        cout<<"-> End Polygon "<<i<<"<-"<<endl;
    }
    cout<<"----- End GRing -----"<<endl;
}
void MODISFileData::set_gring() {
    // load the GRing values
    vector <double> gring_lat; gring_lat.reserve(4); /* the GRing latitude points */
    vector <double> gring_lon; gring_lon.reserve(4); /* the GRing longitude points */
    metadata->get_value((void*)&gring_lat,get_gring_latitude_name().c_str());
    metadata->get_value((void*)&gring_lon,get_gring_longitude_name().c_str());

    // build the GRing polygons
    typedef vector<double> Point;
    typedef vector<Point> Line;
    typedef vector<Point> Polygon;
    // initialize the GRing
    gring.clear();
    //gring.push_back(Polygon());
    gring.push_back(Polygon()); // at least one polygon set
    GRing::iterator current_poly=gring.begin(); // the current polygon that is in the way to be build
    Point p1(2); // first point
    Point p2(2); // point following p1
    Point p_inter(2); // eventual intersection point
    Line l_gring; l_gring.reserve(2); // line defined by the 2 following points of the GRing
    Point p1_lon_axis(2); p1_lon_axis[0]=-10.; p1_lon_axis[1]=180.;
    Point p2_lon_axis(2); p2_lon_axis[0]= 10.;  p2_lon_axis[1]=180.;
    // line defining the longitude axis (-180 or 180 depending of the case)
    Line l_lon_axis; l_lon_axis.reserve(2); l_lon_axis.push_back(p1_lon_axis); l_lon_axis.push_back(p2_lon_axis);
    uint last_point=gring_lat.size()-1;
    for (uint i=0;i<=last_point;++i) { // cross the GRing values
        p1[0]=gring_lat[i]; p1[1]=gring_lon[i];
        if (i==last_point) {
            // point that follows the last point is the first one to close the GRing
            p2[0]=gring_lat[0]; p2[1]=gring_lon[0];
        } else {
            // next point
            p2[0]=gring_lat[i+1]; p2[1]=gring_lon[i+1];
        }

        //cout<<"Segment ["<<p1[0]<<","<<p1[1]<<"]->["<<p2[0]<<","<<p2[1]<<"]"<<endl;

        // check if the segment between p1 and p2 crosses the change date longitude
        l_gring.clear();
        // set the gring segment
        l_gring.push_back(p1);
        l_gring.push_back(p2);
        if (abs(p2[1]-p1[1])>=180.) { // p1 and p2 crosses the change date longitude axis
            // split the segment in 2 : one to the -180 longitude line, the other to the +180 one
            double crossed_lon_axis; // wether the first crossed axis is the +180. or -180. longitude
            crossed_lon_axis=(p1[1]<0.?-180.:180.);
            // set a longitude axis segment
            l_lon_axis[0][1]=crossed_lon_axis;
            l_lon_axis[1][1]=crossed_lon_axis;
            // search intersection with the first crossed axis
            MyTools::get_intersection<double>(l_gring,l_lon_axis,p_inter);
            current_poly->push_back(p1);
            current_poly->push_back(p_inter);
            // change current polygon
            if (current_poly==gring.begin()) {
                // if unexisting 2nd polygon -> create it
                  if (gring.size()==1)
                     gring.resize(2);
                    //gring.push_back(Polygon());
                ++current_poly;
            } else // last polygon must be close
                --current_poly;
            // build intersection with the complementary longitude axis (if first crossed axis = 180, the complementary is -180...)
            p_inter[1]=-crossed_lon_axis;
            current_poly->push_back(p_inter);
        } else { // general case : p1 and p2 do not cross the change date longitude
            current_poly->push_back(p1);
        }
    }
}

void MODISFileData::set_lat_lon_index_max() {
    assert(is_hdf_file_loaded());
    // general case : lat, lon data will be the size of the geolocation datasets as given in the product
    lat_lon_index_max[0] = hdf_file->get_sds(latitude_sds_name.c_str()).get_dimension(0)-1;
    lat_lon_index_max[1] = hdf_file->get_sds(latitude_sds_name.c_str()).get_dimension(1)-1;
}

MODISFileData::~MODISFileData() {
    free_geolocation_data();
    free_hdf_metadata();
    free_hdf_file();
}
void MODISFileData::set_time_coverage() {
    if (  product == "MYD02QKM" || product == "MYD02HKM" || product == "MYD021KM" || product == "MYD03" ||
         product == "MYD04_L2" || product == "MYD05_L2" || product == "MYD06_L2" )
        time_coverage=5.*60.; // those files have a duration of 5 min
    else
        cout<<"MODISFileData::set_time_coverage : the time coverage of "<<product<<" have to be added"<<endl;
}
void MODISFileData::parse_filename( const string& short_filename ) {
    check_filename(short_filename); // throw a parse_MODIS_filename_error if the filename doesn't matches the pattern
    string modis_filename(short_filename);

    // the 3 first characters define the platform
    string platform_code=modis_filename.substr(0,3);
    if (platform_code=="MOD")
        platform=TERRA;
    else if (platform_code=="MYD")
        platform=AQUA;
    else {
        cerr<<"MODISFileData::parse_filename : Unknown MODIS platform "<<platform_code<<" in "<<__FILE__<<" at "<<__LINE__<<endl;
        platform=UNDEFINED_PLATFORM_TYPE;
    }

    // extract the product name
    int first_dot_pos = modis_filename.find_first_of(".",0);
    product=modis_filename.substr(0,first_dot_pos);
    if      (product=="MYD02QKM")
        product_id =   MYD02QKM;
    else if (product=="MYD02HKM")
        product_id =   MYD02HKM;
    else if (product=="MYD021KM")
        product_id =   MYD021KM;
    else if (product=="MYD03")
        product_id =   MYD03;
    else if (product=="MYD02SSH")
        product_id =   MYD02SSH;
    else if (product=="MYD04_L2")
        product_id =   MYD04_L2;
    else if (product=="MYD05_L2")
        product_id =   MYD05_L2;
    else if (product=="MYD06_L2")
        product_id =   MYD06_L2;
    else if (product=="MOD06_L2")
        product_id =   MOD06_L2;
    else // unsupported product
    {
        product_id = UNDEFINED_PRODUCT_ID;
        string msg ( "Unsupported MODIS product " + product ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }

    // set the geolocation dataset resolution, in km
    switch ( product_id ) {
        case ( MYD02QKM ) :
        case ( MYD03    ) :
            resolution=1.;
            break;
        case ( MYD02HKM ) :
            resolution=.5;
            break;
        case ( MYD02SSH ) :
        case ( MYD05_L2 ) :
        case ( MYD06_L2 ) :
        case ( MYD021KM ) :
            resolution=5.;
            break;
        case ( MYD04_L2 ) :
            resolution=10.;
            break;
        default: // should never arrive here
            string msg ( "Unsupported MODIS product " + product ) ;
            g_exception e( __FILE__ , __LINE__ , msg );
            throw e;
            break;
    }

    // date info is contained between first and third dot position
    int second_dot_pos = modis_filename.find_first_of(".",first_dot_pos+1);
    int third_dot_pos = modis_filename.find_first_of(".",second_dot_pos+1);
    string s_date = modis_filename.substr(first_dot_pos+1,third_dot_pos-first_dot_pos-1);
    s_date+= "00"; // add default seconds
    date->set_date_str(s_date,"A%Y%j.%H%M%S");
}
bool MODISFileData::check_filename( const string& short_filename ) const {
    string modis_filename(short_filename);

    string _type = modis_filename.substr(0,3);
    string _ext = modis_filename.substr(modis_filename.size()-4,4);
    size_t first_dot_pos = modis_filename.find_first_of(".",0);
    string prod_id = modis_filename.substr ( 3, first_dot_pos - 3 );

    if  ( ( _type == "MYD" || _type == "MOD" ) &&
          ( prod_id == "02QKM" || prod_id == "02HKM" || prod_id == "021KM" ||
            prod_id == "03" || prod_id == "02SSH" ||
            prod_id == "04_L2" || prod_id == "05_L2" || prod_id == "06_L2" ) &&
          ( _ext == ".hdf" ) ) {
        return true;
    } else {
        invalid_filename e(__FILE__,__LINE__,modis_filename,"MODIS");
        throw e;
        return false; // inutile
    }
}
void MODISFileData::set_lat_lon_min_max() {
    bool already_loaded_metadata = is_hdf_metadata_loaded();
    if (!already_loaded_metadata)
        load_hdf_metadata();

    // read min,max in the metadata
    metadata->get_value((void*)(&lat_max),"NORTHBOUNDINGCOORDINATE");
    metadata->get_value((void*)(&lat_min),"SOUTHBOUNDINGCOORDINATE");
    metadata->get_value((void*)(&lon_max),"EASTBOUNDINGCOORDINATE");
    metadata->get_value((void*)(&lon_min),"WESTBOUNDINGCOORDINATE");

    if (!already_loaded_metadata)
        free_hdf_metadata();
}

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

    int nearest_point_idx[2];
    double nearest_point_distance=-1.;

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


    if (!data_already_loaded)
        free_geolocation_data();

    return nearest_point_distance;
}
const bool MODISFileData::get_index( const float &lat, const float & lon, const double time, int * nearest_pix_idx, float coloc_tolerance) {
    fill_n(nearest_pix_idx, 2, -1); // default output
//     if (!contain_time(time)) // if time = -1, contain_time will be true
//         return false;

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

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

    // set default colocation tolerance if not given
    if (coloc_tolerance<0.)
        coloc_tolerance = colocation_tolerance;

    // conversion of tolerance in degree -> km
    double tolerance_m = 111.11e3 * coloc_tolerance;

    /*** 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-999; // 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+999; // 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 << "--- searched pixel (" << lat << "," << lon << ") ---" << endl;
//     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) {
      if ( abs(i_pix->get_lat()-lat) < coloc_tolerance ) {
        // the pixel is in the colocalisation frame
//         delta = sqrt(pow(delta_lat,2)+pow(delta_lon,2));
        delta = Geocentric::get_haversine_dist (lat, lon, i_pix->get_lat(), i_pix->get_lon());
//         cout<<"   -> Possible point ("<<i_pix->get_lat()<<","<<i_pix->get_lon()<<") ["<<i_pix->get_val()[0]<<","<<i_pix->get_val()[1]<<"] dist=" << delta <<endl;

        if ( delta < tolerance_m && 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();
    copy (_nearest_pix_idx, _nearest_pix_idx + 2, nearest_pix_idx);

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

const bool MODISFileData::get_index(const float &lat, const float& lon, int& nearest_pix_idx, float coloc_tolerance)
{
    int nearest_idx[2];
    nearest_idx[0]=-1;
    nearest_idx[1]=-1;

    get_index(lat, lon, -1.0 , nearest_idx, coloc_tolerance);

    if(nearest_idx[0]==-1 && nearest_idx[1]==-1) return false;

    nearest_pix_idx = nearest_idx[0] + nearest_idx[1] * lat_lon_index_max[0];

    return true;
}

const string MODISFileData::get_radix() {
    string year=MyTools::to_string(date->get_year());
    string month = Date::get_month_string(date->get_month());
    string day = Date::get_month_string(date->get_day());
    string s_date = year+month+day;
    return product+".A"+s_date;
}
const bool MODISFileData::contain_data( const float & lat, const float & lon, const double & time, const double &tolerance ) {
    return contain_location(lat,lon,tolerance) && contain_time(time);
}
const bool MODISFileData::contain_location( const float & lat, const float & lon, const double &tolerance  ) {
    int idx[]={-1,-1};
    get_index(lat, lon, -1., idx, tolerance);
    return (idx[0]!=-1 && idx[1]!=-1);
    //return is_inside_gring(lat,lon);
}
const bool MODISFileData::is_inside_gring( const float & lat, const float & lon ) const {
    //assert(gring_lat.size()==gring_lon.size());
    double my_lat(lat);
    double my_lon(lon);
    double angle=0;
    double p1[2];
    double p2[2];
    uint end=0;
    for (uint polygon=0;polygon<gring.size();polygon++) { // for each polygon that compose the GRing (if many)
        end=gring[polygon].size();
        for (uint point=0;point<end;point++) { // for each point of the polygon
            p1[1] = gring[polygon][point][0] - my_lat;
            p1[0] = gring[polygon][point][1] - my_lon;
            p2[1] = gring[polygon][(point+1)%end][0] - my_lat;
            p2[0] = gring[polygon][(point+1)%end][1] - my_lon;
            angle += MyTools::angle_2D(p1,p2);
        }
    }
    if (abs(angle) < M_PI)
        return false;
    else
        return true;
}
void MODISFileData::get_5_to_1km_area( const int * five_km_idx, int * one_km_min_idx, int * one_km_max_idx ) {
    one_km_min_idx[0]=one_km_min_idx[1]=one_km_max_idx[0]=one_km_max_idx[1]=-1;
    if (five_km_idx[0]<0 || five_km_idx[0]>lat_lon_index_max[0]) {
        bad_index e(five_km_idx[0],lat_lon_index_max[0]);
        throw e;
    } else if (five_km_idx[1]<0 || five_km_idx[1]>lat_lon_index_max[1]) {
        bad_index e(five_km_idx[1],lat_lon_index_max[1]);
        throw e;
    } else {
        one_km_min_idx[0]=five_km_idx[0]*5;
        one_km_max_idx[0]=min(one_km_min_idx[0]+4,static_cast<int>(lat_lon_index_max[0]));
        one_km_min_idx[1]=five_km_idx[1]*5;
        one_km_max_idx[1]=min(one_km_min_idx[1]+4,static_cast<int>(lat_lon_index_max[1]));
    }
}
const int MODISFileData::convert_5_to_1km_index( const int &five_km_idx, const bool is_x_axis ) const{
    int ret=-1;
    if (!is_x_axis && (five_km_idx<0 || five_km_idx>lat_lon_index_max[0]) ) {
        bad_index e(five_km_idx,lat_lon_index_max[0]);
        throw e;
    } else if (is_x_axis && (five_km_idx<0 || five_km_idx>lat_lon_index_max[1]) ) {
        bad_index e(five_km_idx,lat_lon_index_max[1]);
        throw e;
    } else {
        ret=2+5*five_km_idx;
    }
    return ret;
}

void MODISFileData::convert_1km_to_250m_index ( const int  i_pix_1km  [],
                                                      float i_pix_250m []
                                              ) const
{
    if ( i_pix_1km [0] < 0 || i_pix_1km [0] > lat_lon_index_max [0]  )
    {
        bad_index e( i_pix_1km [0], lat_lon_index_max[0] );
        throw e;
    } else if ( i_pix_1km [1] < 0 || i_pix_1km [1] > lat_lon_index_max [1]  )
    {
        bad_index e( i_pix_1km [1], lat_lon_index_max[1] );
        throw e;
    }

    // along track
    i_pix_250m [0] = 1.5 + 4 * i_pix_1km [0];
    // along swath
    i_pix_250m [1] =       4 * i_pix_1km [1];
}

const bool MODISFileData::get_1km_to_250m_pixel( const float32 lat,
                                                 const float32 lon,
                                                 const int i_pix_1km[],
                                                       int i_pix_250m[]
                                                )
{
    assert( i_pix_250m != NULL );

    fill_n ( i_pix_250m, 2, -1 );
    if ( i_pix_1km[0] == -1 || i_pix_1km[1] == -1 )
        return false;

    // seach the indexes bounds of the 250m pixels that are situed in the 1km footprint
    int itk_250m_min, itk_250m_max, isc_250m_min, isc_250m_max;
    get_bounds_1km_to_250m ( i_pix_1km[0], i_pix_1km[1], itk_250m_min, itk_250m_max, isc_250m_min, isc_250m_max );
//     printf ( "pix_1km [%ld,%ld]\t(%f,%f)\n", i_pix_1km[0], i_pix_1km[1], lat, lon );
//     printf ( "itk_250m = [%ld,%ld]\tisc_250m = [%ld,%ld]\n", itk_250m_min, itk_250m_max, isc_250m_min, isc_250m_max );

    // in these 250m pixels, search the nearest to lat, lon
    double best_dist = 99999., dist = 99999.;
    float32 _lat_250m, _lon_250m; // lat, lon of the 250m pixels
    int itk_250m, isc_250m;
    for ( itk_250m = itk_250m_min ; itk_250m <= itk_250m_max ; ++itk_250m ) {
        for ( isc_250m = isc_250m_min ; isc_250m <= isc_250m_max ; ++isc_250m ) {
            get_250m_pix_pos( itk_250m, isc_250m, _lat_250m, _lon_250m );
            dist = Geocentric::get_haversine_dist( lat, lon, _lat_250m, _lon_250m );
// printf ( "pix_250m [%ld,%ld]\t(%f,%f)\tdist=%f\n", itk_250m, isc_250m, _lat_250m, _lon_250m, dist );

            if ( dist < best_dist ) {
                // a nearest point has been found -> store it
                best_dist = dist;
                i_pix_250m [0] = itk_250m;
                i_pix_250m [1] = isc_250m;
//                 lat_250m = _lat_250m;
//                 lon_250m = _lon_250m;
            }
        }
    }

//     printf ( "i_pix_250m = [%ld,%ld]\n", i_pix_250m[0], i_pix_250m[1] );

    return true;
}

const bool MODISFileData::get_5km_to_1km_pixel( const float32 lat,
                                                 const float32 lon,
                                                 const int i_pix_5km[],
                                                       int i_pix_1km[]
                                                )
{
    assert( i_pix_1km != NULL );

    fill_n ( i_pix_1km, 2, -1 );
    if ( i_pix_5km[0] == -1 || i_pix_5km[1] == -1 )
        return false;

    // seach the indexes bounds of the 1km pixels that are situed in the 5km footprint
    int itk_1km_min, itk_1km_max, isc_1km_min, isc_1km_max;
    get_bounds_5km_to_1km ( i_pix_5km[0], i_pix_5km[1], itk_1km_min, itk_1km_max, isc_1km_min, isc_1km_max );
//     printf ( "pix_5km [%ld,%ld]\t(%f,%f)\n", i_pix_5km[0], i_pix_5km[1], lat, lon );
//     printf ( "itk_1km = [%ld,%ld]\tisc_1km = [%ld,%ld]\n", itk_1km_min, itk_1km_max, isc_1km_min, isc_1km_max );

    // in these 1km pixels, search the nearest to lat, lon
    double best_dist = 99999., dist = 99999.;
    float32 _lat_1km, _lon_1km; // lat, lon of the 1km pixels
    int itk_1km, isc_1km;
    for ( itk_1km = itk_1km_min ; itk_1km <= itk_1km_max ; ++itk_1km ) {
        for ( isc_1km = isc_1km_min ; isc_1km <= isc_1km_max ; ++isc_1km ) {
            get_1km_pix_pos( itk_1km, isc_1km, _lat_1km, _lon_1km );
            dist = Geocentric::get_haversine_dist( lat, lon, _lat_1km, _lon_1km );
// printf ( "pix_1km [%ld,%ld]\t(%f,%f)\tdist=%f\n", itk_1km, isc_1km, _lat_1km, _lon_1km, dist );
            if ( dist < best_dist ) {
                // a nearest point has been found -> store it
                best_dist = dist;
                i_pix_1km [0] = itk_1km;
                i_pix_1km [1] = isc_1km;
//                 lat_1km = _lat_1km;
//                 lon_1km = _lon_1km;
            }
        }
    }
//     printf ( "i_pix_1km = [%ld,%ld]\n", i_pix_1km[0], i_pix_1km[1] );
    return true;
}

const bool MODISFileData::get_5_to_1km_index( const float & lat, const float & lon, int * idx_1km, const int * idx_5km ) {
    assert(idx_1km!=NULL);
    assert(idx_5km[0]>=0);
    assert(idx_5km[0]<=lat_lon_index_max[0]);
    assert(idx_5km[1]>=0);
    assert(idx_5km[1]<=lat_lon_index_max[1]);
//     if (!is_inside_gring(lat,lon))
//         return false;

    // limits
    int size_y=lat_lon_index_max[0]+1;
    int size_x=lat_lon_index_max[1]+1;
    int npix_1km[]={5*size_y,5*size_x-1};

    // 4 nearest pixel indexes at 5km resol
    int coarsepix1_idx_5km[]={-1,-1};
    int coarsepix2_idx_5km[]={-1,-1};
    int coarsepix3_idx_5km[]={-1,-1};
    int coarsepix4_idx_5km[]={-1,-1};
    assert(is_geolocation_data_loaded());

    /* Use a bilinear interpolation to find the indexes of the 1km pixel that contains the point (lat,lon)
       To do that :
       1- Find 4 coarse pixels that are neighbours of (lat,lon). The 1st one contains this measure
       2- Compute the indexes (at 1km) of the 1km pixel that is at the center of the coarse one
       3- Read the (lat,lon) of the 1km pixel that is in the center of the coarse one
       4- Use a bilinear interpolation to map the desired (lat,lon) measure to its 1km indexes (using the 4 points as references)
    */

    //--- Set the 5km indexes of the 4 coarse pixels that are neighbours of (lat,lon). The 1st one contains this measure
    // -> 1st pixel -> only a copy of idx_5km
    copy(idx_5km,idx_5km+2,coarsepix1_idx_5km);
    // read its coordinates
    int linear_pix_idx_5km=coarsepix1_idx_5km[0]*size_x+coarsepix1_idx_5km[1];
    float coarsepix1_coord[]={lat_data[linear_pix_idx_5km],lon_data[linear_pix_idx_5km]};
    // along Y axes
    if (idx_5km[0]==lat_lon_index_max[0]) {
      coarsepix2_idx_5km[0]=idx_5km[0];
      coarsepix3_idx_5km[0]=idx_5km[0]-1;
      coarsepix4_idx_5km[0]=idx_5km[0]-1;
    } else if (idx_5km[0]==0) {
      coarsepix2_idx_5km[0]=idx_5km[0];
      coarsepix3_idx_5km[0]=idx_5km[0]+1;
      coarsepix4_idx_5km[0]=idx_5km[0]+1;
    } else { // general case
      if (lat>=coarsepix1_coord[0]) { // point is "above" the central pixel
        coarsepix2_idx_5km[0]=idx_5km[0];
        coarsepix3_idx_5km[0]=idx_5km[0]-1;
        coarsepix4_idx_5km[0]=idx_5km[0]-1;
      } else { // point is "below" the central pixel
        coarsepix2_idx_5km[0]=idx_5km[0];
        coarsepix3_idx_5km[0]=idx_5km[0]+1;
        coarsepix4_idx_5km[0]=idx_5km[0]+1;
      }
    }
    // along X axis
    if (idx_5km[1]==lat_lon_index_max[1]) {
      coarsepix2_idx_5km[1]=idx_5km[1]-1;
      coarsepix3_idx_5km[1]=idx_5km[1];
      coarsepix4_idx_5km[1]=idx_5km[1]-1;
    } else if (idx_5km[1]==0) {
      coarsepix2_idx_5km[1]=idx_5km[1]+1;
      coarsepix3_idx_5km[1]=idx_5km[1];
      coarsepix4_idx_5km[1]=idx_5km[1]+1;
    } else { // general case
      if (lon>=coarsepix1_coord[1]) { // point is "on the right" of the central pixel
        coarsepix2_idx_5km[1]=idx_5km[1]+1;
        coarsepix3_idx_5km[1]=idx_5km[1];
        coarsepix4_idx_5km[1]=idx_5km[1]+1;
      } else { // point is "on the left" of the central pixel
        coarsepix2_idx_5km[1]=idx_5km[1]-1;
        coarsepix3_idx_5km[1]=idx_5km[1];
        coarsepix4_idx_5km[1]=idx_5km[1]-1;
      }
    }

    //--- Compute the indexes (at 1km) of the 1km pixel that is at the center of the coarse one
    int coarsepix1_idx_1km[]={coarsepix1_idx_5km[0]*5+2,coarsepix1_idx_5km[1]*5+2};
    int coarsepix2_idx_1km[]={coarsepix2_idx_5km[0]*5+2,coarsepix2_idx_5km[1]*5+2};
    int coarsepix3_idx_1km[]={coarsepix3_idx_5km[0]*5+2,coarsepix3_idx_5km[1]*5+2};
//     int coarsepix4_idx_1km[]={coarsepix4_idx_5km[0]*5+2,coarsepix4_idx_5km[1]*5+2};

    //--- Read the coordinates of the 1km pixel that is at the centre of each coarse ones
    // -> 1st pixel already read
/*    int linear_pix_idx_5km=coarsepix1_idx_5km[0]*NB_X_5KM+coarsepix1_idx_5km[1];
    float coarsepix1_coord[]={lat_data[linear_pix_idx_5km],lon_data[linear_pix_idx_5km]};*/
    // -> 2nd pixel
    linear_pix_idx_5km=coarsepix2_idx_5km[0]*size_x+coarsepix2_idx_5km[1];
    float coarsepix2_coord[]={lat_data[linear_pix_idx_5km],lon_data[linear_pix_idx_5km]};
    // -> 3rd pixel
    linear_pix_idx_5km=coarsepix3_idx_5km[0]*size_x+coarsepix3_idx_5km[1];
    float coarsepix3_coord[]={lat_data[linear_pix_idx_5km],lon_data[linear_pix_idx_5km]};
    // -> 4th pixel
    linear_pix_idx_5km=coarsepix4_idx_5km[0]*size_x+coarsepix4_idx_5km[1];
//     float coarsepix4_coord[]={lat_data[linear_pix_idx_5km],lon_data[linear_pix_idx_5km]};

    //--- Apply the interpolation
    idx_1km[0]=(int)(linear_interpolation(lat,coarsepix1_coord[0],coarsepix3_coord[0],coarsepix1_idx_1km[0],coarsepix3_idx_1km[0]));
    idx_1km[1]=(int)(linear_interpolation(lon,coarsepix1_coord[1],coarsepix2_coord[1],coarsepix1_idx_1km[1],coarsepix2_idx_1km[1]));

//     cout<<"P1\t("<<coarsepix1_coord[0]<<","<<coarsepix1_coord[1]<<")\t5km ["<<coarsepix1_idx_5km[0]<<","<<coarsepix1_idx_5km[1]<<"]\t1km ["<<coarsepix1_idx_1km[0]<<","<<coarsepix1_idx_1km[1]<<"]"<<endl;
//     cout<<"P2\t("<<coarsepix2_coord[0]<<","<<coarsepix2_coord[1]<<")\t5km ["<<coarsepix2_idx_5km[0]<<","<<coarsepix2_idx_5km[1]<<"]\t1km ["<<coarsepix2_idx_1km[0]<<","<<coarsepix2_idx_1km[1]<<"]"<<endl;
//     cout<<"P3\t("<<coarsepix3_coord[0]<<","<<coarsepix3_coord[1]<<")\t5km ["<<coarsepix3_idx_5km[0]<<","<<coarsepix3_idx_5km[1]<<"]\t1km ["<<coarsepix3_idx_1km[0]<<","<<coarsepix3_idx_1km[1]<<"]"<<endl;
//     cout<<"P4\t("<<coarsepix4_coord[0]<<","<<coarsepix4_coord[1]<<")\t5km ["<<coarsepix4_idx_5km[0]<<","<<coarsepix4_idx_5km[1]<<"]\t1km ["<<coarsepix4_idx_1km[0]<<","<<coarsepix4_idx_1km[1]<<"]"<<endl;
//
//     cout<<"Pixel ("<<lat<<","<<lon<<") ["<<idx_1km[0]<<","<<idx_1km[1]<<"]"<<endl;

    // check the validity of the interpolated indexes
    if (idx_1km[0]>(npix_1km[0]-1)) idx_1km[0]=npix_1km[0]-1;
    if (idx_1km[1]>(npix_1km[1]-1)) idx_1km[1]=npix_1km[1]-1;
    if ( idx_1km[0]<0 || idx_1km[0]>(npix_1km[0]-1) || idx_1km[1]<0 || idx_1km[1]>(npix_1km[1]-1) ) {
      idx_1km[0]=-1,idx_1km[1]=-1;
      return false;
    }

    return true;
}
void MODISFileData::free_geolocation_data() {
    delete[] lat_data;
    lat_data=NULL;
    delete[] lon_data;
    lon_data=NULL;
    v_pixel.clear();
}

void MODISFileData::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();

        // load the data
        if (lat_data==NULL)
            lat_data=static_cast<float32*>(read_data((VOIDP)lat_data,latitude_sds_name.c_str()));
        if (lon_data==NULL)
            lon_data=static_cast<float32*>(read_data((VOIDP)lon_data,longitude_sds_name.c_str()));

        // read the geolocation of the granule corners
        int i_trck_max = lat_lon_index_max[0]; // i track max
        int i_scan_max = lat_lon_index_max[1]; // i scan max
        int sz_scan    = i_scan_max + 1;
        lat_gran_corner [0] = lat_data [0];
        lat_gran_corner [1] = lat_data [ i_scan_max ];
        lat_gran_corner [2] = lat_data [ i_trck_max * sz_scan + i_scan_max ];
        lat_gran_corner [3] = lat_data [ i_trck_max * sz_scan ];
        lon_gran_corner [0] = lon_data [0];
        lon_gran_corner [1] = lon_data [ i_scan_max ];
        lon_gran_corner [2] = lon_data [ i_trck_max * sz_scan + i_scan_max ];
        lon_gran_corner [3] = lon_data [ i_trck_max * sz_scan ];

    //         printf ( "lat_lon_index_max  : %d %d\n", lat_lon_index_max[0], lat_lon_index_max[1] );
    //         printf ( "Latitude of the granule corners  : %f %f %f %f\n", lat_gran_corner[0], lat_gran_corner[1], lat_gran_corner[2], lat_gran_corner[3] );
    //         printf ( "Longitude of the granule corners : %f %f %f %f\n", lon_gran_corner[0], lon_gran_corner[1], lon_gran_corner[2], lon_gran_corner[3] );

        if (!hdf_file_already_loaded)
            free_hdf_file();

        // build the sorted pixels list
        load_v_pixel();
    }
}

void MODISFileData::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;
    size_t 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));
// cout << "pix [" << y_idx << ", " << x_idx << "] " << v_pixel.back() << endl;
        }
    }
    sort(v_pixel.begin(), v_pixel.end());
}

const bool MODISFileData::is_geolocation_data_loaded() const {
    return (lat_data!=NULL && lon_data!=NULL);
}

void MODISFileData::get_sds_calibration( const string & sds_name, double & scale, double & offset, const string & scale_name, const string & offset_name, const int channel_nb )
{
    assert(is_hdf_file_loaded());
    // set default value index
    int ival=max(channel_nb,0);
    // read the scaling
    scale=0.;
    offset=0.;
    // Be careful of the scaling data type
    // - MYD06 -> double
    // - MYD021KM, MYD02QKM -> float
    if ( product == "MYD02QKM" || product == "MYD021KM" ) {
        float my_offset, my_scale;
        hdf_file->get_sds_attr ( sds_name.c_str(), scale_name.c_str(),  & my_scale,  ival );
        hdf_file->get_sds_attr ( sds_name.c_str(), offset_name.c_str(), & my_offset, ival);
        scale  = static_cast<float> ( my_scale );
        offset = static_cast<float> ( my_offset );
    } else {
        hdf_file->get_sds_attr ( sds_name.c_str(), scale_name.c_str(),  & scale,  ival );
        hdf_file->get_sds_attr ( sds_name.c_str(), offset_name.c_str(), & offset, ival);
    }

    //cout<<scale<<" "<<offset<<endl;
//     scale=my_scale;
//     offset=my_offset;
}
void MODISFileData::init_colocation_tolerance()
{
    if ( product_id == MYD02QKM )
        colocation_tolerance=0.01;
    else
        // set the default colocation tolerance : pixel size of 1km * resolution of the product
        colocation_tolerance=0.01 * resolution;
}