Merge pull request #113 from baagaard-usgs/feature-check-model-units

Add check for inconsistent model units

libsrc/geomodelgrids/serial/Query.cc

@@ -42,6 +42,19 @@ class geomodelgrids::serial::_Query {
  geomodelgrids::serial::Query::values_map_type createModelValuesIndex(const geomodelgrids::serial::Model& model,
  const std::vector<std::string>& queryNamesLower);
 
+ /** Check consistency of units in model.
+ *
+ * @param[in] valueUnits Map from index of query value to units.
+ * @param[in] modelValueIndex Map from index of query value to index of model value.
+ * @param[in] modelValues Names of model values.
+ * @param[in] modelUnitsLower Units of model values (lowercase).
+ */
+ static
+ void checkUnits(std::map<size_t, std::string>* valueUnits,
+ std::map<size_t, size_t>& modelValueIndex,
+ const std::vector<std::string>& modelValues,
+ const std::vector<std::string>& modelUnitsLower);
+
  static
  unsigned char tolower(unsigned char c);
 
@@ -89,14 +102,19 @@ geomodelgrids::serial::Query::initialize(const std::vector<std::string>& modelFi
  const size_t numModels = modelFilenames.size();
  _models.resize(numModels);
  _valuesIndex.resize(numModels);
- for (size_t i = 0; i < numModels; ++i) {
- _models[i] = new geomodelgrids::serial::Model();assert(_models[i]);
- _models[i]->setInputCRS(inputCRSString);
- _models[i]->open(modelFilenames[i].c_str(), geomodelgrids::serial::Model::READ);
- _models[i]->loadMetadata();
- _models[i]->initialize();
-
- _valuesIndex[i] = _Query::createModelValuesIndex(*_models[i], _valuesLowercase);
+ std::map<size_t, std::string> valueUnits;
+ for (size_t iModel = 0; iModel < numModels; ++iModel) {
+ _models[iModel] = new geomodelgrids::serial::Model();assert(_models[iModel]);
+ _models[iModel]->setInputCRS(inputCRSString);
+ _models[iModel]->open(modelFilenames[iModel].c_str(), geomodelgrids::serial::Model::READ);
+ _models[iModel]->loadMetadata();
+ _models[iModel]->initialize();
+
+ _valuesIndex[iModel] = _Query::createModelValuesIndex(*_models[iModel], _valuesLowercase);
+
+ const std::vector<std::string>& modelValues = _models[iModel]->getValueNames();
+ const std::vector<std::string>& modelUnitsLower = _Query::toLower(_models[iModel]->getValueUnits());
+ _Query::checkUnits(&valueUnits, _valuesIndex[iModel], modelValues, modelUnitsLower);
  } // for
 } // initialize
 
@@ -276,4 +294,27 @@ geomodelgrids::serial::_Query::createModelValuesIndex(const geomodelgrids::seria
 } // createModelValuesIndex
 
 
+// ------------------------------------------------------------------------------------------------
+void
+geomodelgrids::serial::_Query::checkUnits(std::map<size_t,std::string>* valueUnits,
+ std::map<size_t, size_t>& valuesIndex,
+ const std::vector<std::string>& modelValues,
+ const std::vector<std::string>& modelUnits) {
+ assert(valueUnits);
+
+ for (size_t iValue = 0; iValue < valuesIndex.size(); ++iValue) {
+ const size_t indexModel = valuesIndex[iValue];
+ if (valueUnits->count(iValue) == 0) {
+ (*valueUnits)[iValue] = modelUnits[indexModel];
+ } else if ((*valueUnits)[iValue] != modelUnits[indexModel]) {
+ std::ostringstream msg;
+ msg << "Inconsistent units among models for value '" << modelValues[indexModel] << "'. "
+ << "Inconsistent units are '" << (*valueUnits)[iValue] << "' and '"
+ << modelUnits[indexModel] << "'.";
+ throw std::runtime_error(msg.str());
+ } // if/else
+ } // for
+}
+
+
 // End of file

tests/data/Makefile.am

@@ -11,7 +11,8 @@ dist_noinst_DATA = \
  three-blocks-topo.h5 \
  one-block-topo-bad-topo.h5 \
  three-blocks-topo-bad-blocks.h5 \
- three-blocks-topo-missing-metadata.h5
+ three-blocks-topo-missing-metadata.h5 \
+ three-blocks-topo-inconsistent-units.h5
 
 
 # 'export' the input files by performing a mock install

tests/data/generate.py

@@ -419,6 +419,12 @@ def missing_metadata(self):
  for attr in top.attrs:
  del top.attrs[attr]
 
+ def inconsistent_units(self):
+ self.filename = "three-blocks-topo-inconsistent-units.h5"
+ self.create()
+ with h5py.File(self.filename, "a") as h5:
+ h5.attrs["data_units"] = ["km", "km/s"]
+
 
 # ==============================================================================
 if __name__ == "__main__":
@@ -429,6 +435,7 @@ def missing_metadata(self):
  OneBlockTopo().bad_topo_metadata()
  ThreeBlocksTopo().bad_block_metadata()
  ThreeBlocksTopo().missing_metadata()
+ ThreeBlocksTopo().inconsistent_units()
 
 
 # End of file

tests/libtests/apps/TestQuery.cc

@@ -45,6 +45,7 @@ class geomodelgrids::apps::TestQuery : public CppUnit::TestFixture {
  CPPUNIT_TEST(testRunTwoModels);
  CPPUNIT_TEST(testRunBadInput);
  CPPUNIT_TEST(testRunBadOutput);
+ CPPUNIT_TEST(testRunInconsistentUnits);
 
  CPPUNIT_TEST_SUITE_END();
 
@@ -111,6 +112,9 @@ class geomodelgrids::apps::TestQuery : public CppUnit::TestFixture {
  /// Test run() wth bad output.
  void testRunBadOutput(void);
 
+ /// Test run() wth inconsistent units.
+ void testRunInconsistentUnits(void);
+
 }; // class TestQuery
 CPPUNIT_TEST_SUITE_REGISTRATION(geomodelgrids::apps::TestQuery);
 
@@ -530,6 +534,25 @@ geomodelgrids::apps::TestQuery::testRunBadOutput(void) {
 } // testRunBadOutput
 
 
+// ------------------------------------------------------------------------------------------------
+// Test run() with inconsistent units.
+void
+geomodelgrids::apps::TestQuery::testRunInconsistentUnits(void) {
+ const int nargs = 6;
+ const char* const args[nargs] = {
+ "test",
+ "--models=../../data/one-block-flat.h5,../../data/three-blocks-topo-inconsistent-units.h5",
+ "--points=two-models.in",
+ "--output=two-models.out",
+ "--points-coordsys=EPSG:4326",
+ "--values=two,one",
+ };
+ Query query;
+ CPPUNIT_ASSERT_THROW(query.run(nargs, const_cast<char**>(args)), std::runtime_error);
+
+} // testRunInconsistentUnits
+
+
 // ------------------------------------------------------------------------------------------------
 void
 geomodelgrids::apps::_TestQuery::createPointsFile(std::ostream& sout,