/*************************************************************************** * 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 ECMWFFILEDATA_H #define ECMWFFILEDATA_H #include "gribfiledata.h" #include "meteofiledata.h" #include "interpolation.h" #include "earth_geometry.h" /** * @author Nicolas PASCAL * @class ECMWFFileData manage the reading of ECMWF GRIB files archived at ICARE * @TODO remove FileDataReader interface implementation */ class ECMWFFileData : public MeteoFileData, public GribFileData, public FileDataReader { protected: /** * @enum ECMWFProductType the different types of ECMWF products */ typedef enum { UNDEFINED = 0, ANALYSIS, FORECAST } ECMWFProductType; /** the type of ECMWF product */ ECMWFProductType product_id; /** the product version */ string version; /** number of grid cells along the longitudes */ int sz_x; /** number of grid cells along the latitudes */ int sz_y; /** latitude of the up grid point (border of cell) */ float lat0; /** longitude of the left grid point (border of cell)*/ float lon0; /** * @brief compute the time range covered by this file, depending of the product */ void set_time_coverage(); /** * @brief extract sone informations using the filename * @param short_filename the filename without its path */ void parse_filename(const string & short_filename); /** * @brief test if the filename seems to be a valid ECMWF file * @param short_filename the filename without its path */ bool check_filename(const string & short_filename) const; /** * @brief read the grid characteristics in the file attributes */ void init_grid (); /** * @brief if not already set, read the time levels in the file * Interface implementation purpose. GRIB ECMWF files are composed of only one time */ virtual void load_time_level() {;}; public: /** vector of isobaric pressures, in hPa */ static const unsigned short isobaric_level_pressures[]; /** maximum number of isobaric pressure levels */ static const int nb_max_isobaric_level; /** a, b coefficients used to compute the pressure of the ECMWF 91 vertical hybrid levels * Source : http://www.ecmwf.int/products/data/technical/model_levels/model_def_91.html */ static const double a_hybrid_91[92]; static const double b_hybrid_91[92]; /** number of hybrid levels */ static const int nb_level_hybrid_91; /** * @brief constructor * @param name the filename * @param mode opening mode. only "r" at this time */ ECMWFFileData(const string &name, const string &mode="r"); /** * @brief destructor */ ~ECMWFFileData(); /** * @brief free eventually loaded data. * Not used here, only for interface compatibility */ void free_geolocation_data() {;}; /** * @brief find the index of the nearest point to (lat,lon,time) in the data. * time is ignored for ECMWF files because only one time level is used * @throw g_exception if lat or lon are invalid * @param lat the latitude * @param lon the longitude * @param time the time (using TAI convention) * @param indexes the output indexes : a 3 values array {lat_index,lon_index,time_index} * @return true if the coincidence could have been found. False else */ inline const bool get_index(const float &lat, const float& lon, const double &time, int * indexes); /** * @brief find the index of the nearest point to (lat,lon) in the data. * If (lat,lon,time) is not contained in this file, returned indexes are [-1,-1] * @param lat the latitude * @param lon the longitude * @param indexes the output indexes : a 2 values array {lat_index,lon_index} * @return true if the coincidence could have been found. False else */ inline const bool get_index(const float &lat, const float& lon, int * indexes); /** * @brief find the index of the nearest point to (lat,lon, p_lvl) in the data. * @warning not implemented at this time * If (lat,lon,p_lvl) is not contained in this file, returned indexes are [-1,-1, -1] * @param lat the latitude * @param lon the longitude * @param p_lvl the isobaric pressure level in hPa * @param indexes the output indexes : a 2 values array {lat_index,lon_index} * @return true if the coincidence could have been found. False else */ inline const bool get_index(const float &lat, const float& lon, const float &p_lvl, int * indexes) { UnimplementedMethod e(__FILE__,__LINE__,__PRETTY_FUNCTION__); throw e; }; /** * @brief return the index of the grid point that has the nearest latitude to @a lat * @param lat the latitude * @return the index, or raise an exception if out of range */ inline const int get_lat_index(const float &lat) const; /** * @brief return the index of the grid point that has the nearest longitude to @a lon * @param lon the longitude, in [-180, 180] range * @return the index, or raise an exception if out of range */ inline const int get_lon_index(const float &lon) const; /** * @brief return the latitude of the grid point at indice @a i_y * @param i_y grid point indice along Y * @return the grid point latitude */ inline const float32 get_lat (const int &i_y) const; /** * @brief return the longitude of the grid point at indice @a i_x * @param i_x grid point indice along X * @return the grid point longitude */ inline const float32 get_lon (const int &i_x) const; /** * @brief read a subset of a data set * @param data the data output buffer * @param sds_name the data set to read * @param start the start point * @param stride the step between to read values * @param edges the number of values to read in each direction * @param rank the rank of the sds * @return a pointer to the data buffer. Use it if the allocation is done inside this method */ void * read_data(void * data, const char * sds_name, int * start=NULL, int * stride=NULL, int *edges=NULL, int rank=-1) { return this->GribFileData::read_data(data, sds_name, start, stride, edges, rank); }; /** * @brief accessor to the product type ID * @return the product type ID */ ECMWFProductType get_product_id() const { return product_id; }; /** * @brief accessor to the product type as a string * @return the product type as a string */ const string get_s_product_id () const { string s_prod_id (""); switch (product_id) { case UNDEFINED: s_prod_id = "UNDEFINED"; break; case ANALYSIS: s_prod_id = "ANALYSIS"; break; case FORECAST: s_prod_id = "FORECAST"; break; default : g_exception e(__FILE__ , __LINE__ , "Invalid product type ID " + s_prod_id); throw e; }; return s_prod_id; }; /** *@brief accessor to the number of grid cells along the longitudes *@return the number of grid cells along the longitudes */ int get_sz_x() const { return sz_x; } /** *@brief accessor to the number of grid cells along the latitudes *@return the number of grid cells along the latitudes */ int get_sz_y() const { return sz_y; } /** * @brief return a description of the grid (resolution, number of cells,...) as a string * @return a description of the grid */ const string get_grid_desc (); /** * @brief returns the pressure levels of the given variable * @param var_name the name of the searched variable * @return the pressure levels */ const vector get_pressure_levels(const string & var_name); /** * @brief returns the number of vertical levels of the given variable * @param var_name the name of the searched variable * @return the number of pressure levels */ const int get_nb_pressure_levels (const string & var_name); /** * @brief accessor to the pressures of the 91 model levels in hPa * Source : http://www.ecmwf.int/research/ifsdocs/DYNAMICS/Chap2_Discretization4.html#961180 * @param p_surf surface pressure in hPa * @param v_p_hybrid [OUT] vector of hybrid pressure full levels, from top to bottom. * @return the pressures of the 91 model levels */ void get_hybrid_level_pressures(const double &p_surf, double * v_p_hybrid) const { for (int i = 0 ; i < nb_level_hybrid_91 ; ++i) v_p_hybrid[i] = get_hybrid_level_pressure (i, p_surf); }; /** * @brief compute the pressure of an hybrid level, in the 91 levels model * Source : http://www.ecmwf.int/research/ifsdocs/DYNAMICS/Chap2_Discretization4.html#961180 * @param i_level number of the level. Must ne in range [0,90] * @param p_surf surface pressure in hPa * @return the pressure of the hybrid level */ const double get_hybrid_level_pressure (const int i_level, const double &p_surf) const { if (i_level < 0 || i_level > nb_level_hybrid_91 ) { g_exception e(__FILE__ , __LINE__ , "Invalid hybrid level number " + MyTools::to_string(i_level) + " Must in range [0, 91]"); throw e; } // cout << "i_level=" << i_level << " a=[" << a_hybrid_91 [i_level] << "," << a_hybrid_91 [i_level+1] << "] b=[" << b_hybrid_91 [i_level] << "," << b_hybrid_91 [i_level+1] << "] psurf=" << p_surf << " -> " << .5 * ((a_hybrid_91 [i_level] + a_hybrid_91 [i_level+1]) * 0.01 + (b_hybrid_91 [i_level] + b_hybrid_91 [i_level+1]) * p_surf) << endl; return .5 * ((a_hybrid_91 [i_level] + a_hybrid_91 [i_level+1]) * 0.01 + (b_hybrid_91 [i_level] + b_hybrid_91 [i_level+1]) * p_surf); }; /** * @brief access to the isobaric levels pressures vector, in hPa and sorted in decreasing order * @return the isobaric levels pressures vector */ static const unsigned short* get_isobaric_level_pressures() { return isobaric_level_pressures; } /** * @brief compute the pressure at a list of altitudes * @param v_out_pres [OUT] output pressure in hPa * @param v_alt [IN] altitudes where to compute the pressure in km * @param sz_alt [IN] number of altitudes and output pressures * @param v_geop [IN] geopotential on input levels in m^2.s^-2 * @param v_inpres [IN] pressure on input levels in hPa. Stricly decreasing * @param n_pres [IN] number of input levels * @param interp_type [IN] type of interpolation to use : "linear" , "linear_log" or "spline" */ template void get_v_pres_on_alt ( Tyout * v_out_pres, const Txout * v_alt, const int sz_alt, const Txin * v_geop, const Tyin * v_inpres, const int n_pres, const string interp_type = "linear"); /** * @brief compute the pressure at a given altitude * @warning no extrapolation done. If @a alt is out of the latitudes represented in @a v_ingeop NaN is returned * @param alt [IN] altitude where to compute the pressure in km * @param v_ingeop [IN] geopotential on input levels in m^2.s^-2 * @param v_inpres [IN] pressure on input levels in hPa. Stricly decreasing * @param sz_inpres [IN] number of input levels * @param interp_type [IN] type of interpolation to use : "linear" , "linear_log" or "spline" * @return the output pressure in hPa */ template double get_pres_on_alt ( const Txout alt, const Txin * v_ingeop, const Tyin * v_inpres, const int sz_inpres, const string interp_type = "linear"); /** * @brief compute the geopotential of an altitude in m * The formula used is : geop = g0 * alt. g0 is assumed to be constant in the atmosphere and on the surface on the earth * @param alt an altitude in m * @return the geopotential of this altitude, in m2.s-2 */ double alt_to_geop (const double & alt) { return alt * Earth::g0; }; /** * @brief compute the abscissa of the given pressures in the isobaric levels system. The abscissa is designed to be used lineary with log(P) * @param v_out_abs [OUT] pressures abscissa * @param v_out_pres [IN] searched pressures * @param sz_out_pres [IN] number of searched pressures */ template void get_abscissa_iso (Tyout * v_out_abs, const Txout * v_out_pres, const int sz_out_pres); /** * @brief compute the abscissa of the given pressures in the hybrid levels system. The abscissa is designed to be used lineary with log(P) * @param v_out_abs [OUT] pressures abscissa * @param v_out_pres [IN] searched pressures * @param sz_out_pres [IN] number of searched pressures * @param v_in_pres [IN] hybrid levels pressures * @param sz_in_pres [IN] number of hybrid levels pressures */ template void get_abscissa_hybrid (Tyout * v_out_abs, const Txout * v_out_pres, const int sz_out_pres, const Txin * v_in_pres, const int sz_in_pres); /** * @brief interpolate the variable on the given pressures, using precomputed abscissa based on a linear approach on log(P) * @param v_out_var [OUT] output variable interpolated on abscissa * @param v_out_abs [IN] abscissa of output level, in the reference coordinates system * @param sz_out [IN] number of elements to interpolate * @param v_in_var [IN] value of variable at reference levels * @param v_in_pres [IN] pressure of reference levels * @param sz_in [IN] number of reference levels */ template void interp_on_pres (Tyout * v_out_var, const Txout * v_out_abs, const int sz_out, const Tyin * v_in_var, const Txin * v_in_pres, const int sz_in); /** --- FOR OLD DATA2GRID INTERFACE compatibility - START --- */ void close_data_file(); void open_data_file(); // vector get_dataset_dimension(const string &sds_name); // void get_dataset_fill_value(const string &name,void *fillvalue); int get_n_dataset(); string get_dataset_name(int i); int get_dataset_data_type(string sds_name); string get_values_attr_dataset(string sds_name,string attr); bool has_attr_dataset(string sds_name,string attr); /** --- FOR OLD DATA2GRID INTERFACE compatibility - END --- */ }; const int ECMWFFileData::get_lat_index(const float & lat) const { if (lat < -90. || lat > 90.) { g_exception e(__FILE__ , __LINE__ , "Invalid latitude " + MyTools::to_string(lat)); throw e; } return MyTools::round(fabs(lat - lat0) / delta_lat); } const float32 ECMWFFileData::get_lat(const int & i_y) const { return lat0 - i_y * delta_lat; } const float32 ECMWFFileData::get_lon(const int & i_x) const { float32 lon = lon0 + i_x * delta_lon; if (lon > 180.) lon -= 360.; return lon; } const int ECMWFFileData::get_lon_index(const float & _lon) const { if (_lon < -180. || _lon > 180.) { g_exception e(__FILE__ , __LINE__ , "Invalid longitude " + MyTools::to_string(_lon)); throw e; } // longitudes as [-180, 180] -> [0, 360[ float lon (_lon); if (lon < 0.) lon = 360. + lon; int ilon = MyTools::round((lon - lon0) / delta_lon); // treat extreme case of longitudes in range [-delta_lon/2, 0] if (ilon == sz_x) ilon = 0; return ilon; } const bool ECMWFFileData::get_index(const float & lat, const float & lon, int * indexes) { indexes [0] = get_lat_index(lat); indexes [1] = get_lon_index(lon); return true; } const bool ECMWFFileData::get_index(const float & lat, const float & lon, const double & time, int * indexes) { return get_index (lat, lon, indexes); } template void ECMWFFileData::get_v_pres_on_alt( Tyout * v_outpres, const Txout * v_alt, const int sz_alt, const Txin * v_ingeop, const Tyin * v_inpres, const int sz_inpres, const string interp_type) { // extrapolation forbidden -> NaN returned if out of range double nan = numeric_limits::signaling_NaN(); if (interp_type == "linear" || interp_type == "linear_log") { for (int i = 0 ; i < sz_alt ; ++i) v_outpres [i] = get_pres_on_alt (v_alt[i], v_ingeop, v_inpres, sz_inpres, interp_type) ; } else if (interp_type == "spline") { // computations made in double. Last geopotential is ignored because set on a virtual pressure of 0 hPa. This is the reason of the "- 1" in size of the input pressures vector int sz_in = sz_inpres - 1; double x_in [sz_in]; copy (v_ingeop, v_ingeop + sz_in, x_in); double y_in [sz_in]; copy (v_inpres, v_inpres + sz_in, y_in); int sz_out = sz_alt; double x_out [sz_out]; for (int i = 0 ; i < sz_out ; ++i) x_out[i] = alt_to_geop (v_alt[i] * 1000.); // * 1000 : km -> m; double y_out [sz_out]; // compute spline interp // ------------------------------------------------------- spline_interp (y_out, x_out, sz_out, y_in, x_in, sz_in, nan); copy (y_out, y_out + sz_out, v_outpres); } else { g_exception e(__FILE__ , __LINE__ , "Invalid interpolation type " + interp_type); throw e; } } template double ECMWFFileData::get_pres_on_alt ( const Txout alt, const Txin * v_ingeop, const Tyin * v_inpres, const int sz_inpres, const string interp_type) { // geopotential of searched altitude double outgeop = alt_to_geop (alt * 1000.); // * 1000 : km -> m; // extrapolation forbidden -> NaN returned if out of range double nan = numeric_limits::signaling_NaN(); // last geopotential is ignored because set on a virtual pressure of 0 hPa -> -2 if (outgeop < v_ingeop[0] || outgeop > v_ingeop [sz_inpres - 2]) return nan; // process to linear interpolation of pressure int i = 0; while ((i < (sz_inpres - 2)) && (outgeop > v_ingeop[i])) ++i; // if (i >= (sz_inpres - 2)) if (i >= (sz_inpres - 1)) return nan; else { double outpres = nan; if (interp_type == "linear") { outpres = v_inpres[i - 1] + (outgeop - v_ingeop[i - 1]) * (v_inpres[i] - v_inpres[i - 1])/(v_ingeop[i] - v_ingeop[i - 1]); } else if (interp_type == "linear_log") { outpres = log(v_inpres[i - 1]) + (log(outgeop) - log(v_ingeop[i - 1])) * (log(v_inpres[i]) - log(v_inpres[i - 1]))/(log(v_ingeop[i]) - log(v_ingeop[i - 1])); outpres = exp(outpres); } else { g_exception e(__FILE__ , __LINE__ , "Invalid interpolation type " + interp_type); throw e; } if (outpres < 0) outpres = nan; // cout << alt << " " << outgeop << " v_ingeop=[" << v_ingeop[i - 1] << "," << v_ingeop[i] << "] v_inpres=[" << v_inpres[i - 1] << "," << v_inpres[i] << "] -> " << outpres << endl; return outpres; } } template void ECMWFFileData::get_abscissa_iso(Tyout * v_out_abs, const Txout * v_out_pres, const int sz_out) { MyTools::get_log_abscissa_sequential (v_out_abs, v_out_pres, sz_out, isobaric_level_pressures, nb_max_isobaric_level - 1); } template void ECMWFFileData::get_abscissa_hybrid (Tyout * v_out_abs, const Txout * v_out_pres, const int sz_out, const Txin * v_in_pres, const int sz_in) { MyTools::get_log_abscissa_sequential (v_out_abs, v_out_pres, sz_out, v_in_pres, sz_in); } template void ECMWFFileData::interp_on_pres(Tyout * v_out_var, const Txout * v_out_abs, const int sz_out, const Tyin * v_in_var, const Txin * v_in_pres, const int sz_in) { // extrapolation forbidden -> NaN returned if out of range Tyout nan = numeric_limits::signaling_NaN(); int i1, i2; double y1, y2; double abs; for (int i = 0 ; i < sz_out ; ++i) { abs = v_out_abs [i]; if (isnan(abs)) v_out_var [i] = nan; else { i1 = long(floor(abs)); i2 = i1 + 1; y1 = v_in_var [i1]; y2 = v_in_var [i2]; v_out_var [i] = y1 + (abs - i1) * (y2 - y1); } } } #endif