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