PISM, A Parallel Ice Sheet Model  stable v2.0.5 committed by Constantine Khrulev on 2022-10-14 09:56:26 -0800
VariableMetadata.cc
Go to the documentation of this file.
1 // Copyright (C) 2009--2022 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 "VariableMetadata.hh"
24 #include "pism/util/io/File.hh"
25 #include "pism_options.hh"
26 #include "IceGrid.hh"
27 
28 #include "ConfigInterface.hh"
29 #include "error_handling.hh"
30 #include "pism/util/Logger.hh"
31 #include "pism_utilities.hh"
32 
33 namespace pism {
34 
36  unsigned int ndims)
37  : m_n_spatial_dims(ndims),
38  m_unit_system(std::move(system)),
39  m_short_name(name),
40  m_time_independent(false),
41  m_output_type(PISM_NAT) {
42 
45 
46  // long_name is unset
47  // standard_name is unset
48  // pism_intent is unset
49  // coordinates is unset
50 
51  // valid_min and valid_max are unset
52 }
53 
54 /** Get the number of spatial dimensions.
55  */
57  return m_n_spatial_dims;
58 }
59 
61  m_time_independent = flag;
62 }
63 
65  m_output_type = type;
66 }
67 
69  return m_time_independent;
70 }
71 
73  return m_output_type;
74 }
75 
77  return m_unit_system;
78 }
79 
80 //! @brief Check if the range `[min, max]` is a subset of `[valid_min, valid_max]`.
81 /*! Throws an exception if this check failed.
82  */
83 void VariableMetadata::check_range(const std::string &filename, double min, double max) const {
84 
85  auto units_string = get_string("units");
86  auto name_string = get_name();
87  const char
88  *units = units_string.c_str(),
89  *name = name_string.c_str(),
90  *file = filename.c_str();
91  double eps = 1e-12;
92 
93  if (has_attribute("valid_min") and has_attribute("valid_max")) {
94  double valid_min = get_number("valid_min");
95  double valid_max = get_number("valid_max");
96  if ((min < valid_min - eps) or (max > valid_max + eps)) {
97  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "some values of '%s' in '%s' are outside the valid range [%e, %e] (%s).\n"
98  "computed min = %e %s, computed max = %e %s",
99  name, file,
100  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) {
105  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "some values of '%s' in '%s' are less than the valid minimum (%e %s).\n"
106  "computed min = %e %s, computed max = %e %s",
107  name, file,
108  valid_min, units, min, units, max, units);
109  }
110  } else if (has_attribute("valid_max")) {
111  double valid_max = get_number("valid_max");
112  if (max > valid_max + eps) {
113  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "some values of '%s' in '%s' are greater than the valid maximum (%e %s).\n"
114  "computed min = %e %s, computed max = %e %s",
115  name, file,
116  valid_max, units, min, units, max, units);
117  }
118  }
119 }
120 
122  const std::string &name,
123  const std::vector<double> &zlevels)
124  : VariableMetadata(name, system),
125  m_x("x", system),
126  m_y("y", system),
127  m_z("z", system),
128  m_zlevels(zlevels) {
129 
130  m_x["axis"] = "X";
131  m_x["long_name"] = "X-coordinate in Cartesian system";
132  m_x["standard_name"] = "projection_x_coordinate";
133  m_x["units"] = "m";
134 
135  m_y["axis"] = "Y";
136  m_y["long_name"] = "Y-coordinate in Cartesian system";
137  m_y["standard_name"] = "projection_y_coordinate";
138  m_y["units"] = "m";
139 
140  m_z["axis"] = "Z";
141  m_z["long_name"] = "Z-coordinate in Cartesian system";
142  m_z["units"] = "m";
143  m_z["positive"] = "up";
144 
145  if (m_zlevels.size() > 1) {
146  z().set_name("z"); // default; can be overridden easily
147  m_n_spatial_dims = 3;
148  } else {
149  z().set_name("");
150  m_n_spatial_dims = 2;
151  }
152 }
153 
154 const std::vector<double>& SpatialVariableMetadata::levels() const {
155  return m_zlevels;
156 }
157 
158 //! Report the range of a \b global Vec `v`.
159 void VariableMetadata::report_range(const Logger &log, double min, double max,
160  bool found_by_standard_name) {
161 
162  // units::Converter constructor will make sure that units are compatible.
164  this->get_string("units"),
165  this->get_string("glaciological_units"));
166  min = c(min);
167  max = c(max);
168 
169  std::string spacer(get_name().size(), ' ');
170 
171  if (has_attribute("standard_name")) {
172 
173  if (found_by_standard_name) {
174  log.message(2,
175  " %s / standard_name=%-10s\n"
176  " %s \\ min,max = %9.3f,%9.3f (%s)\n",
177  get_name().c_str(),
178  get_string("standard_name").c_str(), spacer.c_str(), min, max,
179  get_string("glaciological_units").c_str());
180  } else {
181  log.message(2,
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(),
185  get_string("standard_name").c_str(), spacer.c_str(), min, max,
186  get_string("glaciological_units").c_str());
187  }
188 
189  } else {
190  log.message(2,
191  " %s / %-10s\n"
192  " %s \\ min,max = %9.3f,%9.3f (%s)\n",
193  get_name().c_str(),
194  get_string("long_name").c_str(), spacer.c_str(), min, max,
195  get_string("glaciological_units").c_str());
196  }
197 }
198 
200  return m_x;
201 }
202 
204  return m_y;
205 }
206 
208  return m_z;
209 }
210 
212  return m_x;
213 }
214 
216  return m_y;
217 }
218 
220  return m_z;
221 }
222 
223 //! Checks if an attribute is present. Ignores empty strings, except
224 //! for the "units" attribute.
225 bool VariableMetadata::has_attribute(const std::string &name) const {
226 
227  auto j = m_strings.find(name);
228  if (j != m_strings.end()) {
229  if (name != "units" and (j->second).empty()) {
230  return false;
231  }
232 
233  return true;
234  }
235 
236  return (m_doubles.find(name) != m_doubles.end());
237 }
238 
240  return not (this->all_strings().empty() and this->all_doubles().empty());
241 }
242 
243 void VariableMetadata::set_name(const std::string &name) {
244  m_short_name = name;
245 }
246 
247 //! Set a scalar attribute to a single (scalar) value.
248 void VariableMetadata::set_number(const std::string &name, double value) {
249  m_doubles[name] = std::vector<double>(1, value);
250 }
251 
252 //! Set a scalar attribute to a single (scalar) value.
253 void VariableMetadata::set_numbers(const std::string &name, const std::vector<double> &values) {
254  m_doubles[name] = values;
255 }
256 
258  m_doubles.clear();
259 }
260 
262  m_strings.clear();
263 }
264 
265 std::string VariableMetadata::get_name() const {
266  return m_short_name;
267 }
268 
269 //! Get a single-valued scalar attribute.
270 double VariableMetadata::get_number(const std::string &name) const {
271  auto j = m_doubles.find(name);
272  if (j != m_doubles.end()) {
273  return (j->second)[0];
274  }
275 
277  "variable \"%s\" does not have a double attribute \"%s\"",
278  get_name().c_str(), name.c_str());
279 }
280 
281 //! Get an array-of-doubles attribute.
282 std::vector<double> VariableMetadata::get_numbers(const std::string &name) const {
283  auto j = m_doubles.find(name);
284  if (j != m_doubles.end()) {
285  return j->second;
286  }
287 
288  return {};
289 }
290 
292  return m_strings;
293 }
294 
296  return m_doubles;
297 }
298 
299 /*!
300  * Set units that may not be supported by UDUNITS.
301  *
302  * For example: "Pa s^(1/3)" (ice hardness units with Glen exponent n=3).
303  */
304 void VariableMetadata::set_units_without_validation(const std::string &value) {
305  m_strings["units"] = value;
306  m_strings["glaciological_units"] = value;
307 }
308 
309 //! Set a string attribute.
310 void VariableMetadata::set_string(const std::string &name, const std::string &value) {
311 
312  if (name == "units") {
313  // create a dummy object to validate the units string
314  units::Unit tmp(m_unit_system, value);
315 
316  m_strings[name] = value;
317  m_strings["glaciological_units"] = value;
318  } else if (name == "glaciological_units") {
319  m_strings[name] = value;
320 
321  units::Unit internal(m_unit_system, get_string("units"));
322  units::Unit glaciological(m_unit_system, value);
323 
324  if (not units::are_convertible(internal, glaciological)) {
325  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "units \"%s\" and \"%s\" are not compatible",
326  get_string("units").c_str(), value.c_str());
327  }
328  } else if (name == "short_name") {
329  set_name(name);
330  } else {
331  m_strings[name] = value;
332  }
333 }
334 
335 //! Get a string attribute.
336 /*!
337  * Returns an empty string if an attribute is not set.
338  */
339 std::string VariableMetadata::get_string(const std::string &name) const {
340  if (name == "short_name") {
341  return get_name();
342  }
343 
344  auto j = m_strings.find(name);
345  if (j != m_strings.end()) {
346  return j->second;
347  }
348 
349  return {};
350 }
351 
352 void VariableMetadata::report_to_stdout(const Logger &log, int verbosity_threshold) const {
353 
354  const VariableMetadata::StringAttrs &strings = this->all_strings();
355  const VariableMetadata::DoubleAttrs &doubles = this->all_doubles();
356 
357  // Find the maximum name length so that we can pad output below:
358  size_t max_name_length = 0;
359  for (const auto &s : strings) {
360  max_name_length = std::max(max_name_length, s.first.size());
361  }
362  for (const auto &d : doubles) {
363  max_name_length = std::max(max_name_length, d.first.size());
364  }
365 
366  // Print text attributes:
367  for (const auto &s : strings) {
368  std::string name = s.first;
369  std::string value = s.second;
370  std::string padding(max_name_length - name.size(), ' ');
371 
372  if (value.empty()) {
373  continue;
374  }
375 
376  log.message(verbosity_threshold, " %s%s = \"%s\"\n",
377  name.c_str(), padding.c_str(), value.c_str());
378  }
379 
380  // Print double attributes:
381  for (const auto &d : doubles) {
382  std::string name = d.first;
383  std::vector<double> values = d.second;
384  std::string padding(max_name_length - name.size(), ' ');
385 
386  if (values.empty()) {
387  continue;
388  }
389 
390  const double large = 1.0e7;
391  const double small = 1.0e-4;
392  if ((std::fabs(values[0]) >= large) || (std::fabs(values[0]) <= small)) {
393  // use scientific notation if a number is big or small
394  log.message(verbosity_threshold, " %s%s = %12.3e\n",
395  name.c_str(), padding.c_str(), values[0]);
396  } else {
397  log.message(verbosity_threshold, " %s%s = %12.5f\n",
398  name.c_str(), padding.c_str(), values[0]);
399  }
400 
401  }
402 }
403 
404 ConstAttribute::ConstAttribute(const VariableMetadata *var, const std::string &name)
405  : m_name(name), m_var(const_cast<VariableMetadata*>(var)) {
406 }
407 
409  : m_name(std::move(a.m_name)), m_var(a.m_var) {
410  a.m_name.clear();
411  a.m_var = nullptr;
412 }
413 
414 ConstAttribute::operator std::string() const {
415  return m_var->get_string(m_name);
416 }
417 
418 ConstAttribute::operator double() const {
419  auto values = m_var->get_numbers(m_name);
420  if (values.size() == 1) {
421  return values[0];
422  }
424  "%s:%s has more than one value",
425  m_var->get_name().c_str(), m_name.c_str());
426 }
427 
428 ConstAttribute::operator std::vector<double> () const {
429  return m_var->get_numbers(m_name);
430 }
431 
432 void Attribute::operator=(const std::string &value) {
433  m_var->set_string(m_name, value);
434 }
435 void Attribute::operator=(const std::initializer_list<double> &value) {
436  m_var->set_numbers(m_name, value);
437 }
438 
439 void Attribute::operator=(const std::vector<double> &value) {
440  m_var->set_numbers(m_name, value);
441 }
442 
443 } // 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:50
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})
void set_numbers(const std::string &name, const std::vector< double > &values)
Set a scalar attribute to a single (scalar) value.
void report_to_stdout(const Logger &log, int verbosity_threshold) const
std::map< std::string, std::string > StringAttrs
unsigned int n_spatial_dimensions() const
void set_time_independent(bool flag)
VariableMetadata(const std::string &name, units::System::Ptr system, unsigned int ndims=0)
void set_units_without_validation(const std::string &value)
units::System::Ptr m_unit_system
The unit system to use.
void set_name(const std::string &name)
double get_number(const std::string &name) const
Get a single-valued scalar attribute.
const DoubleAttrs & all_doubles() const
void set_number(const std::string &name, double value)
Set a scalar attribute to a single (scalar) value.
const StringAttrs & all_strings() const
std::string get_string(const std::string &name) const
Get a string attribute.
void set_output_type(IO_Type type)
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.
std::map< std::string, std::vector< double > > DoubleAttrs
units::System::Ptr unit_system() const
std::map< std::string, std::vector< double > > m_doubles
scalar and array attributes
bool get_time_independent() const
std::vector< double > get_numbers(const std::string &name) const
Get an array-of-doubles attribute.
IO_Type get_output_type() const
bool has_attribute(const std::string &name) const
std::string get_name() const
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].
void set_string(const std::string &name, const std::string &value)
Set a string attribute.
std::shared_ptr< System > Ptr
Definition: Units.hh:47
#define PISM_ERROR_LOCATION
bool are_convertible(const Unit &u1, const Unit &u2)
Definition: Units.cc:244
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
IO_Type
Definition: IO_Flags.hh:31
@ PISM_NAT
Definition: IO_Flags.hh:32
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.