/***************************************************************************
 *   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"

/**
 * @struct CloudMask_1 facility for cloud mask's flags reading. Equivalent to a int8 type.
Read and uncompress the cloud mask first byte (that's what "_1" means)
*/
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 managing the opening, reading, and accessing to the data of a MODIS File.
MODIS products managed : MYD06_L2 ; MYD021KM ;
@author Nicolas PASCAL
*/
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;

    /** name of the geolocation latitudes dataset */
    static const string latitude_sds_name;
    /** name of the geolocation longitudes dataset */
    static const string longitude_sds_name;
    /** x size to do a fast search when searching for a coincidence and when a previous coincidence in the granule has already been computed. Optimization purpose*/
    static const int fast_srch_frame_x;
    /** y size to do a fast search when searching for a coincidence and when a previous coincidence in the granule has already been computed. Optimization purpose*/
    static const int fast_srch_frame_y;
    /** default maximal distance for which to pixels are considered as coincident. Depends of the product */
    float colocation_tolerance;
    /**
     * @enum PlatformType defines the possible types of MODIS platform
     */
    enum PlatformType {
        UNDEFINED_PLATFORM_TYPE = 0,
        TERRA,
        AQUA,
    };
    /** Platform for this product */
    PlatformType platform;

    /**
     * @enum Product_ID defines an ID for the different supported MODIS product
     */
    enum Product_ID {
        UNDEFINED_PRODUCT_ID = 0,
        MYD02QKM,
        MYD02HKM,
        MYD021KM,
        MYD03,
        MYD02SSH,
        MYD04_L2,
        MYD05_L2,
        MYD06_L2,
        MOD06_L2,
    };
    /** integer that represents the type of MODIS product */
    Product_ID product_id;
    /** name of the product */
    string product;
    /** resolution of the geolocation data in km */
    float resolution;

    /** minimal bounding latitude of the granule */
    float lat_min;
    /** maximal bounding latitude of the granule */
    float lat_max;
    /** minimal bounding longitude of the granule */
    float lon_min;
    /** maximal bounding longitude of the granule */
    float lon_max;
    /** Maximal valid indexes of the geolocation datasets along [Y,X] */
    vector <int> lat_lon_index_max;
    /** list of the polygons that compose the GRing (each polygon is a list of points) */
    GRing gring;
    /** list of pixels data, sorted by increasing latitude, then increasing longitude */
    vector<PixelType> v_pixel;
    /** latitudes of the granule corners as read in the geolocation dataset, in the order [0,0], [0,I_SCAN_MAX], [J_TRACK_MAX, 0], [J_TRACK_MAX, I_SCAN_MAX] */
    vector < float32 > lat_gran_corner;
    /** longitudes of the granule corners as read in the geolocation dataset, in the order [0,0], [0,I_SCAN_MAX], [J_TRACK_MAX, 0], [J_TRACK_MAX, I_SCAN_MAX] */
    vector < float32 > lon_gran_corner;
    /**
     * @brief Initialize some basic attributes
     */
    void init();
    /**
     * @brief sets the time coverage of the product
     */
    void set_time_coverage();
    /**
     * @brief initializes some characteristics based on the product name
     * @param file_basename the basename of the product file
     */
    void parse_filename(const string & file_basename);
    /**
     * @brief check if the product name matches the official naming convention
     * @param file_basename the basename of the product file
     * @return true if the file name matches the convention
     */
    bool check_filename(const string & file_basename) const;
    /**
     * @brief sets the size of the geolocation dataset
     */
    void set_lat_lon_index_max();
    /**
     * @brief set the GRing coordinates
     * Read in the coremetadata
     */
    void set_gring();
    /**
     * @brief test if the given point is inside the GRing
      * The grings which are crossing the change date longitude are treated
     * @param lat the latitude of the point
     * @param lon the longitude of the point
     * @return true if the point is inside the gring
     */
    const bool is_inside_gring( const float & lat, const float & lon ) const;
    /**
     * @brief access to the name of the field that contains the latitude gring
     * @param product the file product
     * @return the gring latitude name. Default : "GRINGPOINTLATITUDE"
     */
    virtual string get_gring_latitude_name(const string &product="") {
        return "GRINGPOINTLATITUDE";
    };
    /**
     * @brief access to the name of the field that contains the longitude gring
     * @param product the file product
     * @return the gring longitude name. Default : "GRINGPOINTLONGITUDE"
     */
    virtual string get_gring_longitude_name(const string &product="") {
        return "GRINGPOINTLONGITUDE";
    };
    /**
     * @brief load the list data pixels
     */
    virtual void load_v_pixel() ;

    /**
     * @brief defines the colocatiuon tolerance, depending of the product type, because the geolocation datasets are not always given at the same resolution
     */
    void init_colocation_tolerance ();

public:
    /**Is the maximum number of pixels the along X at 10Km resolution */
    static const int NB_MAX_X_10KM;
    /**Is the maximum number of pixels the along Y at 10Km resolution */
    static const int NB_MAX_Y_10KM;
    /**Is the maximum number of pixels the along X at 5Km resolution */
    static const int NB_MAX_X_5KM;
    /**Is the maximum number of pixels the along Y at 5Km resolution */
    static const int NB_MAX_Y_5KM;
    /**Is the maximum number of pixels the along X at 1Km resolution */
    static const int NB_MAX_X_1KM;
    /**Is the maximum number of pixels the along Y at 1Km resolution */
    static const int NB_MAX_Y_1KM;
    /**Is the maximum number of pixels the along X at 250m resolution */
    static const int NB_MAX_X_250M;
    /**Is the maximum number of pixels the along Y at 250m resolution */
    static const int NB_MAX_Y_250M;

    /**
     * @brief accessor to the size of the swath for a 10km resolution
     * @return the size of the swath at 10km
     */
    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;
    }
    /**
     * @brief accessor to the size of the track for a 10km resolution
     * @return the size of the track at 10km
     */
    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;
    }

    /**
     * @brief accessor to the size of the swath for a 5km resolution
     * @return the size of the swath at 5km
     */
    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;
    }
    /**
     * @brief accessor to the size of the track for a 5km resolution
     * @return the size of the track at 5km
     */
    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;
    }
    /**
     * @brief accessor to the size of the swath for a 1km resolution
     * @return the size of the swath at 1km
     */
    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;
    }

    /**
     * @brief accessor to the size of the track for a 1km resolution
     * @return the size of the track at 1km
     */
    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;
    }

    /**
     * @brief accessor to the size of the cross-track dimension of the 250m resolution grid
     * @warning only supported for MYD02QKM products at this time
     * @return the size of the cross-track dimension at 250m
     */
    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;
    };

    /**
     * @brief accessor to the size of the along-track dimension of the 250m resolution grid
     * @warning only supported for MYD02QKM products at this time
     * @return the size of the along-track dimension at 250m
     */
    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;
    };

    /**
     * @brief Constructor
     * @param _name the name (path+filename) of the file to be opened
     * @param mode the opening mode. LIMITATION : only "r" mode treated at this time
     */
    MODISFileData(const string &_name, const string &mode="r");
    /**
     * @brief Destructor :
     */
    ~MODISFileData();
    /**
     * @brief prints out the gring values
     */
    void print_gring() const;

    /**
     * @brief access to the radix of the MODIS file. $
     * What is called the radix is the filename ( without path ) with product and date informations only (and without dot also)
     * Ex : "MYD021KM.A2003281.1340.004.2004158015852.hdf" as a radix of "MYD021KM.A2003281"
     * @return the radix
     */
    const string get_radix();
    /**
     * @brief convert indexes at 5km to ones at 1km.
     * Uses a linear interpolation between the data values read at 5km, and the real desired lat, lon.
     * @param lat the latitude
     * @param lon the longitude
     * @param one_km_indexes the indexes at 1 km resolution (output)
     * @param five_km_indexes the indexes at 5 km resolution (input). Not necessary, but if not precised, the search of those indexes is proceed, and it takes time.
     */
    const bool get_5_to_1km_index(const float &lat, const float &lon, int * one_km_indexes, const int * five_km_indexes=NULL);
    /**
     * @brief get 250m indexes of nearest point when geolocation data are given at 1km.
     * @warning this method is supposed to be used only with MYD02QKM products
     * Uses a bilinear interpolation between the geolocation data at 1km, and the real desired lat, lon to obtain the 250m indexes
     * @param lat [IN] the latitude of the point to search
     * @param lon [IN] the longitude of the point to search
     * @param i_pix_1km [IN] the indexes at 1km resolution (input)
     * @param i_pix_250m [OUT] the indexes at 250m resolution
     */
    const bool get_1km_to_250m_pixel( const float32 lat, const float32 lon, const int i_pix_1km[2],
                                      int i_pix_250m[2]);
    /**
     * @brief get 1km indexes of nearest point when geolocation data are given at 5km.
     * @warning this method is tested only with MYD06 and MYD05 products
     * Uses a bilinear interpolation between the geolocation data at 5km, and the real desired lat, lon to obtain the 1km indexes
     * @param lat [IN] the latitude of the point to search
     * @param lon [IN] the longitude of the point to search
     * @param i_pix_1km [IN] the indexes at 5km resolution
     * @param i_pix_250m [OUT] the indexes at 1km resolution
     */
    const bool get_5km_to_1km_pixel( const float32 lat, const float32 lon, const int i_pix_5km[2],
                                      int i_pix_1km[2]);
    /**
     * @brief get the top-left and bottom-right corner of a 5km cell using 1km indexes.
     * The indexes start from 0. The resulting 1km indexes are included in the 5km area.
     * Precondition : the allocations MUST be already done
     * @param five_km_idx the indexes of a cell of 5km side, using 5km indexes
     * @param one_km_min_idx the top-left corner indexes (1km)
     * @param one_km_max_idx  the bottom-right corner indexes (1km)
     */
    void get_5_to_1km_area(const int * five_km_idx, int * one_km_min_idx, int * one_km_max_idx);
    /**
    * @brief convert an index at 1km to the one representing the same point at 5km.
    * @param five_km_index the index at 5 km resolution
    * @param is_x_axis wether the axis we search the index is the x or the y one
    * @return the index of the same point at at 1 km resolution
    */
    const int convert_5_to_1km_index(const int &five_km_index, const bool is_x_axis = false) const;

    /**
    * @brief convert the indexes of a point in the 250m grid to the corresponding indexes in the 1km grid.
    * @param i_pix_1km the (track, swath) indexes at 1km resolution
    * @param i_pix_250m the (track, swath) indexes at 250m resolution. In float type because of non integer offset between the 2 grids along track
    */
    void convert_1km_to_250m_index( const int  i_pix_1km  [2],
                                          float i_pix_250m [2] ) const;
    /**
     * check if this file has data coincident with (lat,lon,time)
     * @param lat the latitude of the event
     * @param lon the longitude of the event
     * @param time the time of the event. If -1., time test will be ignored
     * @return
     */
    const bool contain_data(const float &lat,const float &lon, const double &time, const double &tolerance=0.1);
    /**
     * @brief check if this file has eventually data coincident with (lat,lon)
     * Actually, it only tests if (lat,lon) is contained in the data's bounding rectangle.
     * @param lat the latitude of the event
     * @param lon the longitude of the event
     * @param tolerance acceptable bias between the nearest point in the data and the given (lat,lon) point
     * @return true if a point in the data has been found in the colocation frame
     */
    const bool contain_location(const float &lat,const float &lon, const double &tolerance=0.1) ;
    /**
     * @brief access to the higher usable indexes of the "Latitude" and "Longitude" sds.
     * In other words, it returns the size of those sds arrays, diminuted of 1. Ex : if the sds "Latitude" has a size of (287,3), it will return [286,2]
     * @return an array containing the maximal usable indexes for the geolocation sds
     */
    const vector<int> get_lat_lon_index_max_5km() const {
        return lat_lon_index_max;
    }

    /**
     * @brief access to the higher usable indexes of the "Latitude" and "Longitude" sds.
     * In other words, it returns the size of those sds arrays, diminuted of 1. Ex : if the sds "Latitude" has a size of (287,3), it will return [286,2]
     * @return an array containing the maximal usable indexes for the geolocation sds
     */
    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;
    }

    /**
     * @brief accessor to the product name (MYD06_L2,...)
     * @return the product name
     */
    string get_product() const {
        return product;
    }
    /**
     * @brief accessor to the resolution of th product (in Km)
     * @return the product's resolution (Km)
     */
    float get_resolution() const {
        return resolution;
    }
    /**
     * @brief read the geolocation data and put it in memory.
     * This method is used to make the search of the indexes of a (lat,lon,time) point faster.
     */
    void load_geolocation_data();
    /**
     * @brief free eventually loaded geolocation data
     */
    void free_geolocation_data();
   /**
    * @brief closes the file.
    */
    virtual void close_data_file() {
        free_hdf_file();
    };
    /**
    * @brief opens the file.
    */
    virtual void open_data_file() {
        load_hdf_file();
    };
    /**
     * @brief Tell if the geolocation data have been already loaded
     */
    const bool is_geolocation_data_loaded() const;
    /**
     * @brief read the calibration coefficients of the sds given in parametre
     * In most cases, you
     * @param sds_name the name of the sds where to read the calibration coefficients
     * @param scale the output scale
     * @param offset the output offset
     * @param scale_name the name of the scale attribute as given in the sds
     * @param offset_name the name of the offset attribute as given in the sds
     * @param channel_nb the channel number. Useful only for multiple scaling values in one sds
     */
    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);
    /**
     * @brief read the data of a MODIS file sds, and apply the calibration to it
     * The applied calibration uses the formula : scale*(x-add_offset)
     * The method manages the allocation of the data' array if data is NULL.
     * Here is assumed that once the calibration is applied, the float type fits the data. Perhaps, it will be nexessary to review it
     * The template arguments are :
     * - Value_T : the type of the calibrated values
     * - Count_T : the type of the count (ie uncalibrated) data
     * @warning this method has been used only with Value_T of float.
     * @param data the buffer to fill with the read values
     * @param sds_name the name of the sds (Scientific Data Set) we want to access.
     * @param start begining of the selection. If NULL, start at (0,0) if rank is 2 ; (0,0,0) if rank is 3...
     * @param stride step between 2 interesting values. If NULL, this step is set to 1 in each dimension (ie all values will be read)
     * @param edges number of values to be read in each dimension. if NULL, it will be all data along each dimension.
     * @param rank the dimension of start, stride and edges
     * @param fill_value the fill value to use in the calibrated data
     * @param scale_name the name of the scale attribute as define in the sds
     * @param offset_name the name of the offset attribute as define in the sds
     * @param channel_nb the index of the channel for scaling with multiple values
     * @return a pointer to the read data array. It is useful when you let this method manage the allocation to update the allocated pointer to the data. If you've done it yourself or you are using a static buffer, it isn't necessary to catch this return pointer
     */
    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);
    /**
     * @brief search the index of the nearest point in the data to (@a lat , @a lon )
     * If @a coloc_tolerance is given, it will compute only the distance of the points that have a distance to ( @a lat, @a lon ) inferior to @a coloc_tolerance
     * @param near_point_idx the indexes of the nearest point found (or  {-1,-1} if not found)
     * @param lat the latitude
     * @param lon the longitude
     * @param time the time of the observation. If -1., skipped
     * @param coloc_tolerance the -/+ maximal tolerance for 2 points considered as colocated
     * @return true if found
     */
    const bool get_index(const float &lat, const float& lon, const double time, int * near_point_idx, float coloc_tolerance=-1.);

     /**
     * @brief search the index of the nearest point in the data to (@a lat , @a lon )
     * If @a coloc_tolerance is given, it will compute only the distance of the points that have a distance to ( @a lat, @a lon ) inferior to @a coloc_tolerance
     * @param near_point_idx the indexes of the nearest point found (or  {-1,-1} if not found)
     * @param lat the latitude
     * @param lon the longitude
     * @param time the time of the observation. If -1., skipped
     * @param coloc_tolerance the -/+ maximal tolerance for 2 points considered as colocated
     * @return true if found
     */
    const bool get_index(const float &lat, const float& lon, int & nearest_pix_idx, float coloc_tolerance=-1.);
    /**
     * @brief build the list of indices of pixels that are in colocation tolerance, sorted by increasing distance to (lat,lon)
     * If (lat,lon) is not found or out of the colocalisation_frame, returns an empty vector
     * @param lat the latitude
     * @param lon the longitude
     * @param colocation_tolerance the acceptable bias (in km or degrees. Supposed to be in, a plane approximation) between [lat,lon] and the nearest data point.
     * @return the list of indices of pixels
     */
     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;
     };
    /**
     * @brief compute the distance to (lat,lon) of the nearest point in the data
     * If @a coloc_tolerance is given, it will compute only the distance of the points that have a distance to ( @a lat, @a lon ) inferior to @a coloc_tolerance
     * @param lat the latitude
     * @param lon the longitude
     * @param coloc_tolerance the -/+ maximal tolerance for 2 points considered as colocated
     * @return the distance to the nearest point, or -1 if no point in the colocalisation frame has been found.
     */
    const float get_nearest_point_distance(const float &lat,const float &lon, float coloc_tolerance=-1.);
    /**
     * @brief retrieve the coordinates of a pixel using its index
     * @param ipix [IN] index of the pixel
     * @param lat [OUT] latitude of the pixel
     * @param lon [OUT] longitude of the pixel
     * @param time [OUT] timestamp of the pixel
     */
    virtual void get_pixel_coord ( const vector < int > & ipix, float &lat, float &lon, double &time ) {
        UnimplementedMethod e(__FILE__,__LINE__,__PRETTY_FUNCTION__);
        throw e;
    };
    /**
     * @brief accessor to the latitudes of the granule corners
     * @return the latitudes of the granule corners
     */
    vector< float32 > get_lat_gran_corner() const
    {
            return lat_gran_corner;
    }
    /**
     * @brief accessor to the longitudes of the granule corners
     * @return the longitudes of the granule corners
     */
    vector< float32 > get_lon_gran_corner() const
    {
            return lon_gran_corner;
    }

    /**
     * @brief return the dimensions of the geolocation dataset, in [ npix_track, npix_scan ] order
     * @return the dimensions of the geolocation dataset
     */
    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;
    }

    /**
     * @brief compute the (lat,lon) of a 250m pixel specified with its indexes, when the geolocation datasets are given at 1km resolution
     * @warning at this time, this method is only supported for MYD02QKM products
     * @param itk_250m [IN] grid index of the 250m pixel along-track
     * @param isc_250m [IN] grid index of the 250m pixel cross-track
     * @param lat [OUT] the latitude of the 250m pixel center at [ itk_250m, isc_250m ]
     * @param lon [OUT] the longitude of the 250m pixel center at [ itk_250m, isc_250m ]
     */
    inline void get_250m_pix_pos (  const int itk_250m, const int isc_250m,
                                    float32 & lat, float32 & lon );

    /**
     * @brief compute the (lat,lon) of a 1km pixel specified with its indexes, when the geolocation datasets are given at 1km resolution
     * @warning at this time, this method is only supported for MYD06, MOD06, MYD05 products
     * @param itk_1km [IN] grid index of the 1km pixel along-track
     * @param isc_1km [IN] grid index of the 1km pixel cross-track
     * @param lat [OUT] the latitude of the 1km pixel center at [ itk_1km, isc_1km ]
     * @param lon [OUT] the longitude of the 1km pixel center at [ itk_1km, isc_1km ]
     */
    inline void get_1km_pix_pos (  const int itk_1km, const int isc_1km,
                                   float & lat, float & lon );

    /**
     * @brief return the bound indexes, in the 250m grid, of all the 250m pixels that are included in a 1km pixel
     * The 250m pixels that are situed at the border of a 1km pixel are also returned.
     * The last 250m pixel in the cross-track direction is also considered as included in the last 1km cross-track pixel
     * @warning at this time, this method is only supported for MYD02QKM products
     * @param itk_1km [IN] 1km grid index of the 1km pixel along-track
     * @param isc_1km [IN] 1km grid index of the 1km pixel cross-track
     * @param itk_250m_min [OUT] min 250m grid index along-track of the 250m pixels included in the footprint of the requested 1km pixel
     * @param itk_250m_max [OUT] max 250m grid index along-track of the 250m pixels included in the footprint of the requested 1km pixel
     * @param isc_250m_min [OUT] min 250m grid index cross-track of the 250m pixels included in the footprint of the requested 1km pixel
     * @param isc_250m_max [OUT] max 250m grid index cross-track of the 250m pixels included in the footprint of the requested 1km pixel
     */
    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
                                        );

    /**
     * @brief return the bound indexes, in the 1km grid, of all the 1km pixels that are included in a 5km pixel
     * The 1km pixels that are situed at the border of a 5km pixel are also returned.
     * The last 1km pixel in the cross-track direction is also considered as included in the last 5km cross-track pixel
     * @warning at this time, this method is only supported for MYD02QKM products
     * @param itk_5km [IN] 5km grid index of the 5km pixel along-track
     * @param isc_5km [IN] 5km grid index of the 5km pixel cross-track
     * @param itk_1km_min [OUT] min 1km grid index along-track of the 1km pixels included in the footprint of the requested 5km pixel
     * @param itk_1km_max [OUT] max 1km grid index along-track of the 1km pixels included in the footprint of the requested 5km pixel
     * @param isc_1km_min [OUT] min 1km grid index cross-track of the 1km pixels included in the footprint of the requested 5km pixel
     * @param isc_1km_max [OUT] max 1km grid index cross-track of the 1km pixels included in the footprint of the requested 5km pixel
     */
    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:

    /**
     * @brief sets the granule bounding coordinates
     * Minimal and maximal latitudes and longitudes are read in the coremetadata
     */
    void set_lat_lon_min_max();

    /**
     * @brief return the position of 250m pixel defined by its along-track index between 2 successive 1km pixels on a same cross-track line
     * @warning If the 2 pixels are crossing the changing date meridian, the negative longitude will be returned with an increment of 360 to manage correctly the interpolation
     * @param isc_1km [IN] index of the cross-track line in the 1km grid
     * @param itk_p1_1km [IN] index of the first along-track pixel in the 1km grid
     * @param itk_p2_1km [IN] index of the second along-track pixel in the 1km grid
     * @param itk_250m [IN] index of the along-track 250m along-track pixel to interpolate
     * @param lat [OUT] latitude of the pixel at [itk_250m, isc_1km], interpolated between the 2 successive 1km pixels p1 and p2
     * @param lon [OUT] longitude of the pixel at [itk_250m, isc_1km], interpolated between the 2 successive 1km pixels p1 and p2
     */
    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 );

    /**
     * @brief retrieve the position of a 250m pixel set by its cross-track index and the position of 2 along-track lines set by 2 sucessive 1km pixels
     * @param lat_left_250m [IN] latitude of the 250m along track pixel that is the intersection between cross-track line at index isc_250m and the along-track line that joins the 2 sucessive bounding 1km pixels on the left of the point to interpolate
     * @param lon_left_250m [IN] longitude of the 250m along track pixel that is the intersection between cross-track line at index isc_250m and the along-track line that joins the 2 sucessive bounding 1km pixels on the left of the point to interpolate
     * @param isc_left_1km [IN] cross-track index at 1km of the left border
     * @param lat_right_250m [IN] latitude of the 250m along track pixel that is the intersection between cross-track line at index isc_250m and the along-track line that joins the 2 sucessive bounding 1km pixels on the right of the point to interpolate
     * @param lon_right_250m [IN] longitude of the 250m along track pixel that is the intersection between cross-track line at index isc_250m and  the along-track line that joins the 2 sucessive bounding 1km pixels on the right of the point to interpolate
     * @param isc_right_1km [IN] cross-track index at 1km of the right border
     * @param lat [OUT] latitude of the pixel at [itk_250m, isc_1km], interpolated between the 2 successive 1km pixels p1 and p2
     * @param lon [OUT] longitude of the pixel at [itk_250m, isc_1km], interpolated between the 2 successive 1km pixels p1 and p2
     */
    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 );

    /**
     * @brief return the position of 1km pixel defined by its along-track index between 2 successive 5km pixels on a same cross-track line
     * @warning If the 2 pixels are crossing the changing date meridian, the negative longitude will be returned with an increment of 360 to manage correctly the interpolation
     * @param isc_5km [IN] index of the cross-track line in the 5km grid
     * @param itk_p1_5km [IN] index of the first along-track pixel in the 5km grid
     * @param itk_p2_5km [IN] index of the second along-track pixel in the 5km grid
     * @param itk_1km [IN] index of the along-track 1km along-track pixel to interpolate
     * @param lat [OUT] latitude of the pixel at [itk_1km, isc_5km], interpolated between the 2 successive 5km pixels p1 and p2
     * @param lon [OUT] longitude of the pixel at [itk_1km, isc_5km], interpolated between the 2 successive 5km pixels p1 and p2
     */
    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 );

    /**
     * @brief retrieve the position of a 1km pixel set by its cross-track index and the position of 2 along-track lines set by 2 sucessive 5km pixels
     * @param lat_left_1km [IN] latitude of the 1km along track pixel that is the intersection between cross-track line at index isc_1km and the along-track line that joins the 2 sucessive bounding 5km pixels on the left of the point to interpolate
     * @param lon_left_1km [IN] longitude of the 1km along track pixel that is the intersection between cross-track line at index isc_1km and the along-track line that joins the 2 sucessive bounding 5km pixels on the left of the point to interpolate
     * @param isc_left_5km [IN] cross-track index at 5km of the left border
     * @param lat_right_1km [IN] latitude of the 1km along track pixel that is the intersection between cross-track line at index isc_1km and the along-track line that joins the 2 sucessive bounding 5km pixels on the right of the point to interpolate
     * @param lon_right_1km [IN] longitude of the 1km along track pixel that is the intersection between cross-track line at index isc_1km and  the along-track line that joins the 2 sucessive bounding 5km pixels on the right of the point to interpolate
     * @param isc_right_5km [IN] cross-track index at 5km of the right border
     * @param lat [OUT] latitude of the pixel at [itk_1km, isc_5km], interpolated between the 2 successive 5km pixels p1 and p2
     * @param lon [OUT] longitude of the pixel at [itk_1km, isc_5km], interpolated between the 2 successive 5km pixels p1 and p2
     */
    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