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

/***************************************************************************
 *   Copyright (C) 2011 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 GRIBFILEDATA_H
#define GRIBFILEDATA_H

#include <cstdio>
#include <string>
#include <map>
#include "grib_api.h"

#include "filedatareader.h"
using namespace std;

#define MAX_KEY_LEN  255
#define MAX_VAL_LEN  1024

class t_grib_var_desc {
public:
    string name;
    string long_name;
    int param_id;
    string type_id;
    int sz_x;
    int sz_y;
    string z_level_type;
    double fill_value;
    double scale;
    double offset;
    vector <int> v_msg;
    vector <int> v_lvl;

    t_grib_var_desc (const string & name = "", const string & long_name = "", const int param_id = -1, const string & type_id = "",
        const int sz_x = -1, const int sz_y = -1, const string & z_level_type = "",
        const double fill_value = -DBL_MAX, const double scale = -DBL_MAX, const double offset = -DBL_MAX) :
        name (name),
        long_name (long_name),
        param_id (param_id),
        type_id (type_id),
        sz_x (sz_x),
        sz_y (sz_y),
        z_level_type (z_level_type),
        fill_value (fill_value),
        scale (scale),
        offset (offset) {
        // check pressure level validity
        if (!z_level_type.empty() &&
            z_level_type != string("surface") &&
            z_level_type != string("isobaricInhPa") &&
            z_level_type != string("hybrid") ) {
            string msg ( "Invalid Pressure level type " + z_level_type ) ;
            g_exception e( __FILE__ , __LINE__ , msg );
            throw e;
        }
        v_msg = vector <int> (0);
        v_lvl = vector <int> (0);
    };
    ~t_grib_var_desc () {
        v_msg.clear();
        v_lvl.clear();
    };
    friend std::ostream &operator<< (std::ostream &stream, const t_grib_var_desc &grib_var_desc) {
        stream << grib_var_desc.to_string();
        return stream;
    }
    const string to_string (void) const {
        string desc ("");
        desc += "var name     : " + name + "\n";
        desc += "long name    : " + long_name + "\n";
        desc += "param ID     : " + MyTools::to_string (param_id) + "\n";
        desc += "Data Type    : " + MyTools::to_string (type_id) + "\n";
        desc += "grid size    : " + MyTools::to_string (sz_y) + "x" + MyTools::to_string (sz_x) + "\n";
        desc += "Z level type : " + z_level_type + "\n";
        desc += "Fill value   : " + MyTools::to_string (fill_value) + "\n";
        desc += "Scale factor : " + MyTools::to_string (scale) + "\n";
        desc += "Scale offset : " + MyTools::to_string (offset) + "\n";
        desc += "Msg indices  : ";
        for (vector <int>::const_iterator it_msg = v_msg.begin(); it_msg != v_msg.end() ; ++it_msg )
            desc += MyTools::to_string (*it_msg) + " ";
        desc += "\n";
        desc += "Msg levels   : ";
        for (vector <int>::const_iterator it_lvl = v_lvl.begin(); it_lvl != v_lvl.end() ; ++it_lvl )
            desc += MyTools::to_string (*it_lvl) + " ";
        desc += "\n";
        return desc;
    };
    void add_msg (const int i_msg, const int z_lvl) {
        v_msg.push_back(i_msg);
        v_lvl.push_back(z_lvl);
    }
    vector <int> get_sz_dims (void) {
        vector <int> sz_dims (0);
        sz_dims.reserve (3);

        int sz_z = v_lvl.size();
        if (sz_z > 1) // 3D var
            sz_dims.push_back(sz_z);
        sz_dims.push_back(sz_y);
        sz_dims.push_back(sz_x);
        return sz_dims;
    }
};

class GribFileData : virtual public FileData {
protected:
    // --- parameters attributes IDs --- //
    static const string var_name_id;
    static const string parameter_code_id;
    static const string parameter_name_id;
    static const string fill_value_id;
    static const string scale_factor_id;
    static const string add_offset_id;
    static const string x_dim_id;
    static const string y_dim_id;
    static const string level_type_id;
    static const string level_id;
    static const string data_flag_id;
    static const string data_width_id;
    static const string data_id;

    FILE * fd;
//     grib_handle * handle;
    int n_msg;
//     int icur_msg;
    typedef map < string, t_grib_var_desc > t_var_table;
    t_var_table var_table;
public:
    /*********************
    **   CONSTRUCTORS   **
    *********************/
    GribFileData(const string &name = "", const string &mode = "r", const bool open_file = false);
    GribFileData(const GribFileData &fdata);
    virtual ~GribFileData();
    GribFileData &operator= (const GribFileData &fdata);
    void open_data_file();
    void close_data_file();
    const bool is_file_loaded();
    const bool is_existing_dataset (const string & ds_name) {
        return (var_table.find(ds_name) != var_table.end());
    };
    const t_grib_var_desc & get_dataset_desc (const string & ds_name) {
        return var_table [ds_name];
    };
    string get_msg_key (int i_msg, const string & key);
    virtual void get_dataset_fill_value(const string &dataset_name, void* fill_value);
    vector<int> get_dataset_dimension (const string &dataset_name);
    int get_n_dataset() {
        return var_table.size();
    };
    string get_dataset_name (int i_dataset) {
        return get_msg_key (i_dataset, GribFileData::var_name_id);
    };
    string get_dataset_long_name (int i_dataset) {
        return get_msg_key (i_dataset, GribFileData::parameter_name_id);
    };
    string get_dataset_level (int i_dataset) {
        return get_msg_key (i_dataset, GribFileData::level_id);
    };
    void print_msg_keys (int i_msg, unsigned int keys_filter = GRIB_KEYS_ITERATOR_ALL_KEYS);
    void get_msg_keys (int i_msg, vector < pair < string, string > > & v_keys,
        unsigned int keys_filter = GRIB_KEYS_ITERATOR_ALL_KEYS);
    void print_keys(unsigned int keys_filter = GRIB_KEYS_ITERATOR_ALL_KEYS);
    int get_dataset_data_type(string dataset_name);
    void print_param_table (void);
//     {return this->get_hdf_file()->get_dataset_data_type(dataset_name.c_str());};

// string get_values_attr_dataset(string dataset_name,string attr_name);

//     bool has_attr_dataset(string dataset_name,string attr_name);
/*
    string get_values_attr(string attr_name);

    string get_file_attr(string attr_name);

    bool has_attr(string attr_name);

    bool has_file_attr(string attr_name);*/
    /************************
    ** DATA ACCESS METHODS **
    ************************/

    void * read_data (void * data, const char * dataset_name, int * start = NULL, int * stride = NULL, int * edges = NULL, int rank = -1);

    template <typename T>
    T * get_msg_data (const int i_msg, vector <int> & sz_dim, T * outdata = NULL);
    const void print_msg_data (const int i_msg);
    void init_selector (const char* ds_name, int *&start, int *&stride, int *&edges, int &rank, bool is_selector_allocated[3]);

    void free_selector(const bool *is_selector_allocated, int *start, int *stride, int *edges);
    void check_selector(const char* ds_name, int *start, int *stride, int *edges, const int rank);

//     void get_fillValue(const string &dataset_name, void *fillValue);

//     void get_scaling(const string &dataset_name, float64 &scale, float64 &offset) {
//         string key ("parameterName");
//         return get_msg_key (i_dataset, key);
//     };

    void get_scaling(int i_msg, double &scale, double &offset) {
        string s_scale = get_msg_key (i_msg, GribFileData::scale_factor_id);
        string s_offset = get_msg_key (i_msg, GribFileData::add_offset_id);
        scale = MyTools::to_num <double> (s_scale);
        offset = MyTools::to_num <double> (s_offset);
    };
    const int get_n_msg() const {
        return n_msg;
    };

private :
    void init_n_msg ();
    grib_handle * get_handle (int i_msg);
    void release_handle (grib_handle *handle);
    void init_param_table (void);
    static void print_var_desc ( const pair<string, t_grib_var_desc> & var_table_elt) {
        cout << "---" << var_table_elt.first << "---" << endl;
        cout << var_table_elt.second << endl;
    };

    static string parse_type_id (const unsigned int data_flag, const int data_width) {
        string type_id ("");
        // check if data is stored as float or integers : 3rd bit from right
//         printf ( "flag = %lx width = %ld\n", data_flag, data_width );
//         cout << (data_flag & 0x4) << " " << ((data_flag & 0x4) >> 3) << endl;
        bool is_int = ((data_flag & 0x4) >> 3) == 1;
        if (is_int)
            type_id = "int";
        else
            type_id = "float";
        // add data width
        type_id += MyTools::to_string(data_width);
        return type_id;
    };

};

template <typename T>
T * GribFileData::get_msg_data(const int i_msg, vector <int> & sz_dim, T * outdata) {
    // --- seek message number i_msg
    grib_handle * handle = get_handle (i_msg);
    if (handle == NULL) {
        string msg = "Fail to seek GRIB message number " + MyTools::to_string<int>(i_msg);
        g_exception e( __FILE__ , __LINE__ , msg );
        throw e;
    }
//     cout << get_dataset_name (i_msg) << " " << get_msg_key (i_msg, "level") << endl;
    // --- allocate inner buffer
    size_t sz_buf;
    GRIB_CHECK (grib_get_size (handle, GribFileData::data_id.c_str(), &sz_buf), 0);
    double * inner_buf = NULL;
    // if the output data are requested in double and outdata already allocated, a unique buffer can be used for read and returned data. Else a temporary double buffer is needed.
    if (outdata != NULL && typeid (T) == typeid(double)) {
        inner_buf = outdata;
    } else {
        inner_buf = new double [sz_buf];
        if (inner_buf == NULL) {
            g_exception e(__FILE__ , __LINE__ , "Unable to allocate buffer\n");
            throw e;
        }
    }
    // --- read data
    GRIB_CHECK(grib_get_double_array(handle, GribFileData::data_id.c_str(), inner_buf, &sz_buf), 0);
    // --- read output dimensions size
    sz_dim.resize (2);
    long long_int;
    GRIB_CHECK(grib_get_long(handle, GribFileData::x_dim_id.c_str(), &long_int),0);
    sz_dim[1] = long_int;
    GRIB_CHECK(grib_get_long(handle, GribFileData::y_dim_id.c_str(), &long_int),0);
    sz_dim[0] = long_int;
    release_handle(handle);

    // --- format returned output buffer
    if (typeid (T) == typeid(double)) {
        // double typed output is requested : just return the inner buffer
        return inner_buf;
    } else {
        // innner and returned buffer does not have the same type -> allocate returned one and cast data
        outdata = new T [sz_buf];
        if (outdata == NULL) {
            g_exception e(__FILE__ , __LINE__ , "Unable to allocate buffer\n");
            throw e;
        }
        copy (inner_buf, inner_buf + sz_buf, outdata);
        free (inner_buf);
    }
    return outdata;
}

#endif