WSL/SLF GitLab Repository

Commit 7b639ecc authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A new date format has been added: RFC868. Two MeteoGrids have been renamed:...

A new date format has been added: RFC868. Two MeteoGrids have been renamed: ISW_DIR as ISWR_DIR and ISW_DIFF as ISWR_DIFF while a new grid has been declared (P_SEA). A few new constants have been declared. Some new links into the GRIB documentation as well as support for TAU_CLD. The NetCDF support is much improved, although there are still things left to do (writing grids out has been disabled for the moment, the parsing of the time units is currently not done, there is no support for automatically getting the data out of multiple files in a directory, etc).
parent d7fb729c
......@@ -49,6 +49,7 @@ const int Date::daysNonLeapYear[12] = {31,28,31,30,31,30,31,31,30,31,30,31};
const double Date::DST_shift = 1.0; //in hours
const float Date::MJD_offset = 2400000.5; ///<offset between julian date and modified julian date
const float Date::Unix_offset = 2440587.5; ///<offset between julian date and Unix Epoch time
const float Date::RFC868_offset = 2415020.5; ///< offset between julian date and RFC868 time (ref is 1900-01-01T00:00 GMT)
const float Date::Excel_offset = 2415018.5; ///<offset between julian date and Excel dates (note that excel invented some days...)
const float Date::Matlab_offset = 1721058.5; ///<offset between julian date and Matlab dates
......@@ -272,6 +273,17 @@ void Date::setModifiedJulianDate(const double& julian_in, const double& i_timezo
setDate(tmp_julian, i_timezone, i_dst);
}
/**
* @brief Set date from an RFC868 date (time since 1900-01-01T00:00 GMT).
* @param julian_in julian date to set
* @param i_timezone timezone as an offset to GMT (in hours, optional)
* @param i_dst is it DST? (default: no)
*/
void Date::setRFC868Date(const double& julian_in, const double& i_timezone, const bool& i_dst) {
const double tmp_julian = julian_in + RFC868_offset;
setDate(tmp_julian, i_timezone, i_dst);
}
/**
* @brief Set date from a Unix date.
* @param in_time unix time (ie: as number of seconds since Unix Epoch, always UTC)
......@@ -371,6 +383,24 @@ double Date::getModifiedJulianDate(const bool& gmt) const {
}
}
/**
* @brief Return RFC868 date.
* The RFC868 date is defined as the fractional number of days since 1900-01-01T00:00 GMT
* @param gmt convert returned value to GMT? (default: false)
* @return RFC868 julian date in the current timezone / in GMT depending on the gmt parameter
*/
double Date::getRFC868Date(const bool& gmt) const {
if(undef==true)
throw UnknownValueException("Date object is undefined!", AT);
if(gmt) {
return (gmt_julian - RFC868_offset);
} else {
const double local_julian = GMTToLocal(gmt_julian);
return (local_julian - RFC868_offset);
}
}
/**
* @brief Return truncated julian date (TJD).
* The truncated julian date is defined as the julian day shifted to start at 00:00 and modulo 10000 days.
......
......@@ -78,6 +78,7 @@ class Date {
static const int daysNonLeapYear[];
static const double DST_shift;
static const float MJD_offset;
static const float RFC868_offset;
static const float Unix_offset;
static const float Excel_offset;
static const float Matlab_offset;
......@@ -99,6 +100,7 @@ class Date {
void setDate(const int& year, const unsigned int& month, const unsigned int& day, const unsigned int& hour, const unsigned int& minute, const double& second, const double& in_timezone, const bool& in_dst=false);
void setDate(const time_t& in_time, const bool& in_dst=false);
void setModifiedJulianDate(const double& julian_in, const double& in_timezone, const bool& in_dst=false);
void setRFC868Date(const double& julian_in, const double& in_timezone, const bool& in_dst=false);
void setUnixDate(const time_t& in_time, const bool& in_dst=false);
void setExcelDate(const double excel_in, const double& in_timezone, const bool& in_dst=false);
void setMatlabDate(const double excel_in, const double& in_timezone, const bool& in_dst=false);
......@@ -110,6 +112,7 @@ class Date {
double getJulian(const bool& gmt=false) const;
double getModifiedJulianDate(const bool& gmt=false) const;
double getTruncatedJulianDate(const bool& gmt=false) const;
double getRFC868Date(const bool& gmt=false) const;
time_t getUnixDate() const;
double getExcelDate(const bool& gmt=false) const;
double getMatlabDate(const bool& gmt=false) const;
......
......@@ -43,8 +43,8 @@ bool MeteoGrids::initStaticData()
paramname.push_back("VW_MAX");
paramname.push_back("ISWR");
paramname.push_back("RSWR");
paramname.push_back("ISW_DIFF");
paramname.push_back("ISW_DIR");
paramname.push_back("ISWR_DIFF");
paramname.push_back("ISWR_DIR");
paramname.push_back("ILWR");
paramname.push_back("TAU_CLD");
paramname.push_back("HS");
......@@ -54,6 +54,7 @@ bool MeteoGrids::initStaticData()
paramname.push_back("TSG");
paramname.push_back("TSS");
paramname.push_back("P");
paramname.push_back("P_SEA");
paramname.push_back("U");
paramname.push_back("V");
paramname.push_back("W");
......
......@@ -55,8 +55,8 @@ class MeteoGrids {
VW_MAX, ///< Maximum wind velocity
ISWR, ///< Incoming short wave radiation
RSWR, ///< Reflected short wave radiation
ISW_DIFF, ///< Incoming short wave, diffuse
ISW_DIR, ///< Incoming short wave, direct
ISWR_DIFF, ///< Incoming short wave, diffuse
ISWR_DIR, ///< Incoming short wave, direct
ILWR, ///< Incoming long wave radiation
TAU_CLD, ///< Cloud transmissivity or ISWR/ISWR_clear_sky
HS, ///< Height of snow
......@@ -66,6 +66,7 @@ class MeteoGrids {
TSG, ///< Temperature ground surface
TSS, ///< Temperature snow surface
P, ///< Air pressure
P_SEA, ///< Sea level air pressure
U, ///< East component of wind
V, ///< North component of wind
W, ///< Vertical component of wind
......
......@@ -26,6 +26,7 @@ namespace Cst {
const double stefan_boltzmann = 5.670373e-8; // (W m-2 K-4)
const double gravity = 9.80665; // (m s-2)
const double std_press = 101325.; // (Pa)
const double dry_adiabatique_lapse_rate = 0.0065; // (K/m)
const double gaz_constant_dry_air = 287.058; // (J kg-1 K-1)
const double gaz_constant_water_vapor = 461.9; // (J kg-1 K-1)
......@@ -38,6 +39,7 @@ namespace Cst {
const double l_water_vaporization = 2.504e6; // (J Kg-1)
const double l_water_fusion = 3.34e5; // (J Kg-1)
const double water_molecular_mass = 18.0153e-3; // (Kg)
const double dry_air_mol_mass =0.0289644 ; // (Kg/mol)
const double specific_heat_ice = 2100.0; // (J K-1), at 0C
const double specific_heat_water = 4190.0; // (J K-1) at 0C
......
......@@ -37,7 +37,7 @@ namespace mio {
* This plugin reads GRIB (https://en.wikipedia.org/wiki/GRIB) files as produced by meteorological models.
* Being based on GRIB API (https://software.ecmwf.int/wiki/display/GRIB/Home), it should support both version 1 and 2 of the format (please note that grib_api must be compiled with Position Independent Code ("fPIC" flag)).
* Fields are read based on their marsParam code (this is built as {grib parameter number}.{grib table number} the table being preferably table 2, the parameter being preferably WMO standardized, as in http://rda.ucar.edu/docs/formats/grib/gribdoc/params.html) and levels
* (levels description is available at http://www.nco.ncep.noaa.gov/pmb/docs/on388/). Standard COSMO grids are listed at http://zephyr.ucd.ie/mediawiki/index.php/COSMO_GRIB .
* (levels description is available at http://www.nco.ncep.noaa.gov/pmb/docs/on388/). Standard COSMO grids are listed at http://zephyr.ucd.ie/mediawiki/index.php/COSMO_GRIB while NCEP are listed at http://www.cpc.ncep.noaa.gov/products/wesley/opn_gribtable.html.
*
* Several assumptions/approximations are held/made when reading grids:
* - since models usually use rotated latitude/longitude (see http://www.cosmo-model.org/content/model/documentation/core/default.htm , part I, chapter 3.3), the center of the domain can be approximated by a tangential cartesian coordinate system. We therefore don't re-project the lat/lon grid and use it "as is".
......@@ -578,10 +578,10 @@ void GRIBIO::read2DGrid(const std::string& filename, Grid2DObject& grid_out, con
grid_out.grid2D *= -1.;
} else read2DGrid_indexed(25.201, 1, 0, date, grid_out); //ALWD_S
}
/*if(parameter==MeteoGrids::CLD) { //cloudiness
if(parameter==MeteoGrids::TAU_CLD) { //cloudiness
if(read2DGrid_indexed(74.2, 1, 0, date, grid_out)) //CLCM
grid_out.grid2D /= 100.;
}*/
}
if(parameter==MeteoGrids::ISWR) {
if(read2DGrid_indexed(116.2, 1, 0, date, grid_out)) { //short wave
......
......@@ -18,6 +18,8 @@
#include "NetCDFIO.h"
#include <meteoio/ResamplingAlgorithms2D.h>
#include <meteoio/meteoStats/libinterpol1D.h>
#include <meteoio/meteoLaws/Meteoconst.h>
#include <meteoio/meteoLaws/Atmosphere.h>
#include <meteoio/Timer.h>
#include <meteoio/MathOptim.h>
#include <meteoio/plugins/libncpp.h>
......@@ -80,14 +82,16 @@ const std::string NetCDFIO::cf_longitude = "lon";
const std::string NetCDFIO::cf_altitude = "z";
NetCDFIO::NetCDFIO(const std::string& configfile) : cfg(configfile), in_attributes(), out_attributes(), coordin(), coordinparam(), coordout(), coordoutparam(),
in_dflt_TZ(0.), out_dflt_TZ(0.), in_strict(false), out_strict(false), vecMetaData()
in_dflt_TZ(0.), out_dflt_TZ(0.), in_time_offset(0.), in_time_multiplier(1.), out_time_offset(0.), out_time_multiplier(1.),
in_strict(false), out_strict(false), vecMetaData()
{
IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
parseInputOutputSection();
}
NetCDFIO::NetCDFIO(const Config& cfgreader) : cfg(cfgreader), in_attributes(), out_attributes(), coordin(), coordinparam(), coordout(), coordoutparam(),
in_dflt_TZ(0.), out_dflt_TZ(0.), in_strict(false), out_strict(false), vecMetaData()
in_dflt_TZ(0.), out_dflt_TZ(0.), in_time_offset(0.), in_time_multiplier(1.), out_time_offset(0.), out_time_multiplier(1.),
in_strict(false), out_strict(false), vecMetaData()
{
IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
parseInputOutputSection();
......@@ -102,19 +106,27 @@ void NetCDFIO::parseInputOutputSection()
cfg.getValue("TIME_ZONE", "Input", in_dflt_TZ, IOUtils::nothrow);
cfg.getValue("TIME_ZONE", "Output", out_dflt_TZ, IOUtils::nothrow);
initAttributesMap(cfg.get("NETCDF_SCHEMA", "Input", IOUtils::nothrow), in_attributes);
initAttributesMap(cfg.get("NETCDF_SCHEMA", "Output", IOUtils::nothrow), out_attributes);
/*std::cout << "In in_attributes map:\n";
for (std::map<MeteoGrids::Parameters, NetCDFIO::attributes>::const_iterator it = in_attributes.begin(); it != in_attributes.end(); ++it){
std::cout << MeteoGrids::getParameterName((*it).first) << " = { " << (*it).second.var << " , " << (*it).second.standard_name << " , " << (*it).second.long_name << " , " << (*it).second.units << "}\n";
const string in_schema = IOUtils::strToUpper( cfg.get("NETCDF_SCHEMA", "Input", IOUtils::nothrow) );
initAttributesMap(in_schema, in_attributes);
setTimeTransform(in_schema, in_time_offset, in_time_multiplier);
const string out_schema = IOUtils::strToUpper( cfg.get("NETCDF_SCHEMA", "Output", IOUtils::nothrow) );
initAttributesMap(out_schema, out_attributes);
setTimeTransform(out_schema, out_time_offset, out_time_multiplier);
}
void NetCDFIO::setTimeTransform(const std::string& schema, double &time_offset, double &time_multiplier)
{
if (schema=="CNRM") {
time_offset = Date::MJD_offset;
} else if (schema=="ECMWF") {
time_offset = Date::RFC868_offset;
time_multiplier = 1./24.;
}
exit(0);*/
}
void NetCDFIO::initAttributesMap(std::string schema, std::map<MeteoGrids::Parameters, NetCDFIO::attributes> &attr)
void NetCDFIO::initAttributesMap(const std::string& schema, std::map<MeteoGrids::Parameters, attributes> &attr)
{
if (schema.empty()) return;
IOUtils::toUpper(schema);
if (schema=="CF1") {
attr[MeteoGrids::DEM] = attributes("z", "altitude", "height above mean sea level", "m", IOUtils::nodata);
......@@ -127,24 +139,28 @@ void NetCDFIO::initAttributesMap(std::string schema, std::map<MeteoGrids::Parame
attr[MeteoGrids::AZI] = attributes("aspect", "", "slope aspect", "degrees from north", IOUtils::nodata);
attr[MeteoGrids::TA] = attributes("Tair", "", "Near Surface Air Temperature", "K", IOUtils::nodata);
attr[MeteoGrids::RH] = attributes("HUMREL", "", "Relative Humidity", "%", IOUtils::nodata);
attr[MeteoGrids::QI] = attributes("Qair", "", "", "", IOUtils::nodata);
attr[MeteoGrids::VW] = attributes("Wind", "", "Wind Speed", "m/s", IOUtils::nodata);
attr[MeteoGrids::DW] = attributes("Wind_DIR", "", "Wind Direction", "deg", IOUtils::nodata);
attr[MeteoGrids::QI] = attributes("Qair", "", "", "", IOUtils::nodata);
attr[MeteoGrids::HNW_L] = attributes("Rainf", "", "Rainfall Rate", "kg/m2/s", IOUtils::nodata);
attr[MeteoGrids::HNW_S] = attributes("Snowf", "", "", "", IOUtils::nodata);
attr[MeteoGrids::ISW_DIR] = attributes("DIR_SWdown", "", "Surface Incident Direct Shortwave Radiation", "W/m2", IOUtils::nodata);
attr[MeteoGrids::ISW_DIFF] = attributes("SCA_SWdown", "", "", "", IOUtils::nodata);
attr[MeteoGrids::ISWR_DIR] = attributes("DIR_SWdown", "", "Surface Incident Direct Shortwave Radiation", "W/m2", IOUtils::nodata);
attr[MeteoGrids::ISWR_DIFF] = attributes("SCA_SWdown", "", "", "", IOUtils::nodata);
attr[MeteoGrids::P] = attributes("PSurf", "", "Surface Pressure", "Pa", IOUtils::nodata);
attr[MeteoGrids::ILWR] = attributes("LWdown", "", "Surface Incident Longwave Radiation", "W/m2", IOUtils::nodata);
} else if (schema=="ECMWF") {
attr[MeteoGrids::TA] = attributes("t2m", "", "", "", 2.);
attr[MeteoGrids::P] = attributes("sp", "", "", "", IOUtils::nodata);
attr[MeteoGrids::ISWR] = attributes("ssrd", "", "", "", IOUtils::nodata);
attr[MeteoGrids::ILWR] = attributes("strd", "", "", "", IOUtils::nodata);
attr[MeteoGrids::HNW] = attributes("tp", "", "", "", IOUtils::nodata);
attr[MeteoGrids::TD] = attributes("d2m", "", "", "", 2.);
attr[MeteoGrids::U] = attributes("u10m", "", "", "", 10.);
attr[MeteoGrids::V] = attributes("v10m", "", "", "", 10.);
attr[MeteoGrids::DEM] = attributes("z", "", "Geopotential", "m**2 s**-2", IOUtils::nodata);
attr[MeteoGrids::TA] = attributes("t2m", "", "2 metre temperature", "K", 2.);
attr[MeteoGrids::TD] = attributes("d2m", "", "2 metre dewpoint temperature", "K", 2.);
attr[MeteoGrids::P] = attributes("sp", "surface_air_pressure", "Surface pressure", "Pa", IOUtils::nodata);
attr[MeteoGrids::P_SEA] = attributes("msl", "air_pressure_at_sea_level", "Mean sea level pressure", "Pa", IOUtils::nodata);
attr[MeteoGrids::ISWR] = attributes("ssrd", "surface_downwelling_shortwave_flux_in_air", "Surface solar radiation downwards", "J m**-2", IOUtils::nodata);
attr[MeteoGrids::ILWR] = attributes("strd", "", "Surface thermal radiation downwards", "J m**-2", IOUtils::nodata);
attr[MeteoGrids::HNW] = attributes("tp", "", "Total precipitation", "m", IOUtils::nodata);
attr[MeteoGrids::U] = attributes("u10", "", "10 metre U wind component", "m s**-1", 10.);
attr[MeteoGrids::V] = attributes("v10", "", "10 metre V wind component", "m s**-1", 10.);
attr[MeteoGrids::SWE] = attributes("sd", "lwe_thickness_of_surface_snow_amount", "Snow depth", "m of water equivalent", IOUtils::nodata);
attr[MeteoGrids::TSS] = attributes("skt", "", "Skin temperature", "K", IOUtils::nodata);
} else
throw InvalidArgumentException("Invalid schema selected for NetCDF: \""+schema+"\"", AT);
}
......@@ -161,29 +177,102 @@ void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const std::string& arguments)
}
}
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& /*parameter*/, const Date& date)
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
{
const string filename = cfg.get("GRID2DFILE", "Input");
grid_out.clear();
const string filename = cfg.get("GRID2DFILE", "Input"); //HACK: also allow using GRID2DPATH
string varname="";
read2DGrid_internal(grid_out, filename, varname, date);
if (read2DGrid_internal(grid_out, filename, parameter, date)) return; //schema naming
if (parameter==MeteoGrids::VW || parameter==MeteoGrids::DW) { //VW, DW
Grid2DObject U,V;
const bool hasU = read2DGrid_internal(U, filename, MeteoGrids::U, date);
const bool hasV = read2DGrid_internal(V, filename, MeteoGrids::V, date);
if (hasU==true && hasV==true) {
grid_out.set(U, IOUtils::nodata);
if (parameter==MeteoGrids::VW) {
for (size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
grid_out(ii) = sqrt( Optim::pow2(U(ii)) + Optim::pow2(V(ii)) );
} else {
for (size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
grid_out(ii) = fmod( atan2( U(ii), V(ii) ) * Cst::to_deg + 360., 360.); // turn into degrees [0;360)
}
return;
}
}
if (parameter==MeteoGrids::RH) { //RH
Grid2DObject ta;
DEMObject dem;
bool hasDEM = false;
try {
readDEM(dem);
hasDEM = true;
} catch(...){}
const bool hasQI = read2DGrid_internal(grid_out, filename, MeteoGrids::QI, date);
const bool hasTA = read2DGrid_internal(ta, filename, MeteoGrids::TA, date);
if (hasQI && hasDEM && hasTA) {
for(size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
grid_out(ii) = Atmosphere::specToRelHumidity(dem(ii), ta(ii), grid_out(ii));
return;
}
const bool hasTD = read2DGrid_internal(grid_out, filename, MeteoGrids::TD, date);
if (hasTA && hasTD) {
for(size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
grid_out(ii) = Atmosphere::DewPointtoRh(grid_out(ii), ta(ii), false);
return;
}
}
if (parameter==MeteoGrids::ISWR) { //ISWR
Grid2DObject iswr_diff;
const bool hasISWR_DIFF = read2DGrid_internal(iswr_diff, filename, MeteoGrids::ISWR_DIFF);
const bool hasISWR_DIR = read2DGrid_internal(grid_out, filename, MeteoGrids::ISWR_DIR);
if (hasISWR_DIFF && hasISWR_DIR) {
grid_out += iswr_diff;
return;
}
}
if (parameter==MeteoGrids::HNW) { //HNW
Grid2DObject hnw_s;
const bool hasHNW_S = read2DGrid_internal(hnw_s, filename, MeteoGrids::HNW_S);
const bool hasHNW_L = read2DGrid_internal(grid_out, filename, MeteoGrids::HNW_L);
if (hasHNW_S && hasHNW_L) {
grid_out += hnw_s;
return;
}
}
throw InvalidArgumentException("Parameter \'"+MeteoGrids::getParameterName(parameter)+"\' either not found in file \'"+filename+"\' or not found in current NetCDF schema", AT);
}
void NetCDFIO::readDEM(DEMObject& dem_out)
{
//HACK
const string filename = cfg.get("DEMFILE", "Input");
const string varname = cfg.get("DEMVAR", "Input", IOUtils::nothrow);
if (!varname.empty()) {
if (!read2DGrid_internal(dem_out, filename, varname))
throw InvalidArgumentException("Variable \'"+varname+"\' not found in file \'"+filename+"\'", AT);
} else {
const string dem_var = in_attributes[MeteoGrids::DEM].var;
if (!dem_var.empty() && read2DGrid_internal(dem_out, filename, dem_var)) return;
if (read2DGrid_internal(dem_out, filename, MeteoGrids::DEM)) return; //schema naming
if (read2DGrid_internal(dem_out, filename, "Band1")) return; //ASTER naming
if (read2DGrid_internal(dem_out, filename, "z")) return; //GDAL naming
//last chance: read from pressure grids
Grid2DObject p, ta, p_sea;
if (read2DGrid_internal(p_sea, filename, MeteoGrids::P_SEA, Date(2015,1,1,12,0,0)) &&
read2DGrid_internal(p, filename, MeteoGrids::P, Date(2015,1,1,12,0,0)) &&
read2DGrid_internal(ta, filename, MeteoGrids::TA, Date(2015,1,1,12,0,0))) {
dem_out.set(p, IOUtils::nodata);
const double ex = Cst::gravity*Cst::dry_air_mol_mass / (Cst::gaz_constant*Cst::dry_adiabatique_lapse_rate);
const double ex_inv = 1./ex;
for(size_t ii=0; ii<(dem_out.getNx()*dem_out.getNy()); ii++)
dem_out(ii) = 1./Cst::dry_adiabatique_lapse_rate * (pow(p(ii)/p_sea(ii), -ex_inv) - ta(ii));
return;
}
throw InvalidArgumentException("The variable containing the DEM could not be found. Please specify it using the DEMVAR key.", AT);
}
}
......@@ -226,6 +315,9 @@ void NetCDFIO::readPOI(std::vector<Coords>&)
void NetCDFIO::write2DGrid(const Grid2DObject& grid_in, const std::string& arguments)
{
//Nothing so far
throw IOException("Nothing implemented here", AT);
// arguments is a string of the format filname:varname
vector<string> vec_argument;
IOUtils::readLineToVec(arguments, vec_argument, ':');
......@@ -238,12 +330,36 @@ void NetCDFIO::write2DGrid(const Grid2DObject& grid_in, const std::string& argum
void NetCDFIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameters& parameter, const Date& date)
{
//Nothing so far
throw IOException("Nothing implemented here", AT);
const string filename = cfg.get("GRID2DFILE", "Output");
//const string varname = get_varname(parameter);
const string varname = "";
write2DGrid_internal(grid_in, filename, varname, date);
}
bool NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& filename, const MeteoGrids::Parameters& parameter, const Date& date)
{
const std::map<MeteoGrids::Parameters, attributes>::const_iterator it = in_attributes.find(parameter);
if (it==in_attributes.end()) return false;
const bool status = read2DGrid_internal(grid_out, filename, it->second.var, date);
if (status==true) { //correct for other exotic units (for example, geopotential instead of height)
const string units = it->second.units;
if (units=="m**2 s**-2") grid_out /= Cst::gravity;
if (units=="%") grid_out /= 100.;
if (units=="m" && parameter==MeteoGrids::HNW) grid_out *= 100.;
if (parameter==MeteoGrids::HNW) { //very low precip reset to zero
for (size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
if (grid_out(ii)<1e-3) grid_out(ii)=0.;
}
if (units=="J m**-2") grid_out /= (3600.*3.);
}
return status;
}
bool NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& filename, const std::string& varname, const Date& date)
{
int ncid, varid;
......@@ -256,7 +372,7 @@ bool NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& fi
ncpp::get_variable(ncid, varname, varid);
ncpp::get_dimension(ncid, varname, varid, dimid, dim_varid, dimname, dimlen);
const bool is_record = (date != Date()); //HACK: other possibility: if the file only contains 1 grid and no date -> is_record=false
const bool is_record = (date != Date());
size_t lat_index = 0, lon_index = 1;
if (is_record) { // In case we're reading a record the first index is always the record index
lat_index = 1;
......@@ -283,8 +399,20 @@ bool NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& fi
//read gridded data
double *grid = new double[dimlen[lat_index]*dimlen[lon_index]];
if (is_record) {
const size_t pos = ncpp::find_record(ncid, NetCDFIO::cf_time, dimid[0], date.getModifiedJulianDate());
if (pos == IOUtils::npos) throw IOException("No record for date " + date.toString(Date::ISO), AT);
const double timestamp = (date.getJulian() - in_time_offset) / in_time_multiplier;
const size_t pos = ncpp::find_record(ncid, NetCDFIO::cf_time, dimid[0], timestamp);
if (pos == IOUtils::npos) {
double min, max;
const bool status = ncpp::get_recordMinMax(ncid, NetCDFIO::cf_time, dimid[0], min, max);
if (status) {
Date d_min, d_max;
d_min.setDate(min*in_time_multiplier + in_time_offset, in_dflt_TZ);
d_max.setDate(max*in_time_multiplier + in_time_offset, in_dflt_TZ);
throw IOException("No record for date " + date.toString(Date::ISO) + ". Records avec between " + d_min.toString(Date::ISO) + " and " + d_max.toString(Date::ISO), AT);
}
else
throw IOException("No record for date " + date.toString(Date::ISO), AT);
}
ncpp::read_data(ncid, varname, varid, pos, dimlen[lat_index], dimlen[lon_index], grid);
} else {
ncpp::read_data(ncid, varname, varid, grid);
......
......@@ -71,15 +71,18 @@ class NetCDFIO : public IOInterface {
double height;
} attributes;
void initAttributesMap(std::string schema, std::map<MeteoGrids::Parameters, attributes> &attr);
void initAttributesMap(const std::string& schema, std::map<MeteoGrids::Parameters, attributes> &attr);
void setTimeTransform(const std::string& schema, double &time_offset, double &time_multiplier);
void parseInputOutputSection();
void check_consistency(const int& ncid, const Grid2DObject& grid, double*& lat_array, double*& lon_array,
int& did_lat, int& did_lon, int& vid_lat, int& vid_lon);
int& did_lat, int& did_lon, int& vid_lat, int& vid_lon);
bool read2DGrid_internal(Grid2DObject& grid_out, const std::string& filename, const MeteoGrids::Parameters& parameter, const Date& date=Date());
bool read2DGrid_internal(Grid2DObject& grid_out, const std::string& full_name, const std::string& varname, const Date& date=Date());
void write2DGrid_internal(const Grid2DObject& grid_in, const std::string& filename, const std::string& varname, const Date& date=Date());
//void add_attributes_for_variable(const int& ncid, const int& varid, const std::string& varname);
void create_latlon_dimensions(const int& ncid, const Grid2DObject& grid, int& did_lat, int& did_lon, int& vid_lat, int& vid_lon);
void create_time_dimension(const int& ncid, int& did_time, int& vid_time);
void readWind(const std::string& filename, const Date& date);
// Private variables
static const double plugin_nodata; //plugin specific nodata value, e.g. -999
......@@ -89,7 +92,9 @@ class NetCDFIO : public IOInterface {
const Config cfg;
std::map <MeteoGrids::Parameters, attributes> in_attributes, out_attributes;
std::string coordin, coordinparam, coordout, coordoutparam; //projection parameters
double in_dflt_TZ, out_dflt_TZ; //default time zones
double in_dflt_TZ, out_dflt_TZ; //default time zones
double in_time_offset, in_time_multiplier; //each schema defines its own time specification...
double out_time_offset, out_time_multiplier; //each schema defines its own time specification...
bool in_strict, out_strict;
std::vector<StationData> vecMetaData;
};
......
......@@ -303,6 +303,28 @@ size_t add_record(const int& ncid, const std::string& varname, const int& varid,
return dimlen;
}
bool get_recordMinMax(const int& ncid, const std::string& varname, const int& varid, double &min, double &max)
{
int dimid;
size_t dimlen;
get_dimension(ncid, varname, dimid, dimlen);
//check if record already exists
if (dimlen > 0) {
double *record_value = new double[dimlen];
read_data(ncid, varname, varid, record_value);
min = record_value[0];
max = record_value[dimlen-1];
delete[] record_value;
} else
return false; // data not found
return true;
}
// Finding a certain record variable value (e.g. timestamp) by retrieving all
// record values and then performing a linear search
size_t find_record(const int& ncid, const std::string& varname, const int& varid, const double& data)
......
......@@ -67,6 +67,7 @@ namespace ncpp {
const size_t& pos_start, const double * const data);
//Dealing with variables that have dimension NC_UNLIMITED
bool get_recordMinMax(const int& ncid, const std::string& varname, const int& varid, double &min, double &max);
size_t find_record(const int& ncid, const std::string& varname, const int& varid, const double& data);
size_t add_record(const int& ncid, const std::string& varname, const int& varid, const double& data);
void write_record(const int& ncid, const std::string& varname, const int& varid, const size_t& pos,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment