WSL/SLF GitLab Repository

Commit 7752a3c9 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The NetCDFIO plugin can now properly handle multiple files in one directory...

The NetCDFIO plugin can now properly handle multiple files in one directory specified by "GRID2DPATH" or one single file specified by "GRID2DFILE". Small code cleanup in CNRMIO and CosmoXMLIO.
parent 6d0cf546
......@@ -96,7 +96,6 @@ class CNRMIO : public IOInterface {
static const std::string cf_time, cf_units, cf_days, cf_hours, cf_seconds, cf_latitude, cf_longitude, cf_altitude, cf_ta, cf_rh, cf_p;
static const std::string cnrm_points, cnrm_latitude, cnrm_longitude, cnrm_altitude, cnrm_aspect, cnrm_slope, cnrm_ta, cnrm_rh, cnrm_vw, cnrm_dw, cnrm_qair;
static const std::string cnrm_co2air, cnrm_theorsw, cnrm_neb, cnrm_psum, cnrm_snowf, cnrm_swr_direct, cnrm_swr_diffuse, cnrm_p, cnrm_ilwr, cnrm_timestep;
static const std::string ecmwf_ta, ecmwf_p, ecmwf_iswr, ecmwf_ilwr, ecmwf_psum, ecmwf_u10, ecmwf_v10, ecmwf_td;
static std::map<std::string, size_t> paramname; ///<Associate a name with meteo parameters in Parameters
static std::map<std::string, std::string> map_name; ///Associate MeteoIO parameter names with CNRM parameter names
......
......@@ -211,8 +211,7 @@ void CosmoXMLIO::scanMeteoPath(const std::string& meteopath_in, std::vector< st
{
meteo_files.clear();
std::list<std::string> dirlist;
IOUtils::readDirectory(meteopath_in, dirlist, meteo_ext);
std::list<std::string> dirlist = IOUtils::readDirectory(meteopath_in, meteo_ext);
dirlist.sort();
//Check date in every filename and cache it
......
......@@ -82,7 +82,8 @@ const std::string NetCDFIO::cf_latitude = "latitude";
const std::string NetCDFIO::cf_longitude = "longitude";
const std::string NetCDFIO::cf_altitude = "z";
NetCDFIO::NetCDFIO(const std::string& configfile) : cfg(configfile), in_attributes(), out_attributes(), coordin(), coordinparam(), coordout(), coordoutparam(),
NetCDFIO::NetCDFIO(const std::string& configfile) : cfg(configfile), cache_meteo_files(), in_attributes(), out_attributes(),
coordin(), coordinparam(), coordout(), coordoutparam(),
in_dflt_TZ(0.), out_dflt_TZ(0.), in_time_offset(IOUtils::nodata), in_time_multiplier(IOUtils::nodata),
dem_altimeter(false), in_strict(false), out_strict(false), vecMetaData()
{
......@@ -90,7 +91,8 @@ NetCDFIO::NetCDFIO(const std::string& configfile) : cfg(configfile), in_attribut
parseInputOutputSection();
}
NetCDFIO::NetCDFIO(const Config& cfgreader) : cfg(cfgreader), in_attributes(), out_attributes(), coordin(), coordinparam(), coordout(), coordoutparam(),
NetCDFIO::NetCDFIO(const Config& cfgreader) : cfg(cfgreader), cache_meteo_files(), in_attributes(), out_attributes(),
coordin(), coordinparam(), coordout(), coordoutparam(),
in_dflt_TZ(0.), out_dflt_TZ(0.), in_time_offset(IOUtils::nodata), in_time_multiplier(IOUtils::nodata),
dem_altimeter(false), in_strict(false), out_strict(false), vecMetaData()
{
......@@ -114,6 +116,10 @@ void NetCDFIO::parseInputOutputSection()
string out_schema = "ECMWF";
out_schema = IOUtils::strToUpper( cfg.get("NETCDF_SCHEMA", "Output", IOUtils::nothrow) );
initAttributesMap(out_schema, out_attributes);
string in_grid2d_path;
cfg.getValue("GRID2DPATH", "Input", in_grid2d_path, IOUtils::nothrow);
if (!in_grid2d_path.empty()) scanMeteoPath(in_grid2d_path, cache_meteo_files);
}
void NetCDFIO::initAttributesMap(const std::string& schema, std::map<MeteoGrids::Parameters, attributes> &attr)
......@@ -174,10 +180,30 @@ void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const std::string& arguments)
}
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
{
if (!cache_meteo_files.empty()) {
for (size_t ii=0; ii<cache_meteo_files.size(); ii++) {
const Date date_start = cache_meteo_files[ii].first.first;
const Date date_end = cache_meteo_files[ii].first.second;
if (date>=date_start && date<=date_end) {
const string filename = cache_meteo_files[ii].second;
read2DGrid(grid_out, parameter, date, filename);
return;
}
}
//the date was not found
string in_grid2d_path;
cfg.getValue("GRID2DPATH", "Input", in_grid2d_path);
throw InvalidArgumentException("No Gridded data found for "+date.toString(Date::ISO)+"in '"+in_grid2d_path+"'", AT);
} else {
const string filename = cfg.get("GRID2DFILE", "Input");
read2DGrid(grid_out, parameter, date, filename);
}
}
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date, const std::string& filename)
{
grid_out.clear();
const string filename = cfg.get("GRID2DFILE", "Input"); //HACK: also allow using GRID2DPATH
if (read2DGrid_internal(grid_out, filename, parameter, date)) return; //schema naming
if (parameter==MeteoGrids::VW || parameter==MeteoGrids::DW) { //VW, DW
......@@ -381,6 +407,56 @@ void NetCDFIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parame
}
}
//custom function for sorting cache_meteo_files
struct sort_pred {
bool operator()(const std::pair<std::pair<Date,Date>,string> &left, const std::pair<std::pair<Date,Date>,string> &right) {
if (left.first.first < right.first.first) return true;
if (left.first.first > right.first.first) return false;
return left.first.second < right.first.second; //date_start equallity case
}
};
void NetCDFIO::scanMeteoPath(const std::string& meteopath_in, std::vector< std::pair<std::pair<mio::Date, mio::Date>, std::string> > &meteo_files)
{
meteo_files.clear();
string meteo_ext(".nc");
cfg.getValue("METEO_EXT", "INPUT", meteo_ext, IOUtils::nothrow);
std::list<std::string> dirlist = IOUtils::readDirectory(meteopath_in, meteo_ext);
if (dirlist.empty()) return; //nothing to do if the directory is empty, we will transparently swap to using GRID2DFILE
dirlist.sort();
//Check date range in every filename and cache it
std::list<std::string>::const_iterator it = dirlist.begin();
while ((it != dirlist.end())) {
const std::string& filename = meteopath_in + "/" + *it;
if (!IOUtils::fileExists(filename)) throw AccessException(filename, AT); //prevent invalid filenames
int ncid;
ncpp::open_file(filename, NC_NOWRITE, ncid);
getTimeTransform(ncid, in_time_offset, in_time_multiplier); //always re-read offset and multiplier, it might be different for each file
double min, max;
const bool status = ncpp::get_dimensionMinMax(ncid, NetCDFIO::cf_time, min, max);
ncpp::close_file(filename, ncid); //no need to keep file open anymore
if (!status) throw IOException("Could not get min/max time for file '"+filename+"'", AT);
const Date d_min(min*in_time_multiplier + in_time_offset, in_dflt_TZ);
const Date d_max(max*in_time_multiplier + in_time_offset, in_dflt_TZ);
const std::pair<Date, Date> dateRange(d_min, d_max);
const std::pair<std::pair<Date, Date>,std::string> tmp(dateRange, filename);
meteo_files.push_back(tmp);
it++;
}
std::sort(meteo_files.begin(), meteo_files.end(), sort_pred());
//now handle overlaping files: truncate the end date of the file starting earlier
for (size_t ii=0; ii<(meteo_files.size()-1); ii++) {
if (meteo_files[ii].first.second > meteo_files[ii+1].first.first)
meteo_files[ii].first.second = meteo_files[ii+1].first.first;
}
}
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);
......@@ -441,15 +517,14 @@ bool NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& fi
throw InvalidArgumentException("File \'"+filename+"\' does not contain a time dimension", AT);
} else {
if (date_requested) {
const int timeID = dimid[time_index];
if (in_time_offset==IOUtils::nodata || in_time_multiplier==IOUtils::nodata) {
getTimeTransform(ncid, timeID, in_time_offset, in_time_multiplier);
getTimeTransform(ncid, in_time_offset, in_time_multiplier);
}
const double timestamp = (date.getJulian() - in_time_offset) / in_time_multiplier;
const size_t pos = ncpp::find_record(ncid, NetCDFIO::cf_time, timeID, timestamp);
const size_t pos = ncpp::find_record(ncid, NetCDFIO::cf_time, timestamp);
if (pos == IOUtils::npos) {
double min, max;
const bool status = ncpp::get_recordMinMax(ncid, NetCDFIO::cf_time, timeID, min, max);
const bool status = ncpp::get_dimensionMinMax(ncid, NetCDFIO::cf_time, min, max);
if (status) {
Date d_min, d_max;
d_min.setDate(min*in_time_multiplier + in_time_offset, in_dflt_TZ);
......@@ -620,11 +695,11 @@ void NetCDFIO::write2DGrid_internal(Grid2DObject grid_in, const std::string& fil
delete[] lat_array; delete[] lon_array; delete[] data;
}
void NetCDFIO::getTimeTransform(const int& ncid, const int& varid, double &time_offset, double &time_multiplier) const
void NetCDFIO::getTimeTransform(const int& ncid, double &time_offset, double &time_multiplier) const
{
string time_units;
ncpp::get_attribute(ncid, NetCDFIO::cf_time, varid, "units", time_units);
ncpp::get_DimAttribute(ncid, NetCDFIO::cf_time, "units", time_units);
std::vector<std::string> vecString;
const size_t nrWords = IOUtils::readLineToVec(time_units, vecString);
if (nrWords<3 || nrWords>4) throw InvalidArgumentException("Invalid format for time units: \'"+time_units+"\'", AT);
......
......@@ -72,16 +72,18 @@ class NetCDFIO : public IOInterface {
} attributes;
void initAttributesMap(const std::string& schema, std::map<MeteoGrids::Parameters, attributes> &attr);
void scanMeteoPath(const std::string& meteopath_in, std::vector< std::pair<std::pair<mio::Date, mio::Date>, std::string> > &meteo_files);
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);
//void read2DGrid_meta(const std::string& filename, const std::string& varname, int &ncid, int &varid, std::vector<size_t> &dimlen, size_t &time_index, size_t &lat_index, size_t &lon_index, double &missing_value,double &lat, double &lon);
void read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date, const std::string& filename);
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(), const bool& isPrecip=false);
void write2DGrid_internal(Grid2DObject grid_in, const std::string& filename, const attributes& attr, const Date& date=Date(), const bool& isPrecip=false);
void add_attributes_for_variable(const int& ncid, const int& varid, const attributes& attr, const double& nodata_out);
void getTimeTransform(const int& ncid, const int& varid, double &time_offset, double &time_multiplier) const;
void getTimeTransform(const int& ncid, double &time_offset, double &time_multiplier) const;
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);
......@@ -92,6 +94,7 @@ class NetCDFIO : public IOInterface {
static const std::string cf_time, cf_latitude, cf_longitude, cf_altitude;
const Config cfg;
std::vector< std::pair<std::pair<Date,Date>, std::string> > cache_meteo_files; //cache of meteo files in METEOPATH
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
......
......@@ -66,16 +66,57 @@ void get_dimension(const int& ncid, const std::string& dimname, int& dimid, size
throw IOException("Could not retrieve length for dimension '" + dimname + "': " + nc_strerror(status), AT);
}
void get_attribute(const int& ncid, const std::string& varname, const int& varid, const std::string& attr_name, std::string& attr_value)
void get_DimAttribute(const int& ncid, const std::string& dimname, const std::string& attr_name, std::string& attr_value)
{
int dimid;
get_dimension(ncid, dimname, dimid);
size_t attr_len;
int status = nc_inq_attlen (ncid, dimid, attr_name.c_str(), &attr_len);
if (status != NC_NOERR)
throw IOException("Could not retrieve attribute '" + attr_name + "' for var '" + dimname + "': " + nc_strerror(status), AT);
char* value = new char[attr_len + 1]; // +1 for trailing null
status = nc_get_att_text(ncid, dimid, attr_name.c_str(), value);
if (status != NC_NOERR)
throw IOException("Could not read attribute '" + attr_name + "' for var '" + dimname + "': " + nc_strerror(status), AT);
value[attr_len] = '\0';
attr_value = string(value);
delete[] value;
}
void get_VarAttribute(const int& ncid, const std::string& varname, const std::string& attr_name, std::string& attr_value)
{
int varid;
get_variable(ncid, varname, varid);
size_t attr_len;
int status = nc_inq_attlen (ncid, varid, attr_name.c_str(), &attr_len);
if (status != NC_NOERR)
throw IOException("Could not retrieve attribute '" + attr_name + "' for var '" + varname + "': " + nc_strerror(status), AT);
char* value = new char[attr_len + 1]; // +1 for trailing null
status = nc_get_att_text(ncid, varid, attr_name.c_str(), value);
if (status != NC_NOERR)
throw IOException("Could not read attribute '" + attr_name + "' for var '" + varname + "': " + nc_strerror(status), AT);
value[attr_len] = '\0';
attr_value = string(value);
delete[] value;
}
void get_attribute(const int& ncid, const std::string& varname, const int& varid, const std::string& attr_name, std::string& attr_value)
{
size_t attr_len;
int status = nc_inq_attlen (ncid, varid, attr_name.c_str(), &attr_len);
if (status != NC_NOERR)
throw IOException("Could not retrieve attribute '" + attr_name + "' for var '" + varname + "': " + nc_strerror(status), AT);
char* value = new char[attr_len + 1]; // +1 for trailing null
status = nc_get_att_text(ncid, varid, attr_name.c_str(), value);
if (status != NC_NOERR)
throw IOException("Could not read attribute '" + attr_name + "' for var '" + varname + "': " + nc_strerror(status), AT);
......@@ -323,13 +364,34 @@ 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)
bool get_dimensionMinMax(const int& ncid, const std::string& varname, 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, dimid, record_value);
min = record_value[0];
max = record_value[dimlen-1];
delete[] record_value;
} else
return false; // data not found
return true;
}
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];
......@@ -341,12 +403,36 @@ bool get_recordMinMax(const int& ncid, const std::string& varname, const int& va
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 double& data)
{
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, dimid, record_value);
for (size_t ii=0; ii<dimlen; ii++) {
if (record_value[ii] == data) {
delete[] record_value;
return ii;
}
}
delete[] record_value;
}
return IOUtils::npos; // data not found
}
size_t find_record(const int& ncid, const std::string& varname, const int& varid, const double& data)
{
int dimid;
......
......@@ -45,6 +45,8 @@ namespace ncpp {
//Adding attributes
void add_attribute(const int& ncid, const int& varid, const std::string& attr_name, const std::string& attr_value);
void add_attribute(const int& ncid, const int& varid, const std::string& attr_name, const double& attr_value);
void get_DimAttribute(const int& ncid, const std::string& dimname, const std::string& attr_name, std::string& attr_value);
void get_VarAttribute(const int& ncid, const std::string& varname, const std::string& attr_name, std::string& attr_value);
void get_attribute(const int& ncid, const std::string& varname, const int& varid, const std::string& attr_name, std::string& attr_value);
void get_attribute(const int& ncid, const std::string& varname, const int& varid, const std::string& attr_name, double& attr_value);
bool check_attribute(const int& ncid, const int& varid, const std::string& attr_name);
......@@ -70,7 +72,9 @@ namespace ncpp {
const size_t& pos_start, const int * const data);
//Dealing with variables that have dimension NC_UNLIMITED
bool get_dimensionMinMax(const int& ncid, const std::string& varname, double &min, double &max);
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 double& data);
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,
......
Supports Markdown
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