[cig-commits] r8301 - in cs/spatialdata-0.1/trunk/libsrc: .
spatialdb
brad at geodynamics.org
brad at geodynamics.org
Sun Nov 18 15:59:15 PST 2007
Author: brad
Date: 2007-11-18 15:59:14 -0800 (Sun, 18 Nov 2007)
New Revision: 8301
Added:
cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc
cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.hh
cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc
cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.hh
Modified:
cs/spatialdata-0.1/trunk/libsrc/Makefile.am
cs/spatialdata-0.1/trunk/libsrc/spatialdb/Makefile.am
Log:
Initial implementation of spatial database interface to SCEC CVM-H.
Modified: cs/spatialdata-0.1/trunk/libsrc/Makefile.am
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/Makefile.am 2007-11-17 01:30:59 UTC (rev 8300)
+++ cs/spatialdata-0.1/trunk/libsrc/Makefile.am 2007-11-18 23:59:14 UTC (rev 8301)
@@ -29,7 +29,9 @@
geocoords/CSPicklerAscii.cc \
geocoords/Geoid.cc \
geocoords/Projector.cc \
+ spatialdb/GocadVoxet.cc \
spatialdb/GravityField.cc \
+ spatialdb/SCECCVMH.cc \
spatialdb/SimpleDB.cc \
spatialdb/SpatialDB.cc \
spatialdb/SimpleDBQuery.cc \
Added: cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc 2007-11-17 01:30:59 UTC (rev 8300)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc 2007-11-18 23:59:14 UTC (rev 8301)
@@ -0,0 +1,339 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+//#include <portinfo>
+
+#include "GocadVoxet.hh" // Implementation of class methods
+
+#include "spatialdata/utils/LineParser.hh" // USES LineParser
+
+#include <fstream> // USES std::ifstream
+#include <math.h> // USES round()
+
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringsgream
+#include <assert.h> // USES assert()
+
+#include <iostream> // TEMPORARY
+
+#if defined(WORDS_BIGENDIAN)
+#define NATIVE_BIG_ENDIAN
+#else
+#define NATIVE_LITTLE_ENDIAN
+#endif
+
+// ----------------------------------------------------------------------
+// Constructor
+spatialdata::spatialdb::GocadVoxet::GocadVoxet(void) :
+ _data(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+spatialdata::spatialdb::GocadVoxet::~GocadVoxet(void)
+{ // destructor
+ delete[] _data; _data = 0;
+} // destructor
+
+// ----------------------------------------------------------------------
+// Read data files.
+void
+spatialdata::spatialdb::GocadVoxet::read(const char* dir,
+ const char* filename,
+ const char* property)
+{ // read
+ std::string fullname;
+ fullname = std::string(dir) + std::string("/") + std::string(filename);
+ _readVoxetFile(fullname.c_str(), property);
+
+ fullname = std::string(dir) + std::string("/")
+ + std::string(_property.filename);
+ _readPropertyFile(fullname.c_str());
+} // read
+
+// ----------------------------------------------------------------------
+// Query voxet for value.
+int
+spatialdata::spatialdb::GocadVoxet::query(double* value,
+ const double pt[3]) const
+{ // query
+ // Compute indices of voxet containing pt
+ const int numX = _geometry.n[0];
+ const int numY = _geometry.n[1];
+ const int numZ = _geometry.n[2];
+ const int indexX =
+ round( (pt[0] - (_geometry.o[0]+_geometry.min[0]))*_geometry.scale[0]);
+ const int indexY =
+ round( (pt[1] - (_geometry.o[1]+_geometry.min[1]))*_geometry.scale[1]);
+ const int indexZ =
+ round( (pt[2] - (_geometry.o[2]+_geometry.min[2]))*_geometry.scale[2]);
+
+ // Check if voxet is in range
+ int flag = 0;
+ if (indexX >= 0 && indexX < numX &&
+ indexY >= 0 && indexY < numY &&
+ indexZ >= 0 && indexZ < numZ) {
+ int indexV = indexZ*numY*numX + indexY*numX + indexX;
+
+#if 0
+ std::cout << "pt[0]: " << pt[0]
+ << ", pt[1]: " << pt[1]
+ << ", pt[2]: " << pt[2]
+ << ", numX: " << numX
+ << ", numY: " << numY
+ << ", numZ: " << numZ
+ << ", scaleX: " << _geometry.scale[0]
+ << ", scaleY: " << _geometry.scale[1]
+ << ", scaleZ: " << _geometry.scale[2]
+ << ", indexX: " << indexX
+ << ", indexY: " << indexY
+ << ", indexZ: " << indexZ
+ << ", indexV: " << indexV
+ << ", byteV: " << indexV*sizeof(float)
+ << std::endl;
+#endif
+ *value = _data[indexV];
+ } else {
+ *value = _property.noDataValue;
+ flag = 1;
+ } // if/else
+
+ return flag;
+} // query
+
+// ----------------------------------------------------------------------
+// Read voxet file.
+void
+spatialdata::spatialdb::GocadVoxet::_readVoxetFile(const char* filename,
+ const char* propertyName)
+{ // _readVoxetFile
+ try {
+ std::ifstream vfile(filename);
+ if (!vfile.is_open() || !vfile.good()) {
+ std::ostringstream msg;
+ msg << "Could not open Gocad Voxet file '" << filename
+ << "' for reading.\n";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ utils::LineParser parser(vfile);
+ parser.eatwhitespace(true);
+
+ std::istringstream buffer;
+
+ buffer.str(parser.next());
+ buffer.clear();
+
+ const char* headerE = "GOCAD Voxet 1";
+ const int headerLen = strlen(headerE);
+ std::string hbuffer;
+ hbuffer.resize(headerLen+1);
+ buffer.read((char*) hbuffer.c_str(), sizeof(char)*headerLen);
+ hbuffer[headerLen] = '\0';
+ if (0 != strcmp(headerE, hbuffer.c_str())) {
+ std::ostringstream msg;
+ msg
+ << "Header '" << hbuffer << "' does not match expected header '"
+ << headerE << "' in Gocad Voxet file '" << filename << "'.\n";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ std::string token;
+ const int maxIgnore = 256;
+
+ _resetData();
+ int propertyId = 0;
+ int propertyIdTarget = -1;
+
+ buffer.str(parser.next());
+ buffer.clear();
+ buffer >> token;
+ while (vfile.good() && !vfile.eof()) {
+ if (0 == strcmp(token.c_str(), "AXIS_O"))
+ for (int i=0; i < 3; ++i)
+ buffer >> _geometry.o[i];
+ else if (0 == strcmp(token.c_str(), "AXIS_U"))
+ for (int i=0; i < 3; ++i)
+ buffer >> _geometry.u[i];
+ else if (0 == strcmp(token.c_str(), "AXIS_V"))
+ for (int i=0; i < 3; ++i)
+ buffer >> _geometry.v[i];
+ else if (0 == strcmp(token.c_str(), "AXIS_W"))
+ for (int i=0; i < 3; ++i)
+ buffer >> _geometry.w[i];
+ else if (0 == strcmp(token.c_str(), "AXIS_MIN"))
+ for (int i=0; i < 3; ++i)
+ buffer >> _geometry.min[i];
+ else if (0 == strcmp(token.c_str(), "AXIS_MAX"))
+ for (int i=0; i < 3; ++i)
+ buffer >> _geometry.max[i];
+ else if (0 == strcmp(token.c_str(), "AXIS_N"))
+ for (int i=0; i < 3; ++i)
+ buffer >> _geometry.n[i];
+ else if (0 == strcmp(token.c_str(), "AXIS_NAME"))
+ for (int i=0; i < 3; ++i)
+ buffer >> _geometry.name[i];
+ else if (0 == strcmp(token.c_str(), "AXIS_TYPE"))
+ for (int i=0; i < 3; ++i)
+ buffer >> _geometry.type[i];
+
+ else if (0 == strcmp(token.c_str(), "PROPERTY")) {
+ std::string name;
+ buffer >> propertyId;
+ buffer >> name;
+ if (0 == strcmp(name.c_str(), propertyName)) {
+ propertyIdTarget = propertyId;
+ _property.name = name;
+ } // if
+ } else if (0 == strcmp(token.c_str(), "PROP_NO_DATA_VALUE")) {
+ buffer >> propertyId;
+ if (propertyIdTarget == propertyId)
+ buffer >> _property.noDataValue;
+ } else if (0 == strcmp(token.c_str(), "PROP_ESIZE")) {
+ buffer >> propertyId;
+ if (propertyIdTarget == propertyId)
+ buffer >> _property.esize;
+ } else if (0 == strcmp(token.c_str(), "PROP_ETYPE")) {
+ buffer >> propertyId;
+ if (propertyIdTarget == propertyId)
+ buffer >> _property.etype;
+ } else if (0 == strcmp(token.c_str(), "PROP_OFFSET")) {
+ buffer >> propertyId;
+ if (propertyIdTarget == propertyId)
+ buffer >> _property.offset;
+ } else if (0 == strcmp(token.c_str(), "PROP_FILE")) {
+ buffer >> propertyId;
+ if (propertyIdTarget == propertyId)
+ buffer >> _property.filename;
+ } else if (0 == strcmp(token.c_str(), "END"))
+ break;
+
+ buffer.str(parser.next());
+ buffer.clear();
+ buffer >> token;
+ } // while
+ if (token != "END" || !vfile.good())
+ throw std::runtime_error("I/O error while parsing Gocad Voxet tokens.");
+ } catch (const std::exception& err) {
+ std::ostringstream msg;
+ msg << "Error occurred while reading Gocad Voxet file '"
+ << filename << "'.\n"
+ << err.what();
+ throw std::runtime_error(msg.str());
+ } catch (...) {
+ std::ostringstream msg;
+ msg << "Unknown error occurred while reading Gocad Voxet file '"
+ << filename << "'.\n";
+ throw std::runtime_error(msg.str());
+ } // try/catch
+
+ // Compute resolution
+ assert(0 == _geometry.u[1] && 0 == _geometry.u[2]);
+ assert(0 == _geometry.v[0] && 0 == _geometry.v[2]);
+ assert(0 == _geometry.w[0] && 0 == _geometry.w[1]);
+ const double lenX = _geometry.max[0]*_geometry.u[0] - _geometry.min[0];
+ const double lenY = _geometry.max[1]*_geometry.v[1] - _geometry.min[1];
+ const double lenZ = _geometry.max[2]*_geometry.w[2] - _geometry.min[2];
+ _geometry.scale[0] = (lenX > 0.0) ? (_geometry.n[0]-1) / lenX : 1.0e+30;
+ _geometry.scale[1] = (lenY > 0.0) ? (_geometry.n[1]-1) / lenY : 1.0e+30;
+ _geometry.scale[2] = (lenZ > 0.0) ? (_geometry.n[2]-1) / lenZ : 1.0e+30;
+} // _readVoxetFile
+
+// ----------------------------------------------------------------------
+// Read property data.
+void
+spatialdata::spatialdb::GocadVoxet::_readPropertyFile(const char* filename)
+{ // _readPropertyFile
+ assert(sizeof(float) == _property.esize);
+ const int nvals = _geometry.n[0] * _geometry.n[1] * _geometry.n[2];
+ delete[] _data; _data = new float[nvals];
+
+ try {
+ std::ifstream pfile(filename);
+ if (!pfile.is_open() || !pfile.good()) {
+ std::ostringstream msg;
+ msg << "Could not open Gocad Voxet property file '" << filename
+ << "' for reading.\n";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ pfile.read((char*) _data, sizeof(float)*nvals);
+ _endianBigToNative(_data, nvals);
+ } catch (const std::exception& err) {
+ std::ostringstream msg;
+ msg << "Error occurred while reading Gocad Voxet property file '"
+ << filename << "'.\n"
+ << err.what();
+ throw std::runtime_error(msg.str());
+ } catch (...) {
+ std::ostringstream msg;
+ msg << "Unknown error occurred while reading Gocad Voxet property file '"
+ << filename << "'.\n";
+ throw std::runtime_error(msg.str());
+ } // try/catch
+} // _readPropertyFile
+
+// ----------------------------------------------------------------------
+// Convert array of float values from big endian to native float type.
+void
+spatialdata::spatialdb::GocadVoxet::_endianBigToNative(float* vals,
+ const int nvals) const
+{ // _endianBigToNative
+ assert(0 != vals);
+#if defined(NATIVE_LITTLE_ENDIAN)
+ for (int i=0; i < nvals; ++i) {
+ char* buf = (char*) (vals + i);
+ char tmp = buf[3];
+ buf[3] = buf[0];
+ buf[0] = tmp;
+ tmp = buf[2];
+ buf[2] = buf[1];
+ buf[1] = tmp;
+ } // for
+#endif
+} // _endianBigToNative
+
+// ----------------------------------------------------------------------
+// Reset class data.
+void
+spatialdata::spatialdb::GocadVoxet::_resetData(void)
+{ // _resetData
+ const float zero[] = { 0.0, 0.0, 0.0 };
+ const float one[] = { 1.0, 1.0, 1.0 };
+
+ memcpy(_geometry.o, zero, 3*sizeof(float));
+ memcpy(_geometry.u, zero, 3*sizeof(float));
+ memcpy(_geometry.v, zero, 3*sizeof(float));
+ memcpy(_geometry.w, zero, 3*sizeof(float));
+ memcpy(_geometry.min, zero, 3*sizeof(float));
+ memcpy(_geometry.max, one, 3*sizeof(float));
+ _geometry.name[0] = "";
+ _geometry.name[1] = "";
+ _geometry.name[2] = "";
+ _geometry.type[0] = "";
+ _geometry.type[1] = "";
+ _geometry.type[2] = "";
+
+ _property.name = "";
+ _property.noDataValue = -99999;
+ _property.esize = 4;
+ _property.isSigned = 1;
+ _property.etype = "IEEE";
+ _property.offset = 0;
+ _property.filename = "";
+
+ delete[] _data; _data = 0;
+} // _resetData
+
+
+// End of file
Added: cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.hh
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.hh 2007-11-17 01:30:59 UTC (rev 8300)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.hh 2007-11-18 23:59:14 UTC (rev 8301)
@@ -0,0 +1,135 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/** @file libsrc/spatialdb/GocadVoxet.hh
+ *
+ * @brief C++ implementation of object for managing Gocad Voxets.
+ *
+ */
+
+#if !defined(spatialdata_spatialdb_gocadvoxet_hh)
+#define spatialdata_spatialdb_gocadvoxet_hh
+
+#include <string> // USES std::string
+
+namespace spatialdata {
+ namespace spatialdb {
+ class GocadVoxet;
+ } // spatialdb
+} // spatialdata
+
+class spatialdata::spatialdb::GocadVoxet
+{ // GocadVoxet
+ friend class TestGocadVOxet; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ GocadVoxet(void);
+
+ /// Destructor
+ ~GocadVoxet(void);
+
+ /** Read voxet file and property data.
+ *
+ * @param dir Directory containing voxet data files.
+ * @param filename Name of voxet file.
+ * @param property Name of property in voxet file to read.
+ */
+ void read(const char* dir,
+ const char* filename,
+ const char* property);
+
+ /** Query voxet for value.
+ *
+ * @param value Value for result.
+ * @param pt Location of query.
+ * @returns 0 if pt is inside voxet, 1 if outside voxet.
+ */
+ int query(double* value,
+ const double pt[3]) const;
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ GocadVoxet(const GocadVoxet&); ///< Not implemented
+ const GocadVoxet& operator=(const GocadVoxet&); ///< Not implemented
+
+// PRIVATE STRUCTS //////////////////////////////////////////////////////
+private :
+
+ struct Geometry {
+ float o[3];
+ float u[3];
+ float v[3];
+ float w[3];
+ float min[3];
+ float max[3];
+ int n[3];
+ std::string name[3];
+ std::string type[3];
+ float scale[3];
+ }; // axis
+
+ struct Property {
+ std::string name;
+ float noDataValue;
+ int esize;
+ int isSigned;
+ std::string etype;
+ int offset;
+ std::string filename;
+ }; // property
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /** Read voxet file.
+ *
+ * @param filename Name of Gocad voxet file.
+ * @param propertyName Name of property to read.
+ */
+ void _readVoxetFile(const char* filename,
+ const char* propertyName);
+
+ /** Read property file.
+ *
+ * @param filename Name of voxet data file.
+ */
+ void _readPropertyFile(const char* filename);
+
+ /** Convert array of float values from big endian to native float type.
+ *
+ * @param vals Array of values.
+ * @param nvals Number of values.
+ */
+ void
+ _endianBigToNative(float* vals,
+ const int nvals) const;
+
+ /// Reset class data.
+ void _resetData(void);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ Geometry _geometry; ///< Geometry of voxet data.
+ Property _property; ///< Voxet properties.
+ float* _data; ///< Array with data values
+
+}; // GocadVoxet
+
+#endif // spatialdata_spatialdb_gocadvoxet_hh
+
+
+// End of file
Modified: cs/spatialdata-0.1/trunk/libsrc/spatialdb/Makefile.am
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/Makefile.am 2007-11-17 01:30:59 UTC (rev 8300)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/Makefile.am 2007-11-18 23:59:14 UTC (rev 8301)
@@ -16,8 +16,10 @@
subpkginclude_HEADERS = \
Exception.hh \
Exception.icc \
+ GocadVoxet.hh \
GravityField.hh \
GravityField.icc \
+ SCECCVMH.hh \
SpatialDB.hh \
SpatialDB.icc \
SimpleDB.hh \
Added: cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc 2007-11-17 01:30:59 UTC (rev 8300)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc 2007-11-18 23:59:14 UTC (rev 8301)
@@ -0,0 +1,213 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+//#include <portinfo>
+
+#include "SCECCVMH.hh" // Implementation of class methods
+
+#include "GocadVoxet.hh" // USES GocadVoxet
+
+#include "spatialdata/geocoords/CSGeoProj.hh" // USES CSGeoProj
+#include "spatialdata/geocoords/Projector.hh" // USES Projector
+#include "spatialdata/geocoords/Converter.hh" // USES Converter
+
+#include <math.h> // USES pow()
+#include <sstream> // USES std::ostringstream
+#include <stdexcept> // USES std::logic_error
+#include <assert.h> // USES assert()
+
+#include <iostream> // TEMPORARY
+
+// ----------------------------------------------------------------------
+// Constructor
+spatialdata::spatialdb::SCECCVMH::SCECCVMH(void) :
+ SpatialDB("SCEC CVM-H"),
+ _dataDir("."),
+ _laLowResVp(0),
+ _laLowResTag(0),
+ _laHighResVp(0),
+ _laHighResTag(0),
+ _crustMantleVp(0),
+ _crustMantleVs(0),
+ _crustMantleTag(0),
+ _topoElev(0),
+ _baseDepth(0),
+ _mohoDepth(0),
+ _csUTM(new geocoords::CSGeoProj)
+{ // constructor
+ geocoords::Projector proj;
+ proj.projection("utm");
+ proj.units("m");
+ proj.projOptions("+zone=11");
+
+ _csUTM->datumHoriz("WGS84");
+ _csUTM->datumVert("mean sea level");
+ _csUTM->projector(proj);
+ _csUTM->initialize();
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+spatialdata::spatialdb::SCECCVMH::~SCECCVMH(void)
+{ // destructor
+ delete _laLowResVp; _laLowResVp = 0;
+ delete _laLowResTag; _laLowResTag = 0;
+ delete _laHighResVp; _laHighResVp = 0;
+ delete _laHighResTag; _laHighResTag = 0;
+ delete _crustMantleVp; _crustMantleVp = 0;
+ delete _crustMantleVs; _crustMantleVs = 0;
+ delete _crustMantleTag; _crustMantleTag = 0;
+ delete _topoElev; _topoElev = 0;
+ delete _baseDepth; _baseDepth = 0;
+ delete _mohoDepth; _mohoDepth = 0;
+ delete _csUTM; _csUTM = 0;
+} // destructor
+
+// ----------------------------------------------------------------------
+// Set directory containing SCEC CVM-H data files.
+void
+spatialdata::spatialdb::SCECCVMH::dataDir(const char* dir)
+{ // dataDir
+ _dataDir = dir;
+} // dataDir
+
+// ----------------------------------------------------------------------
+// Open the database and prepare for querying.
+void
+spatialdata::spatialdb::SCECCVMH::open(void)
+{ // open
+ if (0 == _laLowResVp)
+ _laLowResVp = new GocadVoxet;
+ _laLowResVp->read(_dataDir.c_str(), "LA_LR.vo", "\"VINT1D\"");
+
+ if (0 == _laLowResTag)
+ _laLowResTag = new GocadVoxet;
+ _laLowResTag->read(_dataDir.c_str(), "LA_LR.vo", "\"flag\"");
+
+ if (0 == _laHighResVp)
+ _laHighResVp = new GocadVoxet;
+ _laHighResVp->read(_dataDir.c_str(), "LA_HR.vo", "\"vp\"");
+ if (0 == _laHighResTag)
+ _laHighResTag = new GocadVoxet;
+ _laHighResTag->read(_dataDir.c_str(), "LA_HR.vo", "\"tag\"");
+
+ if (0 == _crustMantleVp)
+ _crustMantleVp = new GocadVoxet;
+ _crustMantleVp->read(_dataDir.c_str(), "CM.vo", "\"cvp\"");
+ if (0 == _crustMantleVs)
+ _crustMantleVs = new GocadVoxet;
+ _crustMantleVs->read(_dataDir.c_str(), "CM.vo", "\"cvs\"");
+ if (0 == _crustMantleTag)
+ _crustMantleTag = new GocadVoxet;
+ _crustMantleTag->read(_dataDir.c_str(), "CM.vo", "\"tag\"");
+
+ if (0 == _topoElev)
+ _topoElev = new GocadVoxet;
+ _topoElev->read(_dataDir.c_str(), "topo.vo", "\"topo\"");
+
+ if (0 == _baseDepth)
+ _baseDepth = new GocadVoxet;
+ _baseDepth->read(_dataDir.c_str(), "base.vo", "\"base\"");
+
+ if (0 == _mohoDepth)
+ _mohoDepth = new GocadVoxet;
+ _mohoDepth->read(_dataDir.c_str(), "moho.vo", "\"moho\"");
+} // open
+
+// ----------------------------------------------------------------------
+// Close the database.
+void
+spatialdata::spatialdb::SCECCVMH::close(void)
+{ // close
+ delete _laLowResVp; _laLowResVp = 0;
+ delete _laLowResTag; _laLowResTag = 0;
+ delete _laHighResVp; _laHighResVp = 0;
+ delete _laHighResTag; _laHighResTag = 0;
+ delete _crustMantleVp; _crustMantleVp = 0;
+ delete _crustMantleVs; _crustMantleVs = 0;
+ delete _crustMantleTag; _crustMantleTag = 0;
+ delete _topoElev; _topoElev = 0;
+ delete _baseDepth; _baseDepth = 0;
+ delete _mohoDepth; _mohoDepth = 0;
+} // close
+
+// ----------------------------------------------------------------------
+// Set values to be returned by queries.
+void
+spatialdata::spatialdb::SCECCVMH::queryVals(const char** names,
+ const int numVals)
+{ // queryVals
+} // queryVals
+
+// ----------------------------------------------------------------------
+// Query the database.
+int
+spatialdata::spatialdb::SCECCVMH::query(double* vals,
+ const int numVals,
+ const double* coords,
+ const int numDims,
+ const spatialdata::geocoords::CoordSys* csQuery)
+{ // query
+ // Convert coordinates to UTM
+ for (int i=0; i < numDims; ++i)
+ _xyzUTM[i] = coords[i];
+ spatialdata::geocoords::Converter::convert(_xyzUTM, 1, numDims,
+ _csUTM, csQuery);
+
+ assert(0 != _topoElev);
+ double elev = 0.0;
+ double baseDepth = 0.0;
+ double mohoDepth = 0.0;
+ double vp = 0.0;
+ double vs = 0.0;
+ double density = 0.0;
+ int outsideVoxet =_topoElev->query(&elev, _xyzUTM);
+ outsideVoxet = _baseDepth->query(&baseDepth, _xyzUTM);
+ outsideVoxet = _mohoDepth->query(&mohoDepth, _xyzUTM);
+ outsideVoxet = _laLowResVp->query(&vp, _xyzUTM);
+ if (!outsideVoxet) {
+ std::cout << "Using LR";
+ double vpHR = 0.0;
+ outsideVoxet = _laHighResVp->query(&vpHR, _xyzUTM);
+ if (!outsideVoxet) {
+ std::cout << "Using HR";
+ vp = vpHR;
+ } // if
+ if (vp != 1480.0) {
+ if (vp < 4250.0)
+ vs = (vp > 1500.0) ? (vp-1360.0/1.16) : (1500.0-1360.0)/1.16;
+ else
+ vs = 785.8 - 1.2344*vp + 794.9 * pow(vp/1000.0,2)
+ -123.8 * pow(vp/1000.0,3) + 6.4 * pow(vp/1000.0,4);
+ if (vp > 2160.0)
+ density = vp / 3.0 + 1280.0;
+ else
+ density = 2000.0;
+ } else
+ density = 1000.0;
+ } else {
+ outsideVoxet = _crustMantleVp->query(&vp, _xyzUTM);
+ if (!outsideVoxet)
+ outsideVoxet = _crustMantleVs->query(&vs, _xyzUTM);
+ density = vp / 3.0 + 1280.0;
+ } // if/else
+
+ std::cout << "elev: " << elev << std::endl;
+ std::cout << "basement depth: " << baseDepth << std::endl;
+ std::cout << "moho depth: " << mohoDepth << std::endl;
+ std::cout << "vp: " << vp << std::endl;
+ std::cout << "vs: " << vs << std::endl;
+ std::cout << "density: " << density << std::endl;
+} // query
+
+
+// End of file
Added: cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.hh
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.hh 2007-11-17 01:30:59 UTC (rev 8300)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.hh 2007-11-18 23:59:14 UTC (rev 8301)
@@ -0,0 +1,118 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(spatialdata_spatialdb_sceccvmh_hh)
+#define spatialdata_spatialdb_sceccvmh_hh
+
+#include "spatialdata/spatialdb/SpatialDB.hh" // ISA SpatialDB
+
+#include <string>> // HASA std::string
+
+namespace spatialdata {
+ namespace spatialdb {
+ class SCECCVMH;
+
+ class GocadVoxet; // HASA GocadVoxet
+ class TestSCECCVMH; // unit testing
+ } // spatialdb
+
+ namespace geocoords {
+ class CSGeoProj; // HASA CSGeoProj
+ } // geocoords
+} // spatialdata
+
+class spatialdata::spatialdb::SCECCVMH : SpatialDB
+{ // SCECCVMH
+ friend class TestSCECCVMH; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ SCECCVMH(void);
+
+ /// Destructor
+ ~SCECCVMH(void);
+
+ /** Set directory containing SCEC CVM-H data files.
+ *
+ * @param dir Directory containing data files.
+ */
+ void dataDir(const char* dir);
+
+ /// Open the database and prepare for querying.
+ void open(void);
+
+ /// Close the database.
+ void close(void);
+
+ /** Set values to be returned by queries.
+ *
+ * @pre Must call open() before queryVals()
+ *
+ * @param names Names of values to be returned in queries
+ * @param numVals Number of values to be returned in queries
+ */
+ void queryVals(const char** names,
+ const int numVals);
+
+ /** Query the database.
+ *
+ * @note pVals should be preallocated to accommodate numVals values.
+ *
+ * @pre Must call open() before query()
+ *
+ * @param vals Array for computed values (output from query), must be
+ * allocated BEFORE calling query().
+ * @param numVals Number of values expected (size of pVals array)
+ * @param coords Coordinates of point for query
+ * @param numDims Number of dimensions for coordinates
+ * @param pCSQuery Coordinate system of coordinates
+ *
+ * @returns 0 on success, 1 on failure (i.e., could not interpolate)
+ */
+ int query(double* vals,
+ const int numVals,
+ const double* coords,
+ const int numDims,
+ const spatialdata::geocoords::CoordSys* pCSQuery);
+
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ SCECCVMH(const SCECCVMH&); ///< Not implemented
+ const SCECCVMH& operator=(const SCECCVMH&); ///< Not implemented
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ double _xyzUTM[3];
+ std::string _dataDir;
+ GocadVoxet* _laLowResVp;
+ GocadVoxet* _laLowResTag;
+ GocadVoxet* _laHighResVp;
+ GocadVoxet* _laHighResTag;
+ GocadVoxet* _crustMantleVp;
+ GocadVoxet* _crustMantleVs;
+ GocadVoxet* _crustMantleTag;
+ GocadVoxet* _topoElev;
+ GocadVoxet* _baseDepth;
+ GocadVoxet* _mohoDepth;
+ geocoords::CSGeoProj* _csUTM;
+
+}; // SCECCVMH
+
+#endif // spatialdata_spatialdb_sceccvmh_hh
+
+
+// End of file
More information about the cig-commits
mailing list