/home/pascal/depot/filedata/src/modisfiledata.h

/***************************************************************************
 *   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.             *
 ***************************************************************************/
#ifndef MODISFILEDATA_H
#define MODISFILEDATA_H

#include "hdffiledata.h"
#include "satellitefiledata.h"
#include "pixel.h"
#include "geometry.h"

typedef struct ModisCloudMask_1{
        unsigned int cloud_mask  : 1;
        unsigned int FOV_quality : 2;
        unsigned int day_night   : 1;
        unsigned int sunlight    : 1;
        unsigned int snow_ice    : 1;
        unsigned int land_water  : 2;
    } ModisCloudMask_1;
class MODISFileData : public SatelliteFileData, public HDFFileData {
  protected:
    typedef vector<double>                    GRingPoint;
    typedef vector<GRingPoint>                GRingPolygon;
    typedef vector<GRingPolygon>              GRing;
//     typedef P_Pixel_base< float,vector<int> > PixelType;
//     typedef vector<PixelType>::iterator       I_Pixel;

    static const string latitude_sds_name;
    static const string longitude_sds_name;
    static const int fast_srch_frame_x;
    static const int fast_srch_frame_y;
    float colocation_tolerance;
    enum PlatformType {
        UNDEFINED_PLATFORM_TYPE = 0,
        TERRA,
        AQUA,
    };
    PlatformType platform;

    enum Product_ID {
        UNDEFINED_PRODUCT_ID = 0,
        MYD02QKM,
        MYD02HKM,
        MYD021KM,
        MYD03,
        MYD02SSH,
        MYD04_L2,
        MYD05_L2,
        MYD06_L2,
        MOD06_L2,
    };
    Product_ID product_id;
    string product;
    float resolution;

    float lat_min;
    float lat_max;
    float lon_min;
    float lon_max;
    vector <int> lat_lon_index_max;
    GRing gring;
    vector<PixelType> v_pixel;
    vector < float32 > lat_gran_corner;
    vector < float32 > lon_gran_corner;
    void init();
    void set_time_coverage();
    void parse_filename(const string & file_basename);
    bool check_filename(const string & file_basename) const;
    void set_lat_lon_index_max();
    void set_gring();
    const bool is_inside_gring( const float & lat, const float & lon ) const;
    virtual string get_gring_latitude_name(const string &product="") {
        return "GRINGPOINTLATITUDE";
    };
    virtual string get_gring_longitude_name(const string &product="") {
        return "GRINGPOINTLONGITUDE";
    };
    virtual void load_v_pixel() ;

    void init_colocation_tolerance ();

public:
    static const int NB_MAX_X_10KM;
    static const int NB_MAX_Y_10KM;
    static const int NB_MAX_X_5KM;
    static const int NB_MAX_Y_5KM;
    static const int NB_MAX_X_1KM;
    static const int NB_MAX_Y_1KM;
    static const int NB_MAX_X_250M;
    static const int NB_MAX_Y_250M;

    const int get_x_size_10km() {
        switch ( product_id ) {
            case ( MYD04_L2 ) :
                return ( lat_lon_index_max[1] + 1 );
                break;
            default : // TODO set it for all products
                string msg ( "Unimplemented method for product "  + product ) ;
                g_exception e( __FILE__ , __LINE__ , msg );
                throw e;
                break;
        }
        return - 1;
    }
    const int get_y_size_10km() {
        switch ( product_id ) {
            case ( MYD04_L2 ) :
                return ( lat_lon_index_max[0] + 1 );
                break;
            default : // TODO set it for all products
                string msg ( "Unimplemented method for product "  + product ) ;
                g_exception e( __FILE__ , __LINE__ , msg );
                throw e;
                break;
        }
        return - 1;
    }

    const int get_x_size_5km() {
        switch ( product_id ) {
            case ( MYD06_L2 ) :
            case ( MOD06_L2 ) :
            case ( MYD05_L2 ) :
            case ( MYD021KM ) :
                return ( lat_lon_index_max[1] + 1 );
                break;
            default : // TODO set it for all products
                string msg ( "Unimplemented method for product "  + product ) ;
                g_exception e( __FILE__ , __LINE__ , msg );
                throw e;
                break;
        }
        return - 1;
    }
    const int get_y_size_5km() {
        switch ( product_id ) {
            case ( MYD06_L2 ) :
            case ( MOD06_L2 ) :
            case ( MYD05_L2 ) :
            case ( MYD021KM ) :
                return ( lat_lon_index_max[0] + 1 );
                break;
            default : // TODO set it for all products
                string msg ( "Unimplemented method for product "  + product ) ;
                g_exception e( __FILE__ , __LINE__ , msg );
                throw e;
                break;
        }
        return - 1;
    }
    const int get_x_size_1km() {
        switch ( product_id ) {
            case ( MYD06_L2 ) :
            case ( MOD06_L2 ) :
            case ( MYD05_L2 ) :
                // 270 5km pixels
                return ( get_x_size_5km () * 5 ) + 4;
            case ( MYD021KM ) :
                // 271 5km pixels
                return ( get_x_size_5km () * 5 ) - 1;
            case ( MYD02QKM ) :
            case ( MYD03    ) :
                return ( lat_lon_index_max[1] + 1 );
                break;
            default : // TODO set it for all products
                string msg ( "Unimplemented method for product "  + product ) ;
                g_exception e( __FILE__ , __LINE__ , msg );
                throw e;
                break;
        }
        return - 1;
    }

    const int get_y_size_1km() {
        switch ( product_id ) {
            case ( MYD06_L2 ) :
            case ( MOD06_L2 ) :
            case ( MYD05_L2 ) :
            case ( MYD021KM ) :
                return ( get_y_size_5km () * 5 );
            case ( MYD02QKM ) :
            case ( MYD03    ) :
                return ( lat_lon_index_max[0] + 1 );
                break;
            default : // TODO set it for all products
                string msg ( "Unimplemented method for product "  + product ) ;
                g_exception e( __FILE__ , __LINE__ , msg );
                throw e;
                break;
        }
        return - 1;
    }

    const int get_x_size_250m() {
        switch ( product_id ) {
            case ( MYD02QKM ) :
                return 4 * get_x_size_1km();
                break;
            default : // TODO set it for all products
                string msg ( "Unimplemented method for product "  + product ) ;
                g_exception e( __FILE__ , __LINE__ , msg );
                throw e;
                return - 1;
                break;
        }
        return -1;
    };

    const int get_y_size_250m() {
        switch ( product_id ) {
            case ( MYD02QKM ) :
                return 4 * get_y_size_1km();
                break;
            default : // TODO set it for all products
                string msg ( "At this time, this method is only supported for MYD02QKM products" ) ;
                g_exception e( __FILE__ , __LINE__ , msg );
                throw e;
                break;
        }
        return -1;
    };

    MODISFileData(const string &_name, const string &mode="r");
    ~MODISFileData();
    void print_gring() const;

    const string get_radix();
    const bool get_5_to_1km_index(const float &lat, const float &lon, int * one_km_indexes, const int * five_km_indexes=NULL);
    const bool get_1km_to_250m_pixel( const float32 lat, const float32 lon, const int i_pix_1km[2],
                                      int i_pix_250m[2]);
    const bool get_5km_to_1km_pixel( const float32 lat, const float32 lon, const int i_pix_5km[2],
                                      int i_pix_1km[2]);
    void get_5_to_1km_area(const int * five_km_idx, int * one_km_min_idx, int * one_km_max_idx);
    const int convert_5_to_1km_index(const int &five_km_index, const bool is_x_axis = false) const;

    void convert_1km_to_250m_index( const int  i_pix_1km  [2],
                                          float i_pix_250m [2] ) const;
    const bool contain_data(const float &lat,const float &lon, const double &time, const double &tolerance=0.1);
    const bool contain_location(const float &lat,const float &lon, const double &tolerance=0.1) ;
    const vector<int> get_lat_lon_index_max_5km() const {
        return lat_lon_index_max;
    }

    const vector<int> get_lat_lon_index_max_1km() const {
        vector<int> vect;
        vect.push_back(convert_5_to_1km_index(lat_lon_index_max[0],false)+2);
        vect.push_back(convert_5_to_1km_index(lat_lon_index_max[1],true)+1);
        return vect;
    }

    string get_product() const {
        return product;
    }
    float get_resolution() const {
        return resolution;
    }
    void load_geolocation_data();
    void free_geolocation_data();
    virtual void close_data_file() {
        free_hdf_file();
    };
    virtual void open_data_file() {
        load_hdf_file();
    };
    const bool is_geolocation_data_loaded() const;
    void get_sds_calibration(const string &sds_name, double & scale, double & offset, const string &scale_name="scale_factor", const string &offset_name="add_offset", const int channel_nb=-1);
    template<class Value_T,class Count_T>
    Value_T* read_calibrated_data(Value_T* data,const char* sds_name = "Height", int start[]=NULL, int stride[]=NULL, int edges[]=NULL, int rank=-1, Value_T fill_value =Value_T(0), const string &scale_name="scale_factor", const string &offset_name="add_offset", const int channel_nb=-1);
    const bool get_index(const float &lat, const float& lon, const double time, int * near_point_idx, float coloc_tolerance=-1.);

    const bool get_index(const float &lat, const float& lon, int & nearest_pix_idx, float coloc_tolerance=-1.);
     virtual void get_vindex(vector < vector < int > > &v_index, const float &lat, const float& lon,
                                            const float colocation_tolerance ) {
        UnimplementedMethod e(__FILE__,__LINE__,__PRETTY_FUNCTION__);
        throw e;
     };
    const float get_nearest_point_distance(const float &lat,const float &lon, float coloc_tolerance=-1.);
    virtual void get_pixel_coord ( const vector < int > & ipix, float &lat, float &lon, double &time ) {
        UnimplementedMethod e(__FILE__,__LINE__,__PRETTY_FUNCTION__);
        throw e;
    };
    vector< float32 > get_lat_gran_corner() const
    {
            return lat_gran_corner;
    }
    vector< float32 > get_lon_gran_corner() const
    {
            return lon_gran_corner;
    }

    vector< int > get_geolocation_dim() const
    {
        vector < int > v_dim (2);
        v_dim [0] = lat_lon_index_max [0] + 1;
        v_dim [1] = lat_lon_index_max [1] + 1;
        return v_dim;
    }

    inline void get_250m_pix_pos (  const int itk_250m, const int isc_250m,
                                    float32 & lat, float32 & lon );

    inline void get_1km_pix_pos (  const int itk_1km, const int isc_1km,
                                   float & lat, float & lon );

    inline void get_bounds_1km_to_250m ( const int itk_1km, const int isc_1km,
                                         int & itk_250m_min, int & itk_250m_max,
                                         int & isc_250m_min, int & isc_250m_max
                                        );

    inline void get_bounds_5km_to_1km ( const int itk_5km, const int isc_5km,
                                         int & itk_1km_min, int & itk_1km_max,
                                         int & isc_1km_min, int & isc_1km_max
                                        );

protected:

    void set_lat_lon_min_max();

    inline void get_y_pos_1km_to_250m ( const int isc_1km,
                                        const int itk_p1_1km, const int itk_p2_1km,
                                        const int itk_250m,
                                        double & lat, double & lon );

    inline void get_x_pos_1km_to_250m ( const double lat_left_250m,  const double lon_left_250m,
                                        const int isc_left_1km,
                                        const double lat_right_250m, const double lon_right_250m,
                                        const int isc_right_1km,
                                        const int isc_250m,
                                        double & lat, double & lon );

    inline void get_y_pos_5km_to_1km (  const int isc_5km,
                                        const int itk_p1_5km, const int itk_p2_5km,
                                        const int itk_1km,
                                        double & lat, double & lon );

    inline void get_x_pos_5km_to_1km (  const double lat_left_1km,  const double lon_left_1km,
                                        const int isc_left_5km,
                                        const double lat_right_1km, const double lon_right_1km,
                                        const int isc_right_5km,
                                        const int isc_1km,
                                        double & lat, double & lon );

};

template<class Value_T,class Count_T>
Value_T * MODISFileData::read_calibrated_data( Value_T * data, const char * sds_name, int start[], int stride[], int edges[], int rank, Value_T fill_value, const string & scale_name, const string & offset_name, const int channel_nb )
{
    bool are_limits_initialized[3]; // check if start, stride or edges had been initialized in this method
    bool hdf_file_already_loaded = is_hdf_file_loaded();
    try {
        if (!hdf_file_already_loaded)
            load_hdf_file();
        // --- initialize NULL input parametres ---
        init_read_write_null_input_param(sds_name, start, stride, edges,rank, are_limits_initialized);
        // check the selection limits
        check_read_write_limits(sds_name, start,stride,edges,rank); // throw bad_index in case of problem
        int nb_data = 1;
        int i=0;
        for (i = 0 ; i<rank ; ++i)
            nb_data*=edges[i];
        // read the uncalibrated dataS
        Count_T* count_data=NULL;
        count_data=static_cast<Count_T*>(read_data(count_data,sds_name,start,stride,edges,rank));
        // eventually allocate the calibrated data
        if (data==NULL) {
            // read the size of the sds data
            size_t data_type_size = hdf_file->get_sds_sizeof(sds_name);
            data = static_cast<Value_T*>(malloc(data_type_size*nb_data));
        }
        // read the calibration
        double cal,offset;
        get_sds_calibration(sds_name,cal,offset,scale_name,offset_name,channel_nb);
        // read the fill value
        Count_T count_fill_value;
        get_dataset_fill_value(sds_name,static_cast<void*>(&count_fill_value));
        // apply the calibration
        for (i=0;i<nb_data;++i) {
            if (count_data[i]!=count_fill_value)
                // a valid value
                data[i]=static_cast<Value_T>(cal)*(static_cast<Value_T>(count_data[i])-static_cast<Value_T>(offset));
            else
                // set a fill value
                data[i]=fill_value;
        }

       // free the ressources allocated by this method
        free_read_write_allocations(are_limits_initialized, start, stride, edges);
        delete[] count_data;
        if (!hdf_file_already_loaded)
            free_hdf_file();
        return data;

    } catch(exception &e) {
        // free the ressources allocated in this method
        //             free_read_write_allocations(are_limits_initialized , start, stride,edges);
        if (!hdf_file_already_loaded)
            free_hdf_file();
        cerr<<e.what()<<endl;
        throw;
    } catch(...) {
        // free the ressources allocated in this method
        //             free_read_write_allocations(are_limits_initialized , start, stride,edges);
        cout<< "Unknown exception in "<<__FILE__<<" line : "<<__LINE__<<endl;
        throw;
    }

}

void MODISFileData::get_y_pos_1km_to_250m ( const int isc_1km,
                                            const int itk_p1_1km, const int itk_p2_1km,
                                            const int itk_250m,
                                            double & lat, double & lon )
{
    // make sure P1 and P1 are 2 successive pixels along-track
    if ( fabs ( itk_p1_1km - itk_p2_1km ) != 1. )
    {
        string msg ( "P1 and P2 are on the same cross-track line and must be successive" ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }

    // lat, lon of the 1km bounding pixels
    int sz_sc_1km = get_x_size_1km();
    int i_buf = itk_p1_1km * sz_sc_1km + isc_1km;
    float32 p1_1km_lat = lat_data [ i_buf ];
    float32 p1_1km_lon = lon_data [ i_buf ];
    i_buf = itk_p2_1km * sz_sc_1km + isc_1km;
    float32 p2_1km_lat = lat_data [ i_buf ];
    float32 p2_1km_lon = lon_data [ i_buf ];

    // check for change date meridian particular case

    // change date meridian case sentinel
//     bool is_crossing_change_date = false;
    if ( fabs ( p1_1km_lon - p2_1km_lon ) > 180. )
    {
        if ( p1_1km_lon < 0. )
            p1_1km_lon += 360.;
        else if ( p2_1km_lon < 0. )
            p2_1km_lon += 360.;
//         is_crossing_change_date = true;
    }

    // coordinates of p1, p2 in the 250m grid
    double itk_p1_250m = 1.5 + 4. * itk_p1_1km;
    double itk_p2_250m = 1.5 + 4. * itk_p2_1km;

//     printf ( "itk_p1_250m=%f itk_p2_250m=%f\n", itk_p1_250m, itk_p2_250m );

    // linear interpolation on the position
    double div = 1. / ( itk_p2_250m - itk_p1_250m ); // for speed
    double dtk = ( itk_250m - itk_p1_250m );
    double alpha_lon = ( p2_1km_lon - p1_1km_lon ) * div;
    double alpha_lat = ( p2_1km_lat - p1_1km_lat ) * div;
    lon = p1_1km_lon + alpha_lon * dtk;
    lat = p1_1km_lat + alpha_lat * dtk;

    // in case of change date crossing, turn values > 180. to negative
//     if ( is_crossing_change_date and lon > 180. )
//         lon -= 360.;
}

void MODISFileData::get_x_pos_1km_to_250m ( const double lat_left_250m,  const double lon_left_250m,
                                            const int isc_left_1km,
                                            const double lat_right_250m, const double lon_right_250m,
                                            const int isc_right_1km,
                                            const int isc_250m,
                                            double & lat, double & lon )
{
    // make sure P1 and P1 are 2 successive pixels along-track
    if ( fabs ( isc_left_1km - isc_right_1km ) != 1. )
    {
        string msg ( "The 2 borders are on the same along-track line and must be successive" ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }

    // coordinates of left and right border in the 250m grid
    int isc_left_250m  = 4 * isc_left_1km;
    int isc_right_250m = 4 * isc_right_1km;

    //print "isc_left_1km=%d isc_right_1km=%d"%(isc_left_1km, isc_right_1km)
    //print "isc_left_250m=%d isc_right_250m=%d"%(isc_left_250m, isc_right_250m)
    //print "isc_250m_min=%d isc_250m_max=%d"%(isc_250m_min,isc_250m_max)

    // linear interpolation on the position
    double div = 1. / double ( isc_right_250m - isc_left_250m );
    double dsc = isc_250m - isc_left_250m;
    double alpha_lon = ( lon_right_250m - lon_left_250m ) * div;
    double alpha_lat = ( lat_right_250m - lat_left_250m ) * div;

    lat = lat_left_250m + alpha_lat * dsc;
    lon = lon_left_250m + alpha_lon * dsc;
}

void MODISFileData::get_250m_pix_pos (  const int itk_250m, const int isc_250m,
                                        float32 & lat, float32 & lon )
{
    if ( product_id != MYD02QKM )
    {
        string msg ( "At this time, this method is only supported for MYD02QKM products" ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }

    // change date meridian case sentinel
//     bool is_crossing_change_date = false;

    // check 250m indexes validity
    int sz_tk_250m = get_y_size_250m ();
    int sz_sc_250m = get_x_size_250m ();
    if ( ( isc_250m < 0 ) or ( isc_250m > ( sz_sc_250m - 1 ) ) )
    {
        string msg ( "Invalid cross-track index " + MyTools::to_string ( isc_250m ) +
                     ". Must in range [0," + MyTools::to_string ( sz_sc_250m - 1 ) ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }
    if ( ( itk_250m < 0 ) or ( itk_250m > ( sz_tk_250m - 1 ) ) )
    {
        string msg ( "Invalid along-track index " + MyTools::to_string ( itk_250m ) +
                     ". Must in range [0," + MyTools::to_string ( sz_tk_250m - 1 ) ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }

    // --- set the bounding 1km pixels to take for interpolation
    // set the (track,scan) indexes of the 1km pixel in the 250m grid ( can be a float grid index )
    double itk_1km = ( itk_250m - 1.5 ) / 4.;
    double isc_1km = isc_250m / 4.;

//     printf ( "i_250m=[%ld, %ld] -> i_1km=[%.2f, %.2f]", itk_250m, isc_250m, itk_1km, isc_1km );

    // the width of one scan, in number of pixels, at 250m resolution
    int w_scan_250m = 40;
    // - extra/interpolation 1km pixels along track
    double itk_top_1km    = -1.;
    double itk_bottom_1km = -1.;
    if (   ( itk_250m % w_scan_250m ) <= 1 ) { // extrapolate start of scan
        itk_top_1km    = ceil ( itk_1km );
        itk_bottom_1km = itk_top_1km + 1.;
    } else if ( ( itk_250m % w_scan_250m ) >= 38 ) { // extrapolate end of scan
        itk_top_1km    = floor ( itk_1km );
        itk_bottom_1km = itk_top_1km - 1.;
    } else { // general case : middle of scan
        itk_top_1km    = floor ( itk_1km );
        itk_bottom_1km = itk_top_1km + 1.;
    }

    // - extra/interpolation 1km pixels cross track
    double isc_left_1km  = -1.;
    double isc_right_1km = -1.;
    if ( isc_1km >= 1353. ) { // extrapolate end of scan line
        isc_left_1km  = floor ( isc_1km ) - 1.;
        isc_right_1km = floor ( isc_1km );
    } else { // general case : interpolation
        isc_left_1km  = floor ( isc_1km );
        isc_right_1km = isc_left_1km + 1.;
    }

//     printf ( "itk_top_1km=%f itk_bottom_1km=%f isc_left_1km=%f isc_right_1km=%f\n", itk_top_1km, itk_bottom_1km, isc_left_1km, isc_right_1km );

    // --- set the 1km track lines position ; left border ---
    double lat_left_250m, lon_left_250m;
    get_y_pos_1km_to_250m (
                ( int ) ( isc_left_1km ), ( int ) ( itk_top_1km ), ( int ) ( itk_bottom_1km ),
                ( int ) ( itk_250m ),
                lat_left_250m, lon_left_250m );
    // --- set the 1km track lines position ; right border ---
    double lat_right_250m, lon_right_250m;
    get_y_pos_1km_to_250m (
                ( int ) ( isc_right_1km ), ( int ) ( itk_top_1km ), ( int ) ( itk_bottom_1km ),
                ( int ) ( itk_250m ),
                lat_right_250m, lon_right_250m );

//     printf ( "left_250m=[%f,%f] right_250m=[%f,%f]\n",lat_left_250m, lon_left_250m, lat_right_250m, lon_right_250m);

    // check for change date meridian case
    if  ( fabs ( lon_right_250m  - lon_left_250m  ) > 180. ) {
//         is_crossing_change_date = true;
        // all negative longitudes will be incremented of 360 before interpolation
        if ( lon_left_250m < 0. )
            lon_left_250m += 360.;
        if ( lon_right_250m < 0. )
            lon_right_250m += 360.;
    }

    // for each track line position, interpolate along scan to retrieve the 250m geolocation
    double dbl_lat, dbl_lon;
    get_x_pos_1km_to_250m ( lat_left_250m,  lon_left_250m,  ( int ) ( isc_left_1km ),
                            lat_right_250m, lon_right_250m, ( int ) ( isc_right_1km ),
                            isc_250m,
                            dbl_lat, dbl_lon );

    // in case of change date crossing, turn values > 180. to negative
    if ( dbl_lon > 180. )
        dbl_lon -= 360.;
    if ( dbl_lon < -180. )
        dbl_lon += 360.;

    lat = ( float32 ) ( dbl_lat );
    lon = ( float32 ) ( dbl_lon );
//     printf ( "geolocation = [%f, %f]\n",lat,lon);

}

void MODISFileData::get_bounds_1km_to_250m( const int itk_1km, const int isc_1km,
                                            int & itk_250m_min, int & itk_250m_max,
                                            int & isc_250m_min, int & isc_250m_max )
{
    // set the (track,scan) indexes of the 1km pixel in the 250m grid
    float itk_250m = 1.5 + 4 * itk_1km;
    int  isc_250m =       4 * isc_1km;

    // the width of one scan, at 1km resolution
    int sz_sc_1km = get_x_size_1km();
    // set the interpolation quarters characteristics
    itk_250m_min    = ( int ) ( itk_250m - 1.5 );
    itk_250m_max    = ( int ) ( itk_250m + 1.5 );
    if ( isc_1km == 0 ) {
        isc_250m_min  = 0;
        isc_250m_max  = isc_250m + 2;
    } else if ( isc_1km == sz_sc_1km - 1 ) { // end of scan -> extrapolation along scan
        isc_250m_min  = isc_250m - 2;
        isc_250m_max  = isc_250m + 3;
    } else { // general case : 2 interpolations done along scan : [isc-1, isc] then [isc, isc+1]
        isc_250m_min  = isc_250m - 2;
        isc_250m_max  = isc_250m + 2;
    }
}

void MODISFileData::get_1km_pix_pos (  const int itk_1km, const int isc_1km,
                                       float & lat, float & lon )
{
    if ( product_id != MYD06_L2 && product_id != MOD06_L2 &&
         product_id != MYD05_L2 && product_id != MYD021KM )
    {
        string msg ( "Method not implemented for products " + product ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }

    // check 1km indexes validity
    int sz_tk_1km = get_y_size_1km ();
    int sz_sc_1km = get_x_size_1km ();
    if ( ( isc_1km < 0 ) || ( isc_1km > ( sz_sc_1km - 1 ) ) )
    {
        string msg ( "Invalid cross-track index " + MyTools::to_string ( isc_1km ) +
                     ". Must in range [0," + MyTools::to_string ( sz_sc_1km - 1 ) ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }
    if ( ( itk_1km < 0 ) || ( itk_1km > ( sz_tk_1km - 1 ) ) )
    {
        string msg ( "Invalid along-track index " + MyTools::to_string ( itk_1km ) +
                     ". Must in range [0," + MyTools::to_string ( sz_tk_1km - 1 ) ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }

    // --- set the bounding 5km pixels to take for interpolation
    // set the (track,scan) indexes of the 5km pixel in the 1km grid ( can be a float grid index )
    double itk_5km = ( itk_1km - 2. ) / 5.;
    double isc_5km = ( isc_1km - 2. ) / 5.;

//     printf ( "i_1km=[%ld, %ld] -> i_5km=[%.2f, %.2f]", itk_1km, isc_1km, itk_5km, isc_5km );

    // the width of one scan, in number of pixels, at 1km resolution
    int w_scan_1km = 10;
    // - extra/interpolation 5km pixels along track
    double itk_top_5km    = -1.;
    double itk_bottom_5km = -1.;

    if ( ( itk_1km % w_scan_1km ) <= 2 ) { // extrapolate start of scan
        itk_top_5km    = ceil ( itk_5km );
        itk_bottom_5km = itk_top_5km + 1;
    } else if ( ( itk_1km % w_scan_1km ) >= 7 ) { // extrapolate end of scan
        itk_top_5km    = floor ( itk_5km );
        itk_bottom_5km = itk_top_5km - 1;
    } else { // general case { middle of scan
        itk_top_5km    = floor ( itk_5km );
        itk_bottom_5km = itk_top_5km + 1;
    }
    // - extra/interpolation 5km pixels cross track
    double isc_left_5km  = -1.;
    double isc_right_5km = -1.;
    int sz_sc_5km = get_x_size_5km();
    if ( isc_1km <= 2 ) { // extrapolate start of scan line
        isc_left_5km  = 0;
        isc_right_5km = 1;
    } else if ( isc_5km >= ( sz_sc_5km - 1 )  ) { // extrapolate end of scan line
        isc_left_5km  = sz_sc_5km - 2;
        isc_right_5km = sz_sc_5km - 1;
    } else { // general case { interpolation
        isc_left_5km  = floor ( isc_5km );
        isc_right_5km = isc_left_5km + 1;
    }

//     printf ( "itk_top_5km=%f itk_bottom_5km=%f isc_left_5km=%f isc_right_5km=%f\n", itk_top_5km, itk_bottom_5km, isc_left_5km, isc_right_5km );

    // --- set the 5km track lines position ; left border ---
    double lat_left_1km, lon_left_1km;
    get_y_pos_5km_to_1km (
                ( int ) ( isc_left_5km ), ( int ) ( itk_top_5km ), ( int ) ( itk_bottom_5km ),
                ( int ) ( itk_1km ),
                lat_left_1km, lon_left_1km );
    // --- set the 5km track lines position ; right border ---
    double lat_right_1km, lon_right_1km;
    get_y_pos_5km_to_1km (
                ( int ) ( isc_right_5km ), ( int ) ( itk_top_5km ), ( int ) ( itk_bottom_5km ),
                ( int ) ( itk_1km ),
                lat_right_1km, lon_right_1km );

//     printf ( "left_1km=[%f,%f] right_1km=[%f,%f]\n",lat_left_1km, lon_left_1km, lat_right_1km, lon_right_1km);

    // check for change date meridian case
    if  ( fabs ( lon_right_1km  - lon_left_1km  ) > 180. ) {
        // all negative longitudes will be incremented of 360 before interpolation
        if ( lon_left_1km < 0. )
            lon_left_1km += 360.;
        if ( lon_right_1km < 0. )
            lon_right_1km += 360.;
    }

    // for each track line position, interpolate along scan to retrieve the 1km geolocation
    double dbl_lat, dbl_lon;
    get_x_pos_5km_to_1km ( lat_left_1km,  lon_left_1km,  ( int ) ( isc_left_5km ),
                            lat_right_1km, lon_right_1km, ( int ) ( isc_right_5km ),
                            isc_1km,
                            dbl_lat, dbl_lon );
    lat = ( float ) ( dbl_lat );
    lon = ( float ) ( dbl_lon );

    // in case of change date crossing, turn values > 180. to negative
    if ( lon > 180. )
        lon = lon - 360.;
    else if ( lon < -180. )
        lon += 360.;
}

void MODISFileData::get_y_pos_5km_to_1km (  const int isc_5km,
                                            const int itk_p1_5km, const int itk_p2_5km,
                                            const int itk_1km,
                                            double & lat, double & lon )
{
    // make sure P1 and P1 are 2 successive pixels along-track
    if ( fabs ( itk_p1_5km - itk_p2_5km ) != 1. )
    {
        string msg ( "P1 and P2 are on the same cross-track line and must be successive" ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }

    // lat, lon of the 5km bounding pixels
    int sz_sc_5km = get_x_size_5km();
    int i_buf = itk_p1_5km * sz_sc_5km + isc_5km;
    float32 p1_5km_lat = lat_data [ i_buf ];
    float32 p1_5km_lon = lon_data [ i_buf ];
    i_buf = itk_p2_5km * sz_sc_5km + isc_5km;
    float32 p2_5km_lat = lat_data [ i_buf ];
    float32 p2_5km_lon = lon_data [ i_buf ];

    // check for change date meridian particular case

    // change date meridian case sentinel
    if ( fabs ( p1_5km_lon - p2_5km_lon ) > 180. )
    {
        if ( p1_5km_lon < 0. )
            p1_5km_lon += 360.;
        if ( p2_5km_lon < 0. )
            p2_5km_lon += 360.;
    }

    // coordinates of p1, p2 in the 1km grid
    double itk_p1_1km = 2. + 5. * itk_p1_5km;
    double itk_p2_1km = 2. + 5. * itk_p2_5km;

//     printf ( "itk_p1_1km=%f itk_p2_1km=%f\n", itk_p1_1km, itk_p2_1km );

    // linear interpolation on the position
    double div = 1. / ( itk_p2_1km - itk_p1_1km ); // for speed
    double dtk = ( itk_1km - itk_p1_1km );
    double alpha_lon = ( p2_5km_lon - p1_5km_lon ) * div;
    double alpha_lat = ( p2_5km_lat - p1_5km_lat ) * div;
    lon = p1_5km_lon + alpha_lon * dtk;
    lat = p1_5km_lat + alpha_lat * dtk;

    // in case of change date crossing, turn values > 180. to negative and < 180. to positive
    if ( lon > 180. )
        lon -= 360.;
    else if ( lon < -180. )
        lon += 360.;
}

void MODISFileData::get_x_pos_5km_to_1km (  const double lat_left_1km,  const double lon_left_1km,
                                            const int isc_left_5km,
                                            const double lat_right_1km, const double lon_right_1km,
                                            const int isc_right_5km,
                                            const int isc_1km,
                                            double & lat, double & lon )
{
    // make sure P1 and P1 are 2 successive pixels along-track
    if ( fabs ( isc_left_5km - isc_right_5km ) != 1. )
    {
        string msg ( "The 2 borders are on the same along-track line and must be successive" ) ;
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }

    // coordinates of left and right border in the 1km grid
    int isc_left_1km  = 2 + 5 * isc_left_5km;
    int isc_right_1km = 2 + 5 * isc_right_5km;

    //print "isc_left_5km=%d isc_right_5km=%d"%(isc_left_5km, isc_right_5km)
    //print "isc_left_1km=%d isc_right_1km=%d"%(isc_left_1km, isc_right_1km)
    //print "isc_1km_min=%d isc_1km_max=%d"%(isc_1km_min,isc_1km_max)

    // linear interpolation on the position
    double div = 1. / double ( isc_right_1km - isc_left_1km );
    double dsc = isc_1km - isc_left_1km;
    double alpha_lon = ( lon_right_1km - lon_left_1km ) * div;
    double alpha_lat = ( lat_right_1km - lat_left_1km ) * div;

    lat = lat_left_1km + alpha_lat * dsc;
    lon = lon_left_1km + alpha_lon * dsc;
}

void MODISFileData::get_bounds_5km_to_1km(  const int itk_5km, const int isc_5km,
                                            int & itk_1km_min, int & itk_1km_max,
                                            int & isc_1km_min, int & isc_1km_max )
{
    // set the (track,scan) indexes of the 5km pixel in the 1km grid
    int itk_1km = 2 + 5 * itk_5km;
    int isc_1km = 2 + 5 * isc_5km;

    // the width of one scan, at 5km resolution
    int sz_sc_5km = get_x_size_5km();
    // set the interpolation quarters characteristics
    itk_1km_min    = itk_1km - 2;
    itk_1km_max    = itk_1km + 2;

    // general case : 2 interpolations done along scan : [isc-1, isc] then [isc, isc+1]
    isc_1km_min  = isc_1km - 2;
    isc_1km_max  = isc_1km + 2;
    // treat scan borders cases
    if ( isc_5km == 0 ) {  // start of scan -> extrapolation along scan
        isc_1km_min  = 0;
        isc_1km_max  = isc_1km + 2;
    } else if ( isc_5km == sz_sc_5km - 1 ) { // end of scan -> extrapolation along scan
        if ( product_id == MYD021KM )
            // This product have 271 5km pixels along scan, and so the last along scan pixel is composed of 4 1km pixels
            isc_1km_max  = isc_1km + 1;
        else
            // The others products are supposed by default to have 270 5km pixels along scan, and so the last along scan pixel is composed of 9 1km pixels to interpolate. This is the case for MYD06 and MYD05 products
            isc_1km_max  = isc_1km + 6;
    }
}

#endif