diff --git a/src/core.cpp b/src/core.cpp index 93635fea..170207ec 100644 --- a/src/core.cpp +++ b/src/core.cpp @@ -540,6 +540,24 @@ void core::config_forcing(pt::ptree &value) } nstations = _metdata->nstations(); + + // we need to estimate some + if(_metdata->missing_z()) + { + SPDLOG_WARN("No geopotential height field found in the netcdf file. Using the height of nearest triangle to estimate. This is almost certainly NOT what you want"); + + for(size_t i = 0; i < _mesh->size_faces(); i++) + { + for (auto face = _mesh->face(i); auto& s : face->stations()) + { + if(s->z() == -9999) + { + s->z(face->get_z()); + } + } + + } + } } else { std::vector ascii_data; diff --git a/src/metdata.cpp b/src/metdata.cpp index 6483963e..56725a4c 100644 --- a/src/metdata.cpp +++ b/src/metdata.cpp @@ -32,6 +32,7 @@ metdata::metdata(std::string mesh_proj4) is_first_timestep = true; _is_multipart_nc = false; _just_loaded_nc = false; + _missing_z = false; OGRSpatialReference srs; srs.importFromProj4(_mesh_proj4.c_str()); @@ -147,6 +148,8 @@ void metdata::load_from_netcdf(const std::string& path, const triangulation::bou _nc_ignored_variables.clear(); _nc->open_GEM(path); + + _missing_z = _nc->missing_z(); } catch(netCDF::exceptions::NcException& e) { OGRCoordinateTransformation::DestroyCT(coordTrans); @@ -213,8 +216,8 @@ void metdata::load_from_netcdf(const std::string& path, const triangulation::bou } else { - _start_time = _nc->get_start(); - _end_time = _nc->get_end(); + _file_start_time = _start_time = _nc->get_start(); + _file_end_time = _end_time = _nc->get_end(); _n_timesteps = _nc->get_ntimesteps(); } @@ -242,7 +245,9 @@ void metdata::load_from_netcdf(const std::string& path, const triangulation::bou SPDLOG_DEBUG("Initializing datastructure"); - auto e = _nc->get_z(); + netcdf::data e; + if(!missing_z()) + e = _nc->get_z(); // #pragma omp parallel for // hangs, unclear why, critical sections around the json and gdal calls @@ -271,8 +276,6 @@ void metdata::load_from_netcdf(const std::string& path, const triangulation::bou longitude = (*std::get(lon))[x]; } - z = (*e)[y][x]; - // Some Netcdf files have NaN grid squares, For these cases we will just insert a nullptr station and // don't add the station to the dD list which is the only way it ever gets to modules if ( std::isnan(latitude) || @@ -309,9 +312,17 @@ void metdata::load_from_netcdf(const std::string& path, const triangulation::bou } } + + if(missing_z()) + z = -9999; // estimate this later + else + z = (*e)[y][x]; + + std::shared_ptr s = std::make_shared(station_name, longitude, latitude, z, _variables); + //holds the corresponding x,y grid cell of the netcdf file s->_nc_x = x; s->_nc_y = y; @@ -836,4 +847,9 @@ void metdata::write_stations_to_shp(const std::string& fname) xy.emplace_back(itr->x(), itr->y()); } gis::xy2shp(xy, fname, _mesh_proj4); +} + +bool metdata::missing_z() +{ + return _missing_z; } \ No newline at end of file diff --git a/src/metdata.hpp b/src/metdata.hpp index ae83401f..4c9882f0 100644 --- a/src/metdata.hpp +++ b/src/metdata.hpp @@ -224,6 +224,12 @@ class metdata std::vector< std::shared_ptr>& stations(); + /** + * Netcdf was missing a geopotential height and we need to estimate it + * @return + */ + bool missing_z(); + private: struct ascii_data @@ -323,6 +329,8 @@ class metdata std::string _mesh_proj4; bool _is_geographic; // geographic mesh that requires further reprojection? + bool _missing_z; // missing a geopotential and need to estimate it from the triangles + // spatial searching data structure typedef CGAL::Simple_cartesian Kernel; typedef boost::tuple > Point_and_station; diff --git a/src/timeseries/netcdf.cpp b/src/timeseries/netcdf.cpp index 52f82f50..4ca4211c 100644 --- a/src/timeseries/netcdf.cpp +++ b/src/timeseries/netcdf.cpp @@ -29,6 +29,7 @@ netcdf::netcdf() { _is_open = false; _datetime_field=""; + _missing_z = false; } netcdf::~netcdf() { @@ -412,7 +413,18 @@ void netcdf::open_GEM(const std::string &file) // CF convention assumes that dim and coord have the same name, _lat_field = find_dim_by_standard_name("latitude"); _lon_field = find_dim_by_standard_name("longitude"); - _elevation_field = find_var_by_attr("geopotential_height"); + + try + { + _elevation_field = find_var_by_attr("geopotential_height"); + _missing_z = false; + }catch(...) + { + SPDLOG_WARN("No geopotential height field found. Using the height of nearest triangle. This is almost certainly NOT what you want"); + _elevation_field=""; + _missing_z = true; + } + xgrid = _data.getDim(_lon_field).getSize(); ygrid = _data.getDim(_lat_field).getSize(); @@ -438,6 +450,11 @@ void netcdf::open_GEM(const std::string &file) } +bool netcdf::missing_z() +{ + return _missing_z; +} + size_t netcdf::get_ntimesteps() { return _datetime_length; diff --git a/src/timeseries/netcdf.hpp b/src/timeseries/netcdf.hpp index 2ff2f4d6..2b0c8587 100644 --- a/src/timeseries/netcdf.hpp +++ b/src/timeseries/netcdf.hpp @@ -139,6 +139,9 @@ boost::posix_time::time_duration get_dt(); double get_var2D(std::string var, size_t x, size_t y); std::string get_var_standard_name(const std::string& variable) const; + // do we have a z data variable? + bool missing_z(); + /** * Returns the dimensionality of the spatial coordinates. 1D or 2D * @return @@ -166,6 +169,9 @@ int get_coord_dimensionality() const; boost::posix_time::time_duration _delta_t; // model timestep size. whole integer values, in seconds. e.g, 3600s bool _is_open; + // if we are missing z coords that need to be estimated later + bool _missing_z; + //if we are creating variables std::vector _dimVector; //we need this dimension var to create new variables