Merge pull request #126 from baagaard-usgs/feature-info-bbox-units

Add bounding box and coordinate system units to Info app

docs/user/apps/info.md

@@ -50,6 +50,7 @@ geomodelgrids_info --all --models=tests/data/three-blocks-topo.h5
 Model: tests/data/three-blocks-topo.h5
  Verification
  Verifying metadata...OK
+ Verifying model coordinate system ...OK
  Verifying surface 'top surface'...OK
  Verifying surface 'topography/bathymetry'...OK
  Verifying resolution of top surface matches resolution of topography/bathymetry...OK
@@ -73,14 +74,17 @@ Model: tests/data/three-blocks-topo.h5
  DOI: this.is.a.doi
  Version: 1.0.0
  License: CC0
- Dimensions of model (m): x=60000, y=120000, z=45000
+ Dimensions of model: x=60000, y=120000, z=45000
+ Bounding box (WGS84): (34.3954, -117.8241) (34.6535, -117.2496) (35.603, -117.8789) (35.342, -118.4583)
  Coordinate system:
  CRS (PROJ, EPSG, WKT): EPSG:3311
+ Coordinate system units: x=meter, y=meter, z=meter (assumed)
  Origin: x=200000, y=-400000
  Azimuth (degrees) of y axis from north: 330
  Values stored in model:
  0: one (m)
  1: two (m/s)
+ Vertex-based data
  Surfaces
  Top surface:
  Number of points: x=13, y=25
@@ -91,17 +95,16 @@ Model: tests/data/three-blocks-topo.h5
  Blocks (3)
  Block 'top'
  Resolution: x=10000, y=10000, z=5000
- Elevation (m) of top of block in logical space: 0
+ Elevation of top of block in logical space: 0
  Number of points: x=7, y=13, z=2
- Dimensions (m): x=60000, y=120000, z=5000
+ Dimensions: x=60000, y=120000, z=5000
  Block 'middle'
  Resolution: x=20000, y=20000, z=10000
- Elevation (m) of top of block in logical space: -5000
+ Elevation of top of block in logical space: -5000
  Number of points: x=4, y=7, z=3
- Dimensions (m): x=60000, y=120000, z=20000
+ Dimensions: x=60000, y=120000, z=20000
  Block 'bottom'
  Resolution: x=30000, y=30000, z=10000
- Elevation (m) of top of block in logical space: -25000
+ Elevation of top of block in logical space: -25000
  Number of points: x=3, y=5, z=3
- Dimensions (m): x=60000, y=120000, z=20000
-``` 
+ Dimensions: x=60000, y=120000, z=20000``` 

libsrc/geomodelgrids/apps/Info.cc

@@ -6,6 +6,7 @@
 #include "geomodelgrids/serial/ModelInfo.hh" // USES ModelInfo
 #include "geomodelgrids/serial/Block.hh" // USES Block
 #include "geomodelgrids/serial/Surface.hh" // USES Surface
+#include "geomodelgrids/utils/CRSTransformer.hh" // USES CRSTransformer
 
 #include <cmath> // USES fabs()
 
@@ -25,6 +26,8 @@ namespace geomodelgrids {
  std::string indent(const size_t level,
  const size_t width=4);
 
+ void verifyCoordSys(const geomodelgrids::serial::Model* model);
+
  void verifySurface(const geomodelgrids::serial::Surface* topography,
  const geomodelgrids::serial::Model* model,
  const char* const name);
@@ -35,6 +38,8 @@ namespace geomodelgrids {
  void verifyBlocksZ(const std::vector<geomodelgrids::serial::Block*>& blocks,
  const double zBottom);
 
+ const double* computeBoundingBox(const geomodelgrids::serial::Model* model);
+
  } // _Info
  } // apps
 } // geomodelgrids
@@ -236,7 +241,21 @@ geomodelgrids::apps::Info::_printDescription(geomodelgrids::serial::Model* const
  } // if
 
  const double* dims = model->getDims();
- std::cout << _Info::indent(1) << "Dimensions of model (m): x=" << dims[0] << ", y="<< dims[1] << ", z=" << dims[2] << "\n";
+ std::cout << _Info::indent(1) << "Dimensions of model: x=" << dims[0] << ", y="<< dims[1] << ", z=" << dims[2] << "\n";
+
+ const double* bbox = _Info::computeBoundingBox(model);
+ std::cout << _Info::indent(1) << "Bounding box (WGS84):" << std::resetiosflags(std::ios::fixed);
+ for (size_t i = 0; i < 4; ++i) {
+ const size_t dim = 2;
+ std::cout << " ("
+ << std::setprecision(6)
+ << bbox[i*dim+0]
+ << ", "
+ << std::setprecision(7)
+ << bbox[i*dim+1]
+ << ")";
+ } // for
+ std::cout << "\n";
 } // _printDescription
 
 
@@ -250,6 +269,13 @@ geomodelgrids::apps::Info::_printCoordSys(geomodelgrids::serial::Model* const mo
 
  std::cout << _Info::indent(2) << "CRS (PROJ, EPSG, WKT): " << model->getCRSString() << "\n";
 
+ std::string xUnit;
+ std::string yUnit;
+ std::string zUnit;
+ geomodelgrids::utils::CRSTransformer::getCRSUnits(&xUnit, &yUnit, &zUnit, model->getCRSString().c_str());
+ std::cout << _Info::indent(2) << "Coordinate system units:"
+ << " x=" << xUnit << ", y=" << yUnit << ", z=" << zUnit << "\n";
+
  const double* origin = model->getOrigin();
  std::cout << _Info::indent(2) << "Origin: x=" << origin[0] <<", y=" << origin[1] << "\n";
 
@@ -274,13 +300,13 @@ geomodelgrids::apps::Info::_printValues(geomodelgrids::serial::Model* const mode
  } // for
 
  switch (model->getDataLayout()) {
-  case geomodelgrids::serial::Model::VERTEX:
+ case geomodelgrids::serial::Model::VERTEX:
  std::cout << _Info::indent(2) << "Vertex-based data\n";
- break; 
-  case geomodelgrids::serial::Model::CELL:
+ break;
+ case geomodelgrids::serial::Model::CELL:
  std::cout << _Info::indent(2) << "Cell-based data\n";
- break; 
-  default:
+ break;
+ default:
  std::cout << _Info::indent(2) << "Unknown data layout\n";
  }
 } // _printValues
@@ -327,15 +353,15 @@ geomodelgrids::apps::Info::_printBlocks(geomodelgrids::serial::Model* const mode
  << "x=" << blocks[i]->getResolutionX()
  << ", y=" << blocks[i]->getResolutionY()
  << ", z=" << blocks[i]->getResolutionZ() << "\n";
- std::cout << _Info::indent(3) << "Elevation (m) of top of block in logical space: " << blocks[i]->getZTop() << "\n";
+ std::cout << _Info::indent(3) << "Elevation of top of block in logical space: " << blocks[i]->getZTop() << "\n";
 
  const size_t* dims = blocks[i]->getDims();
  std::cout << _Info::indent(3) << "Number of points: x=" << dims[0] << ", y=" << dims[1] << ", z=" << dims[2] << "\n";
 
  const double dim_x = blocks[i]->getResolutionX() * (dims[0] - 1);
  const double dim_y = blocks[i]->getResolutionY() * (dims[1] - 1);
  const double dim_z = blocks[i]->getResolutionZ() * (dims[2] - 1);
- std::cout << _Info::indent(3) << "Dimensions (m): x=" << dim_x << ", y=" << dim_y << ", z=" << dim_z << "\n";
+ std::cout << _Info::indent(3) << "Dimensions: x=" << dim_x << ", y=" << dim_y << ", z=" << dim_z << "\n";
  } // for
 } // _printBlocks
 
@@ -357,6 +383,7 @@ geomodelgrids::apps::Info::_verify(geomodelgrids::serial::Model* const model) {
  ok = false;
  } // try/catch
  if (ok) { std::cout << "OK\n"; }
+ _Info::verifyCoordSys(model);
 
  const geomodelgrids::serial::Surface* surfaceTop = model->getTopSurface();
  if (surfaceTop) { _Info::verifySurface(surfaceTop, model, "top surface"); }
@@ -419,6 +446,28 @@ geomodelgrids::apps::_Info::indent(const size_t level,
 } // indent
 
 
+// ------------------------------------------------------------------------------------------------
+void
+geomodelgrids::apps::_Info::verifyCoordSys(const geomodelgrids::serial::Model* model) {
+ std::cout << _Info::indent(2) << "Verifying model coordinate system ...";
+
+ bool ok = true;
+
+ std::string xUnit;
+ std::string yUnit;
+ geomodelgrids::utils::CRSTransformer::getCRSUnits(&xUnit, &yUnit, NULL,
+ model->getCRSString().c_str());
+ if (xUnit != yUnit) {
+ if (ok) { std::cout << "FAIL\n"; }
+ std::cout << _Info::indent(3) << "Units for x (" << xUnit << ") and y (" << yUnit
+ << ") in model coordinate system do not match.\n";
+ ok = false;
+ } // if
+
+ if (ok) { std::cout << "OK\n"; }
+}
+
+
 // ------------------------------------------------------------------------------------------------
 void
 geomodelgrids::apps::_Info::verifySurface(const geomodelgrids::serial::Surface* surface,
@@ -552,4 +601,43 @@ geomodelgrids::apps::_Info::verifyBlocksZ(const std::vector<geomodelgrids::seria
 }
 
 
+// ------------------------------------------------------------------------------------------------
+const double*
+geomodelgrids::apps::_Info::computeBoundingBox(const geomodelgrids::serial::Model* model) {
+ assert(model);
+
+ const double* dims = model->getDims();assert(dims);
+ const double* origin = model->getOrigin();assert(origin);
+ const double yazimuth = model->getYAzimuth();
+ const std::string& crsString = model->getCRSString();
+
+ const size_t npts = 4;
+ const size_t dim = 2;
+ const double thetaR = (360.0 - yazimuth) / 180.0 * M_PI;
+ const double cosaz = cos(thetaR);
+ const double sinaz = sin(thetaR);
+
+ static double bbox_model[npts*dim];
+ bbox_model[0*dim+0] = origin[0];
+ bbox_model[0*dim+1] = origin[1];
+ bbox_model[1*dim+0] = origin[0] + dims[0]*cosaz;
+ bbox_model[1*dim+1] = origin[1] + dims[0]*sinaz;
+ bbox_model[2*dim+0] = origin[0] + dims[0]*cosaz - dims[1]*sinaz;
+ bbox_model[2*dim+1] = origin[1] + dims[0]*sinaz + dims[1]*cosaz;
+ bbox_model[3*dim+0] = origin[0] - dims[1]*sinaz;
+ bbox_model[3*dim+1] = origin[1] + dims[1]*cosaz;
+
+ static double bbox_geo[npts*dim];
+ geomodelgrids::utils::CRSTransformer transformer;
+ transformer.setSrc(crsString.c_str());
+ transformer.setDest("EPSG:4326");
+ transformer.initialize();
+ for (size_t i = 0; i < npts; ++i) {
+ transformer.transform(&bbox_geo[i*dim+0], &bbox_geo[i*dim+1], NULL, bbox_model[i*dim+0], bbox_model[i*dim+1], 0.0);
+ } // for
+
+ return bbox_geo;
+}
+
+
 // End of file

libsrc/geomodelgrids/utils/CRSTransformer.cc

@@ -4,6 +4,29 @@
 
 #include <sstream> // USES std::ostringstream
 #include <cassert> // USES assert()
+#include <cstring> // USES strlen()
+#include <strings.h> // USES stdcasecmp()
+
+namespace geomodelgrids {
+ namespace utils {
+ class _CRSTransformer {
+public:
+
+ static
+ void getUnits(std::string* xUnit,
+ std::string* yUnit,
+ std::string* zUnit,
+ const PJ* projCoordSys);
+
+ static
+ void getAxisUnit(std::string* unit,
+ const PJ* projCoordSys,
+ const int axisIndex);
+
+ };
+
+ }
+}
 
 // ------------------------------------------------------------------------------------------------
 // Default constructor.
@@ -136,4 +159,87 @@ geomodelgrids::utils::CRSTransformer::createGeoToXYAxisOrder(const char* crsStri
 }
 
 
+// ------------------------------------------------------------------------------------------------
+// Get units for CRS.
+void
+geomodelgrids::utils::CRSTransformer::getCRSUnits(std::string* xUnit,
+ std::string* yUnit,
+ std::string* zUnit,
+ const char* crsString) {
+ if (xUnit) { *xUnit = "unknown"; }
+ if (yUnit) { *yUnit = "unknown"; }
+ if (zUnit) { *zUnit = "meter (assumed)"; }
+ if (!crsString || (0 == strlen(crsString))) { return; }
+
+ PJ* proj = proj_create(NULL, crsString);assert(proj);
+ PJ* projCoordSys = proj_crs_get_coordinate_system(NULL, proj);
+ if (projCoordSys) {
+ _CRSTransformer::getUnits(xUnit, yUnit, zUnit, projCoordSys);
+ proj_destroy(proj);
+ proj_destroy(projCoordSys);
+ return;
+ } else if (proj_get_type(proj) == PJ_TYPE_BOUND_CRS) {
+ PJ* projTmp = proj_get_source_crs(NULL, proj);
+ proj_destroy(proj);
+ if (!projTmp) {
+ return;
+ } // if
+ const char* srcWKT = proj_as_wkt(NULL, projTmp, PJ_WKT2_2019, NULL);
+ PJ* projSrc = proj_create(NULL, srcWKT);
+ projCoordSys = proj_crs_get_coordinate_system(NULL, projSrc);
+ proj_destroy(projTmp);
+ proj_destroy(projSrc);
+ if (!projCoordSys) {
+ return;
+ }
+ _CRSTransformer::getUnits(xUnit, yUnit, zUnit, projCoordSys);
+ proj_destroy(projCoordSys);
+ return;
+ }
+}
+
+
+// ------------------------------------------------------------------------------------------------
+// Get units for CRS.
+void
+geomodelgrids::utils::_CRSTransformer::getUnits(std::string* xUnit,
+ std::string* yUnit,
+ std::string* zUnit,
+ const PJ* projCoordSys) {
+ const int numAxes = proj_cs_get_axis_count(NULL, projCoordSys);
+ if (xUnit && (numAxes >= 1)) {
+ _CRSTransformer::getAxisUnit(xUnit, projCoordSys, 0);
+ }
+ if (yUnit && (numAxes >= 2)) {
+ _CRSTransformer::getAxisUnit(yUnit, projCoordSys, 1);
+ }
+ if (zUnit && (numAxes >= 3)) {
+ _CRSTransformer::getAxisUnit(zUnit, projCoordSys, 2);
+ }
+}
+
+
+// ------------------------------------------------------------------------------------------------
+// Get units for CRS.
+void
+geomodelgrids::utils::_CRSTransformer::getAxisUnit(std::string* unit,
+ const PJ* projCoordSys,
+ const int axisIndex) {
+ assert(unit);
+ assert(projCoordSys);
+
+ const char* buffer = NULL;
+ proj_cs_get_axis_info(NULL, projCoordSys, 0, NULL, NULL, NULL, NULL, &buffer, NULL, NULL);
+ if (strcasecmp("metre", buffer) == 0) {
+ *unit = "meter";
+ } else if (strcasecmp("meter", buffer) == 0) {
+ *unit = "meter";
+ } else if (strcasecmp("kilometre", buffer) == 0) {
+ *unit = "kilometer";
+ } else {
+ *unit = buffer;
+ }
+}
+
+
 // End of file

libsrc/geomodelgrids/utils/CRSTransformer.hh

@@ -81,6 +81,19 @@ public:
  static
  CRSTransformer* createGeoToXYAxisOrder(const char* crsString);
 
+ /** Get units for CRS.
+ *
+ * @param[out] xUnit Units for x axis (can be NULL).
+ * @param[out] yUnit Units for y axis (can be NULL).
+ * @param[out] zUnit Units for z axis (can be NULL).
+ * @param[in] crsString CRS for coordinate system.
+ */
+ static
+ void getCRSUnits(std::string* xUnit,
+ std::string* yUnit,
+ std::string* zUnit,
+ const char* crsString);
+
  // PRIVATE MEMBERS ----------------------------------------------------------------------------
 private: