tVariableMetadata.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
 (HTM) git clone git://src.adamsgaard.dk/pism
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       tVariableMetadata.cc (13987B)
       ---
            1 // Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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 
           22 #include "VariableMetadata.hh"
           23 #include "pism/util/io/File.hh"
           24 #include "pism_options.hh"
           25 #include "IceGrid.hh"
           26 
           27 #include "ConfigInterface.hh"
           28 #include "error_handling.hh"
           29 #include "pism/util/Logger.hh"
           30 #include "pism_utilities.hh"
           31 
           32 namespace pism {
           33 
           34 VariableMetadata::VariableMetadata(const std::string &name, units::System::Ptr system,
           35                                    unsigned int ndims)
           36   : m_n_spatial_dims(ndims),
           37     m_unit_system(system),
           38     m_short_name(name),
           39     m_time_independent(false),
           40     m_output_type(PISM_NAT) {
           41 
           42   clear_all_strings();
           43   clear_all_doubles();
           44 
           45   // long_name is unset
           46   // standard_name is unset
           47   // pism_intent is unset
           48   // coordinates is unset
           49 
           50   // valid_min and valid_max are unset
           51 }
           52 
           53 VariableMetadata::~VariableMetadata() {
           54   // empty
           55 }
           56 
           57 /** Get the number of spatial dimensions.
           58  */
           59 unsigned int VariableMetadata::get_n_spatial_dimensions() const {
           60   return m_n_spatial_dims;
           61 }
           62 
           63 void VariableMetadata::set_time_independent(bool flag) {
           64   m_time_independent = flag;
           65 }
           66 
           67 void VariableMetadata::set_output_type(IO_Type type) {
           68   m_output_type = type;
           69 }
           70 
           71 bool VariableMetadata::get_time_independent() const {
           72   return m_time_independent;
           73 }
           74 
           75 IO_Type VariableMetadata::get_output_type() const {
           76   return m_output_type;
           77 }
           78 
           79 units::System::Ptr VariableMetadata::unit_system() const {
           80   return m_unit_system;
           81 }
           82 
           83 //! @brief Check if the range `[min, max]` is a subset of `[valid_min, valid_max]`.
           84 /*! Throws an exception if this check failed.
           85  */
           86 void VariableMetadata::check_range(const std::string &filename, double min, double max) {
           87 
           88   auto units_string = get_string("units");
           89   auto name_string  = get_name();
           90   const char
           91     *units = units_string.c_str(),
           92     *name  = name_string.c_str(),
           93     *file  = filename.c_str();
           94 
           95   if (has_attribute("valid_min") and has_attribute("valid_max")) {
           96     double
           97       valid_min = get_number("valid_min"),
           98       valid_max = get_number("valid_max");
           99     if ((min < valid_min) or (max > valid_max)) {
          100       throw RuntimeError::formatted(PISM_ERROR_LOCATION, "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,
          103                                     valid_min, valid_max, units, min, units, max, units);
          104     }
          105   } else if (has_attribute("valid_min")) {
          106     double valid_min = get_number("valid_min");
          107     if (min < valid_min) {
          108       throw RuntimeError::formatted(PISM_ERROR_LOCATION, "some values of '%s' in '%s' are less than the valid minimum (%e %s).\n"
          109                                     "computed min = %e %s, computed max = %e %s",
          110                                     name, file,
          111                                     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) {
          116       throw RuntimeError::formatted(PISM_ERROR_LOCATION, "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,
          119                                     valid_max, units, min, units, max, units);
          120     }
          121   }
          122 }
          123 
          124 //! 3D version
          125 SpatialVariableMetadata::SpatialVariableMetadata(units::System::Ptr system, const std::string &name,
          126                                                  const std::vector<double> &zlevels)
          127   : VariableMetadata("unnamed", system),
          128     m_x("x", system),
          129     m_y("y", system),
          130     m_z("z", system) {
          131 
          132   init_internal(name, zlevels);
          133 }
          134 
          135 //! 2D version
          136 SpatialVariableMetadata::SpatialVariableMetadata(units::System::Ptr system, const std::string &name)
          137   : VariableMetadata("unnamed", system),
          138     m_x("x", system),
          139     m_y("y", system),
          140     m_z("z", system) {
          141 
          142   std::vector<double> z(1, 0.0);
          143   init_internal(name, z);
          144 }
          145 
          146 void SpatialVariableMetadata::init_internal(const std::string &name,
          147                                             const std::vector<double> &z_levels) {
          148   m_x.set_string("axis", "X");
          149   m_x.set_string("long_name", "X-coordinate in Cartesian system");
          150   m_x.set_string("standard_name", "projection_x_coordinate");
          151   m_x.set_string("units", "m");
          152 
          153   m_y.set_string("axis", "Y");
          154   m_y.set_string("long_name", "Y-coordinate in Cartesian system");
          155   m_y.set_string("standard_name", "projection_y_coordinate");
          156   m_y.set_string("units", "m");
          157 
          158   m_z.set_string("axis", "Z");
          159   m_z.set_string("long_name", "Z-coordinate in Cartesian system");
          160   m_z.set_string("units", "m");
          161   m_z.set_string("positive", "up");
          162 
          163   set_name(name);
          164 
          165   m_zlevels = z_levels;
          166 
          167   this->set_time_independent(false);
          168 
          169   if (m_zlevels.size() > 1) {
          170     get_z().set_name("z");      // default; can be overridden easily
          171     m_n_spatial_dims = 3;
          172   } else {
          173     get_z().set_name("");
          174     m_n_spatial_dims = 2;
          175   }
          176 }
          177 
          178 SpatialVariableMetadata::~SpatialVariableMetadata() {
          179   // empty
          180 }
          181 
          182 void SpatialVariableMetadata::set_levels(const std::vector<double> &levels) {
          183   if (levels.size() < 1) {
          184     throw RuntimeError(PISM_ERROR_LOCATION, "argument \"levels\" has to have length 1 or greater");
          185   }
          186   m_zlevels = levels;
          187 }
          188 
          189 
          190 const std::vector<double>& SpatialVariableMetadata::get_levels() const {
          191   return m_zlevels;
          192 }
          193 
          194 //! Report the range of a \b global Vec `v`.
          195 void VariableMetadata::report_range(const Logger &log, double min, double max,
          196                                     bool found_by_standard_name) {
          197 
          198   // units::Converter constructor will make sure that units are compatible.
          199   units::Converter c(m_unit_system,
          200                      this->get_string("units"),
          201                      this->get_string("glaciological_units"));
          202   min = c(min);
          203   max = c(max);
          204 
          205   std::string spacer(get_name().size(), ' ');
          206 
          207   if (has_attribute("standard_name")) {
          208 
          209     if (found_by_standard_name) {
          210       log.message(2, 
          211                   " %s / standard_name=%-10s\n"
          212                   "         %s \\ min,max = %9.3f,%9.3f (%s)\n",
          213                   get_name().c_str(),
          214                   get_string("standard_name").c_str(), spacer.c_str(), min, max,
          215                   get_string("glaciological_units").c_str());
          216     } else {
          217       log.message(2, 
          218                   " %s / WARNING! standard_name=%s is missing, found by short_name\n"
          219                   "         %s \\ min,max = %9.3f,%9.3f (%s)\n",
          220                   get_name().c_str(),
          221                   get_string("standard_name").c_str(), spacer.c_str(), min, max,
          222                   get_string("glaciological_units").c_str());
          223     }
          224 
          225   } else {
          226     log.message(2, 
          227                 " %s / %-10s\n"
          228                 "         %s \\ min,max = %9.3f,%9.3f (%s)\n",
          229                 get_name().c_str(),
          230                 get_string("long_name").c_str(), spacer.c_str(), min, max,
          231                 get_string("glaciological_units").c_str());
          232   }
          233 }
          234 
          235 VariableMetadata& SpatialVariableMetadata::get_x() {
          236   return m_x;
          237 }
          238 
          239 VariableMetadata& SpatialVariableMetadata::get_y() {
          240   return m_y;
          241 }
          242 
          243 VariableMetadata& SpatialVariableMetadata::get_z() {
          244   return m_z;
          245 }
          246 
          247 const VariableMetadata& SpatialVariableMetadata::get_x() const {
          248   return m_x;
          249 }
          250 
          251 const VariableMetadata& SpatialVariableMetadata::get_y() const {
          252   return m_y;
          253 }
          254 
          255 const VariableMetadata& SpatialVariableMetadata::get_z() const {
          256   return m_z;
          257 }
          258 
          259 //! Checks if an attribute is present. Ignores empty strings, except
          260 //! for the "units" attribute.
          261 bool VariableMetadata::has_attribute(const std::string &name) const {
          262 
          263   auto j = m_strings.find(name);
          264   if (j != m_strings.end()) {
          265     if (name != "units" && (j->second).empty()) {
          266       return false;
          267     }
          268 
          269     return true;
          270   }
          271 
          272   if (m_doubles.find(name) != m_doubles.end()) {
          273     return true;
          274   }
          275 
          276   return false;
          277 }
          278 
          279 bool VariableMetadata::has_attributes() const {
          280   return not (this->get_all_strings().empty() and this->get_all_doubles().empty());
          281 }
          282 
          283 void VariableMetadata::set_name(const std::string &name) {
          284   m_short_name = name;
          285 }
          286 
          287 //! Set a scalar attribute to a single (scalar) value.
          288 void VariableMetadata::set_number(const std::string &name, double value) {
          289   m_doubles[name] = std::vector<double>(1, value);
          290 }
          291 
          292 //! Set a scalar attribute to a single (scalar) value.
          293 void VariableMetadata::set_numbers(const std::string &name, const std::vector<double> &values) {
          294   m_doubles[name] = values;
          295 }
          296 
          297 void VariableMetadata::clear_all_doubles() {
          298   m_doubles.clear();
          299 }
          300 
          301 void VariableMetadata::clear_all_strings() {
          302   m_strings.clear();
          303 }
          304 
          305 std::string VariableMetadata::get_name() const {
          306   return m_short_name;
          307 }
          308 
          309 //! Get a single-valued scalar attribute.
          310 double VariableMetadata::get_number(const std::string &name) const {
          311   auto j = m_doubles.find(name);
          312   if (j != m_doubles.end()) {
          313     return (j->second)[0];
          314   } else {
          315     throw RuntimeError::formatted(PISM_ERROR_LOCATION, "variable \"%s\" does not have a double attribute \"%s\"",
          316                                   get_name().c_str(), name.c_str());
          317   }
          318 }
          319 
          320 //! Get an array-of-doubles attribute.
          321 std::vector<double> VariableMetadata::get_numbers(const std::string &name) const {
          322   auto j = m_doubles.find(name);
          323   if (j != m_doubles.end()) {
          324     return j->second;
          325   } else {
          326     return std::vector<double>();
          327   }
          328 }
          329 
          330 const VariableMetadata::StringAttrs& VariableMetadata::get_all_strings() const {
          331   return m_strings;
          332 }
          333 
          334 const VariableMetadata::DoubleAttrs& VariableMetadata::get_all_doubles() const {
          335   return m_doubles;
          336 }
          337 
          338 //! Set a string attribute.
          339 void VariableMetadata::set_string(const std::string &name, const std::string &value) {
          340 
          341   if (name == "units") {
          342     // create a dummy object to validate the units string
          343     units::Unit tmp(m_unit_system, value);
          344 
          345     m_strings[name] = value;
          346     m_strings["glaciological_units"] = value;
          347   } else if (name == "glaciological_units") {
          348     m_strings[name] = value;
          349 
          350     units::Unit internal(m_unit_system, get_string("units")),
          351       glaciological(m_unit_system, value);
          352 
          353     if (not units::are_convertible(internal, glaciological)) {
          354       throw RuntimeError::formatted(PISM_ERROR_LOCATION, "units \"%s\" and \"%s\" are not compatible",
          355                                     get_string("units").c_str(), value.c_str());
          356     }
          357   } else if (name == "short_name") {
          358     set_name(name);
          359   } else {
          360     m_strings[name] = value;
          361   }
          362 }
          363 
          364 //! Get a string attribute.
          365 /*!
          366  * Returns an empty string if an attribute is not set.
          367  */
          368 std::string VariableMetadata::get_string(const std::string &name) const {
          369   if (name == "short_name") {
          370      return get_name();
          371   }
          372 
          373   auto j = m_strings.find(name);
          374   if (j != m_strings.end()) {
          375     return j->second;
          376   } else {
          377     return std::string();
          378   }
          379 }
          380 
          381 void VariableMetadata::report_to_stdout(const Logger &log, int verbosity_threshold) const {
          382 
          383   const VariableMetadata::StringAttrs &strings = this->get_all_strings();
          384   const VariableMetadata::DoubleAttrs &doubles = this->get_all_doubles();
          385 
          386   // Find the maximum name length so that we can pad output below:
          387   size_t max_name_length = 0;
          388   for (auto s : strings) {
          389     max_name_length = std::max(max_name_length, s.first.size());
          390   }
          391   for (auto d : doubles) {
          392     max_name_length = std::max(max_name_length, d.first.size());
          393   }
          394 
          395   // Print text attributes:
          396   for (auto s : strings) {
          397     std::string name  = s.first;
          398     std::string value = s.second;
          399     std::string padding(max_name_length - name.size(), ' ');
          400 
          401     if (value.empty()) {
          402       continue;
          403     }
          404 
          405     log.message(verbosity_threshold, "  %s%s = \"%s\"\n",
          406                 name.c_str(), padding.c_str(), value.c_str());
          407   }
          408 
          409   // Print double attributes:
          410   for (auto d : doubles) {
          411     std::string name  = d.first;
          412     std::vector<double> values = d.second;
          413     std::string padding(max_name_length - name.size(), ' ');
          414 
          415     if (values.empty()) {
          416       continue;
          417     }
          418 
          419     if ((fabs(values[0]) >= 1.0e7) || (fabs(values[0]) <= 1.0e-4)) {
          420       // use scientific notation if a number is big or small
          421       log.message(verbosity_threshold, "  %s%s = %12.3e\n",
          422                   name.c_str(), padding.c_str(), values[0]);
          423     } else {
          424       log.message(verbosity_threshold, "  %s%s = %12.5f\n",
          425                   name.c_str(), padding.c_str(), values[0]);
          426     }
          427 
          428   }
          429 }
          430 
          431 bool set_contains(const std::set<std::string> &S, const VariableMetadata &variable) {
          432   return member(variable.get_name(), S);
          433 }
          434 
          435 TimeseriesMetadata::TimeseriesMetadata(const std::string &name, const std::string &dimension_name,
          436                            units::System::Ptr system)
          437   : VariableMetadata(name, system, 0) {
          438   m_dimension_name = dimension_name;
          439 }
          440 
          441 TimeseriesMetadata::~TimeseriesMetadata()
          442 {
          443   // empty
          444 }
          445 
          446 std::string TimeseriesMetadata::get_dimension_name() const {
          447   return m_dimension_name;
          448 }
          449 
          450 /// TimeBoundsMetadata
          451 
          452 TimeBoundsMetadata::TimeBoundsMetadata(const std::string &var_name, const std::string &dim_name,
          453                            units::System::Ptr system)
          454   : TimeseriesMetadata(var_name, dim_name, system) {
          455   m_bounds_name    = "nv";      // number of vertexes
          456 }
          457 
          458 TimeBoundsMetadata::~TimeBoundsMetadata() {
          459   // empty
          460 }
          461 
          462 std::string TimeBoundsMetadata::get_bounds_name() const {
          463   return m_bounds_name;
          464 }
          465 
          466 } // end of namespace pism