PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
VariableMetadata.cc
Go to the documentation of this file.
1 // Copyright (C) 2009--2023 Constantine Khroulev and Ed Bueler
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 
19 #include <set>
20 #include <algorithm>
21 #include <cmath>
22 
23 #include "pism/util/VariableMetadata.hh"
24 
25 #include "pism/util/Logger.hh"
26 #include "pism/util/error_handling.hh"
27 #include "pism/util/io/File.hh"
28 #include "pism/util/io/IO_Flags.hh"
29 
30 namespace pism {
31 
33  unsigned int ndims)
34  : m_n_spatial_dims(ndims),
35  m_unit_system(std::move(system)),
36  m_short_name(name),
37  m_time_independent(false),
38  m_output_type(io::PISM_NAT) {
39 
40  clear();
41 
42  // long_name is unset
43  // standard_name is unset
44  // coordinates is unset
45 
46  // valid_min and valid_max are unset
47 }
48 
49 /** Get the number of spatial dimensions.
50  */
52  return m_n_spatial_dims;
53 }
54 
55 /** A "time independent" variable will be saved to a NetCDF
56  variable which does not depend on the "time" dimension.
57  */
59  m_time_independent = flag;
60 
61  return *this;
62 }
63 
65  m_output_type = type;
66 
67  return *this;
68 }
69 
71  return m_time_independent;
72 }
73 
75  return m_output_type;
76 }
77 
79  return m_unit_system;
80 }
81 
82 //! @brief Check if the range `[min, max]` is a subset of `[valid_min, valid_max]`.
83 /*! Throws an exception if this check failed.
84  */
85 void VariableMetadata::check_range(const std::string &filename, double min, double max) const {
86 
87  auto units_string = get_string("units");
88  auto name_string = get_name();
89  const char *units = units_string.c_str(), *name = name_string.c_str(), *file = filename.c_str();
90  double eps = 1e-12;
91 
92  if (has_attribute("valid_min") and has_attribute("valid_max")) {
93  double valid_min = get_number("valid_min");
94  double valid_max = get_number("valid_max");
95  if ((min < valid_min - eps) or (max > valid_max + eps)) {
98  "some values of '%s' in '%s' are outside the valid range [%e, %e] (%s).\n"
99  "computed min = %e %s, computed max = %e %s",
100  name, file, valid_min, valid_max, units, min, units, max, units);
101  }
102  } else if (has_attribute("valid_min")) {
103  double valid_min = get_number("valid_min");
104  if (min < valid_min - eps) {
107  "some values of '%s' in '%s' are less than the valid minimum (%e %s).\n"
108  "computed min = %e %s, computed max = %e %s",
109  name, file, valid_min, units, min, units, max, units);
110  }
111  } else if (has_attribute("valid_max")) {
112  double valid_max = get_number("valid_max");
113  if (max > valid_max + eps) {
116  "some values of '%s' in '%s' are greater than the valid maximum (%e %s).\n"
117  "computed min = %e %s, computed max = %e %s",
118  name, file, valid_max, units, min, units, max, units);
119  }
120  }
121 }
122 
124  const std::vector<double> &zlevels)
125  : VariableMetadata(name, system),
126  m_x("x", system),
127  m_y("y", system),
128  m_z("z", system),
129  m_zlevels(zlevels) {
130 
131  m_x["axis"] = "X";
132  m_x["long_name"] = "X-coordinate in Cartesian system";
133  m_x["standard_name"] = "projection_x_coordinate";
134  m_x["units"] = "m";
135 
136  m_y["axis"] = "Y";
137  m_y["long_name"] = "Y-coordinate in Cartesian system";
138  m_y["standard_name"] = "projection_y_coordinate";
139  m_y["units"] = "m";
140 
141  if (m_zlevels.size() > 1) {
142  m_z.set_name("z"); // default; can be overridden easily
143  m_z["axis"] = "Z";
144  m_z["long_name"] = "Z-coordinate in Cartesian system";
145  m_z["units"] = "m";
146  m_z["positive"] = "up";
147 
148  m_n_spatial_dims = 3;
149  } else {
150  z().set_name("").clear();
151  m_n_spatial_dims = 2;
152  }
153 }
154 
155 const std::vector<double> &SpatialVariableMetadata::levels() const {
156  return m_zlevels;
157 }
158 
159 //! Report the range of a \b global Vec `v`.
160 void VariableMetadata::report_range(const Logger &log, double min, double max,
161  bool found_by_standard_name) {
162 
163  // units::Converter constructor will make sure that units are compatible.
164  units::Converter c(m_unit_system, this->get_string("units"), this->get_string("output_units"));
165  min = c(min);
166  max = c(max);
167 
168  std::string spacer(get_name().size(), ' ');
169 
170  if (has_attribute("standard_name")) {
171 
172  if (found_by_standard_name) {
173  log.message(2,
174  " FOUND "
175  " %s / standard_name=%-10s\n"
176  " %s \\ min,max = %9.3f,%9.3f (%s)\n",
177  get_name().c_str(), get_string("standard_name").c_str(), spacer.c_str(), min, max,
178  get_string("output_units").c_str());
179  } else {
180  log.message(2,
181  " FOUND "
182  " %s / WARNING! standard_name=%s is missing, found by short_name\n"
183  " %s \\ min,max = %9.3f,%9.3f (%s)\n",
184  get_name().c_str(), get_string("standard_name").c_str(), spacer.c_str(), min, max,
185  get_string("output_units").c_str());
186  }
187 
188  } else {
189  log.message(2,
190  " FOUND "
191  " %s / %-10s\n"
192  " %s \\ min,max = %9.3f,%9.3f (%s)\n",
193  get_name().c_str(), get_string("long_name").c_str(), spacer.c_str(), min, max,
194  get_string("output_units").c_str());
195  }
196 }
197 
199  return m_x;
200 }
201 
203  return m_y;
204 }
205 
207  return m_z;
208 }
209 
211  return m_x;
212 }
213 
215  return m_y;
216 }
217 
219  return m_z;
220 }
221 
222 //! Checks if an attribute is present. Ignores empty strings, except
223 //! for the "units" attribute.
224 bool VariableMetadata::has_attribute(const std::string &name) const {
225 
226  auto j = m_strings.find(name);
227  if (j != m_strings.end()) {
228  if (name != "units" and (j->second).empty()) {
229  return false;
230  }
231 
232  return true;
233  }
234 
235  return (m_doubles.find(name) != m_doubles.end());
236 }
237 
239  return not(this->all_strings().empty() and this->all_doubles().empty());
240 }
241 
242 VariableMetadata &VariableMetadata::set_name(const std::string &name) {
243  m_short_name = name;
244 
245  return *this;
246 }
247 
248 //! Set a scalar attribute to a single (scalar) value.
249 VariableMetadata &VariableMetadata::set_number(const std::string &name, double value) {
250  return set_numbers(name, {value});
251 }
252 
253 //! Set a scalar attribute to a single (scalar) value.
254 VariableMetadata &VariableMetadata::set_numbers(const std::string &name, const std::vector<double> &values) {
255  m_doubles[name] = values;
256 
257  return *this;
258 }
259 
261  m_doubles.clear();
262  m_strings.clear();
263 
264  return *this;
265 }
266 
267 std::string VariableMetadata::get_name() const {
268  return m_short_name;
269 }
270 
271 //! Get a single-valued scalar attribute.
272 double VariableMetadata::get_number(const std::string &name) const {
273  auto j = m_doubles.find(name);
274  if (j != m_doubles.end()) {
275  return (j->second)[0];
276  }
277 
279  "variable \"%s\" does not have a double attribute \"%s\"",
280  get_name().c_str(), name.c_str());
281 }
282 
283 //! Get an array-of-doubles attribute.
284 std::vector<double> VariableMetadata::get_numbers(const std::string &name) const {
285  auto j = m_doubles.find(name);
286  if (j != m_doubles.end()) {
287  return j->second;
288  }
289 
290  return {};
291 }
292 
293 const std::map<std::string, std::string> &VariableMetadata::all_strings() const {
294  return m_strings;
295 }
296 
297 const std::map<std::string, std::vector<double> > &VariableMetadata::all_doubles() const {
298  return m_doubles;
299 }
300 
301 /*!
302  * Set units that may not be supported by UDUNITS.
303  *
304  * For example: "Pa s^(1/3)" (ice hardness units with Glen exponent n=3).
305  */
307  m_strings["units"] = value;
308  m_strings["output_units"] = value;
309 
310  return *this;
311 }
312 
313 //! Set a string attribute.
314 VariableMetadata &VariableMetadata::set_string(const std::string &name, const std::string &value) {
315 
316  if (name == "units") {
317  // create a dummy object to validate the units string
318  units::Unit tmp(m_unit_system, value);
319 
320  m_strings[name] = value;
321  m_strings["output_units"] = value;
322  } else if (name == "output_units") {
323  m_strings[name] = value;
324 
325  units::Unit internal(m_unit_system, get_string("units"));
326  units::Unit output(m_unit_system, value);
327 
328  if (not units::are_convertible(internal, output)) {
329  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "units \"%s\" and \"%s\" are not compatible",
330  get_string("units").c_str(), value.c_str());
331  }
332  } else if (name == "short_name") {
333  set_name(name);
334  } else {
335  m_strings[name] = value;
336  }
337 
338  return *this;
339 }
340 
341 //! Get a string attribute.
342 /*!
343  * Returns an empty string if an attribute is not set.
344  */
345 std::string VariableMetadata::get_string(const std::string &name) const {
346  if (name == "short_name") {
347  return get_name();
348  }
349 
350  auto j = m_strings.find(name);
351  if (j != m_strings.end()) {
352  return j->second;
353  }
354 
355  return {};
356 }
357 
358 void VariableMetadata::report_to_stdout(const Logger &log, int verbosity_threshold) const {
359 
360  const auto &strings = this->all_strings();
361  const auto &doubles = this->all_doubles();
362 
363  // Find the maximum name length so that we can pad output below:
364  size_t max_name_length = 0;
365  for (const auto &s : strings) {
366  max_name_length = std::max(max_name_length, s.first.size());
367  }
368  for (const auto &d : doubles) {
369  max_name_length = std::max(max_name_length, d.first.size());
370  }
371 
372  // Print text attributes:
373  for (const auto &s : strings) {
374  std::string name = s.first;
375  std::string value = s.second;
376  std::string padding(max_name_length - name.size(), ' ');
377 
378  if (value.empty()) {
379  continue;
380  }
381 
382  log.message(verbosity_threshold, " %s%s = \"%s\"\n",
383  name.c_str(), padding.c_str(), value.c_str());
384  }
385 
386  // Print double attributes:
387  for (const auto &d : doubles) {
388  std::string name = d.first;
389  std::vector<double> values = d.second;
390  std::string padding(max_name_length - name.size(), ' ');
391 
392  if (values.empty()) {
393  continue;
394  }
395 
396  const double large = 1.0e7;
397  const double small = 1.0e-4;
398  if ((std::fabs(values[0]) >= large) || (std::fabs(values[0]) <= small)) {
399  // use scientific notation if a number is big or small
400  log.message(verbosity_threshold, " %s%s = %12.3e\n",
401  name.c_str(), padding.c_str(), values[0]);
402  } else {
403  log.message(verbosity_threshold, " %s%s = %12.5f\n",
404  name.c_str(), padding.c_str(), values[0]);
405  }
406 
407  }
408 }
409 
410 ConstAttribute::ConstAttribute(const VariableMetadata *var, const std::string &name)
411  : m_name(name), m_var(const_cast<VariableMetadata*>(var)) {
412 }
413 
415  : m_name(std::move(a.m_name)), m_var(a.m_var) {
416  a.m_name.clear();
417  a.m_var = nullptr;
418 }
419 
420 ConstAttribute::operator std::string() const {
421  return m_var->get_string(m_name);
422 }
423 
424 ConstAttribute::operator double() const {
425  auto values = m_var->get_numbers(m_name);
426  if (values.size() == 1) {
427  return values[0];
428  }
430  "%s:%s has more than one value",
431  m_var->get_name().c_str(), m_name.c_str());
432 }
433 
434 ConstAttribute::operator std::vector<double> () const {
435  return m_var->get_numbers(m_name);
436 }
437 
438 void Attribute::operator=(const std::string &value) {
439  m_var->set_string(m_name, value);
440 }
441 void Attribute::operator=(const std::initializer_list<double> &value) {
442  m_var->set_numbers(m_name, value);
443 }
444 
445 void Attribute::operator=(const std::vector<double> &value) {
446  m_var->set_numbers(m_name, value);
447 }
448 
449 } // end of namespace pism
void operator=(const std::string &value)
ConstAttribute(const ConstAttribute &)=delete
VariableMetadata * m_var
void message(int threshold, const char format[],...) const __attribute__((format(printf
Print a message to the log.
Definition: Logger.cc:49
A basic logging class.
Definition: Logger.hh:40
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
std::vector< double > m_zlevels
const std::vector< double > & levels() const
SpatialVariableMetadata(units::System::Ptr system, const std::string &name, const std::vector< double > &zlevels={0.0})
VariableMetadata & clear()
Clear all attributes.
void report_to_stdout(const Logger &log, int verbosity_threshold) const
VariableMetadata & set_name(const std::string &name)
unsigned int n_spatial_dimensions() const
VariableMetadata & set_number(const std::string &name, double value)
Set a scalar attribute to a single (scalar) value.
VariableMetadata(const std::string &name, units::System::Ptr system, unsigned int ndims=0)
units::System::Ptr m_unit_system
The unit system to use.
double get_number(const std::string &name) const
Get a single-valued scalar attribute.
std::string get_string(const std::string &name) const
Get a string attribute.
const std::map< std::string, std::string > & all_strings() const
std::map< std::string, std::string > m_strings
string and boolean attributes
void report_range(const Logger &log, double min, double max, bool found_by_standard_name)
Report the range of a global Vec v.
units::System::Ptr unit_system() const
io::Type get_output_type() const
std::map< std::string, std::vector< double > > m_doubles
scalar and array attributes
bool get_time_independent() const
VariableMetadata & set_units_without_validation(const std::string &value)
std::vector< double > get_numbers(const std::string &name) const
Get an array-of-doubles attribute.
VariableMetadata & set_output_type(io::Type type)
bool has_attribute(const std::string &name) const
VariableMetadata & units(const std::string &input)
VariableMetadata & set_numbers(const std::string &name, const std::vector< double > &values)
Set a scalar attribute to a single (scalar) value.
std::string get_name() const
VariableMetadata & set_string(const std::string &name, const std::string &value)
Set a string attribute.
void check_range(const std::string &filename, double min, double max) const
Check if the range [min, max] is a subset of [valid_min, valid_max].
const std::map< std::string, std::vector< double > > & all_doubles() const
VariableMetadata & set_time_independent(bool flag)
std::shared_ptr< System > Ptr
Definition: Units.hh:47
#define PISM_ERROR_LOCATION
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
double min(const array::Scalar &input)
Finds minimum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:193
@ PISM_NAT
Definition: IO_Flags.hh:46
bool are_convertible(const Unit &u1, const Unit &u2)
Definition: Units.cc:242