Skip to content

Commit

Permalink
If a netcdf is missing geopotential height, estimate the cell centre …
Browse files Browse the repository at this point in the history
…elevation as that of the closest triangle's elevation. Should almost never be used
  • Loading branch information
Chrismarsh committed Dec 17, 2024
1 parent dc7db80 commit 397d3fd
Show file tree
Hide file tree
Showing 5 changed files with 71 additions and 6 deletions.
18 changes: 18 additions & 0 deletions src/core.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<metdata::ascii_metdata> ascii_data;
Expand Down
26 changes: 21 additions & 5 deletions src/metdata.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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();
}

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -271,8 +276,6 @@ void metdata::load_from_netcdf(const std::string& path, const triangulation::bou
longitude = (*std::get<netcdf::vec>(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) ||
Expand Down Expand Up @@ -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<station> s = std::make_shared<station>(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;

Expand Down Expand Up @@ -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;
}
8 changes: 8 additions & 0 deletions src/metdata.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -224,6 +224,12 @@ class metdata

std::vector< std::shared_ptr<station>>& stations();

/**
* Netcdf was missing a geopotential height and we need to estimate it
* @return
*/
bool missing_z();

private:

struct ascii_data
Expand Down Expand Up @@ -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<double> Kernel;
typedef boost::tuple<Kernel::Point_2, std::shared_ptr<station> > Point_and_station;
Expand Down
19 changes: 18 additions & 1 deletion src/timeseries/netcdf.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ netcdf::netcdf()
{
_is_open = false;
_datetime_field="";
_missing_z = false;
}
netcdf::~netcdf()
{
Expand Down Expand Up @@ -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();
Expand All @@ -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;
Expand Down
6 changes: 6 additions & 0 deletions src/timeseries/netcdf.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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<netCDF::NcDim> _dimVector; //we need this dimension var to create new variables
Expand Down

0 comments on commit 397d3fd

Please sign in to comment.