[cig-commits] r8983 - in cs/spatialdata-0.1/trunk: libsrc/spatialdb
modulesrc/spatialdb spatialdata/spatialdb
tests/libtests/spatialdb tests/pytests/spatialdb
brad at geodynamics.org
brad at geodynamics.org
Thu Jan 3 14:45:00 PST 2008
Author: brad
Date: 2008-01-03 14:44:59 -0800 (Thu, 03 Jan 2008)
New Revision: 8983
Added:
cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.icc
cs/spatialdata-0.1/trunk/tests/pytests/spatialdb/TestSCECCVMH.py
Modified:
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/Makefile.am
cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc
cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.hh
cs/spatialdata-0.1/trunk/modulesrc/spatialdb/spatialdb.pyxe.src
cs/spatialdata-0.1/trunk/spatialdata/spatialdb/SCECCVMH.py
cs/spatialdata-0.1/trunk/tests/libtests/spatialdb/TestSCECCVMH.cc
cs/spatialdata-0.1/trunk/tests/libtests/spatialdb/TestSCECCVMH.hh
cs/spatialdata-0.1/trunk/tests/pytests/spatialdb/TestSCECCVMH.py.in
Log:
Added queryNearest() for SCEC CVM-H and squashing flag (and implementation of squashing).
Modified: cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc 2008-01-03 22:44:59 UTC (rev 8983)
@@ -92,6 +92,49 @@
} // query
// ----------------------------------------------------------------------
+// Query voxet for value at nearest location.
+int
+spatialdata::spatialdb::GocadVoxet::queryNearest(double* value,
+ const double pt[3]) const
+{ // queryNearest
+ // Compute indices of voxet containing pt
+ const int numX = _geometry.n[0];
+ const int numY = _geometry.n[1];
+ const int numZ = _geometry.n[2];
+ int indexX =
+ round( (pt[0] - (_geometry.o[0]+_geometry.min[0]))*_geometry.scale[0]);
+ int indexY =
+ round( (pt[1] - (_geometry.o[1]+_geometry.min[1]))*_geometry.scale[1]);
+ int indexZ =
+ round( (pt[2] - (_geometry.o[2]+_geometry.min[2]))*_geometry.scale[2]);
+
+ // Correct to range of voxet, if necessary.
+ int flag = 0;
+ if (indexX < 0)
+ indexX = 0;
+ else if (indexX >= numX)
+ indexX = numX - 1;
+ if (indexY < 0)
+ indexY = 0;
+ else if (indexY >= numY)
+ indexY = numY - 1;
+ if (indexZ < 0)
+ indexZ = 0;
+ else if (indexZ >= numZ)
+ indexZ = numZ - 1;
+
+ // Make sure voxet is in range
+ assert(indexX >= 0 && indexX < numX);
+ assert(indexY >= 0 && indexY < numY);
+ assert(indexZ >= 0 && indexZ < numZ);
+
+ int indexV = indexZ*numY*numX + indexY*numX + indexX;
+ *value = _data[indexV];
+
+ return flag;
+} // queryNearest
+
+// ----------------------------------------------------------------------
// Read voxet file.
void
spatialdata::spatialdb::GocadVoxet::_readVoxetFile(const char* filename,
Modified: cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.hh
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.hh 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.hh 2008-01-03 22:44:59 UTC (rev 8983)
@@ -59,6 +59,15 @@
int query(double* value,
const double pt[3]) const;
+ /** Query voxet for value at nearest location.
+ *
+ * @param value Value for result.
+ * @param pt Location of query.
+ * @returns 0 if pt is inside voxet, 1 if outside voxet.
+ */
+ int queryNearest(double* value,
+ const double pt[3]) const;
+
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: cs/spatialdata-0.1/trunk/libsrc/spatialdb/Makefile.am
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/Makefile.am 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/Makefile.am 2008-01-03 22:44:59 UTC (rev 8983)
@@ -20,6 +20,7 @@
GravityField.hh \
GravityField.icc \
SCECCVMH.hh \
+ SCECCVMH.icc \
SpatialDB.hh \
SpatialDB.icc \
SimpleDB.hh \
Modified: cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc 2008-01-03 22:44:59 UTC (rev 8983)
@@ -41,8 +41,10 @@
_baseDepth(0),
_mohoDepth(0),
_csUTM(new geocoords::CSGeoProj),
+ _squashLimit(-2000.0),
_queryVals(0),
- _querySize(0)
+ _querySize(0),
+ _squashTopo(false)
{ // constructor
geocoords::Projector proj;
proj.projection("utm");
@@ -72,19 +74,13 @@
delete _mohoDepth; _mohoDepth = 0;
delete _csUTM; _csUTM = 0;
+ _squashLimit = 0.0;
delete[] _queryVals; _queryVals = 0;
_querySize = 0;
+ _squashTopo = 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)
@@ -187,6 +183,7 @@
} // for
} // queryVals
+#include <iostream>
// ----------------------------------------------------------------------
// Query the database.
int
@@ -218,6 +215,15 @@
spatialdata::geocoords::Converter::convert(_xyzUTM, 1, numDims,
_csUTM, csQuery);
+ bool haveTopo = false;
+ double topoElev = 0;
+ if (_squashTopo && _xyzUTM[2] > _squashLimit) {
+ double z = 0;
+ _topoElev->queryNearest(&topoElev, _xyzUTM);
+ haveTopo = true;
+ _xyzUTM[2] += topoElev;
+ } // if
+
int outsideVoxet = 0;
int queryFlag = 0;
bool haveVp = false;
@@ -256,20 +262,23 @@
vals[iVal] = _calcVs(vp);
break;
case QUERY_TOPOELEV :
- assert(0 != _topoElev);
- outsideVoxet = _topoElev->query(&vals[iVal], _xyzUTM);
- if (outsideVoxet)
- queryFlag = outsideVoxet;
+ if (!haveTopo) {
+ assert(0 != _topoElev);
+ outsideVoxet = _topoElev->queryNearest(&vals[iVal], _xyzUTM);
+ if (outsideVoxet)
+ queryFlag = outsideVoxet;
+ } else
+ vals[iVal] = topoElev;
break;
case QUERY_BASEDEPTH :
assert(0 != _baseDepth);
- outsideVoxet = _baseDepth->query(&vals[iVal], _xyzUTM);
+ outsideVoxet = _baseDepth->queryNearest(&vals[iVal], _xyzUTM);
if (outsideVoxet)
queryFlag = outsideVoxet;
break;
case QUERY_MOHODEPTH :
assert(0 != _mohoDepth);
- outsideVoxet = _mohoDepth->query(&vals[iVal], _xyzUTM);
+ outsideVoxet = _mohoDepth->queryNearest(&vals[iVal], _xyzUTM);
if (outsideVoxet)
queryFlag = outsideVoxet;
break;
@@ -294,18 +303,9 @@
if (!outsideVoxet)
*vp = vpHR;
else
- outsideVoxet = 0; // reset outsideVoxet flag to low-res value
- } else {
- outsideVoxet = _crustMantleVp->query(vp, _xyzUTM);
- if (outsideVoxet) {
- const double vpBg = _backgroundVp();
- if (vpBg > 0.0) {
- outsideVoxet = 0;
- *vp = vpBg;
- } else
- outsideVoxet = 1;
- } // if
- } // else
+ outsideVoxet = 0; // use low-res value
+ } else
+ outsideVoxet = _crustMantleVp->queryNearest(vp, _xyzUTM);
return outsideVoxet;
} // _queryVp
@@ -324,9 +324,9 @@
if (!outsideVoxet)
*tag = tagHR;
else
- outsideVoxet = 0; // reset outsideVoxet flag to low-res value
+ outsideVoxet = 0; // use low-res value
} else
- outsideVoxet = _crustMantleTag->query(tag, _xyzUTM);
+ outsideVoxet = _crustMantleTag->queryNearest(tag, _xyzUTM);
return outsideVoxet;
} // _queryTag
@@ -369,33 +369,5 @@
return vs;
} // _calcVs
-// ----------------------------------------------------------------------
-// Compute vp for background model.
-double
-spatialdata::spatialdb::SCECCVMH::_backgroundVp(void)
-{ // _backgroundVp
- const double z = _xyzUTM[2];
- double vp = -99999.0;
- assert(0 != _topoElev);
- double elev = 0.0;
- const int outsideVoxet = _topoElev->query(&elev, _xyzUTM);
- if (outsideVoxet || z < -10000.0) {
- // if outside horizontal extent of domain or deeper; assume domain
- // is at leat 10km deep (z < -10000.0)
-
- if (z < -35000.0)
- vp = 7800.0;
- else if (z < -15000.0)
- vp = 7800.0 - (7800.0 - 7000.0) / 20000.0 * (z + 35000.0);
- else if (z < 0.0)
- vp = 7000.0-8.8888889e-06*(z+15000.0)*(z+15000.0);
- else
- vp = 5000.0;
- } // if
-
- return vp;
-} // _backgroundVp
-
-
// End of file
Modified: cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.hh
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.hh 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.hh 2008-01-03 22:44:59 UTC (rev 8983)
@@ -49,6 +49,15 @@
*/
void dataDir(const char* dir);
+ /** Set squashed topography/bathymetry flag and minimum elevation of
+ * squashing.
+ *
+ * @param flag True if squashing, false otherwise.
+ * @param limit Minimum elevation for squashing.
+ */
+ void squash(const bool flag,
+ const double limit =-2000.0);
+
/// Open the database and prepare for querying.
void open(void);
@@ -139,17 +148,6 @@
static
double _calcVs(const double vp);
- /** Compute vp for background model.
- *
- * Distribution of Vp with depth that provides approximate match to
- * distrbution of Vp on edges of SCEC CVM-H. This is similar to the
- * starting model in Hauksson's tomographic model (Hauksson, JGR,
- * 2000).
- *
- * @returns Vp in m/s.
- */
- double _backgroundVp(void);
-
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
@@ -167,11 +165,15 @@
GocadVoxet* _mohoDepth;
geocoords::CSGeoProj* _csUTM;
+ double _squashLimit; ///< Elevation above which topography is squashed.
int* _queryVals; ///< Indices of values to be returned in queries.
int _querySize; ///< Number of values requested to be returned in queries.
+ bool _squashTopo; ///< Squash topography/bathymetry to sea level.
}; // SCECCVMH
+#include "SCECCVMH.icc" // inline methods
+
#endif // spatialdata_spatialdb_sceccvmh_hh
Added: cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.icc
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.icc 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.icc 2008-01-03 22:44:59 UTC (rev 8983)
@@ -0,0 +1,34 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// <LicenseText>
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(spatialdata_spatialdb_sceccvmh_hh)
+#error "SCECCVMH.icc must only be included from SCECCVMH.hh"
+#endif
+
+// Set directory containing SCEC CVM-H data files.
+inline
+void
+spatialdata::spatialdb::SCECCVMH::dataDir(const char* dir) {
+ _dataDir = dir;
+}
+
+// Set squashed topography/bathymetry flag and minimum elevation of
+inline
+void
+spatialdata::spatialdb::SCECCVMH::squash(const bool flag,
+ const double limit) {
+ _squashTopo = flag;
+ _squashLimit = limit;
+}
+
+
+// End of file
Modified: cs/spatialdata-0.1/trunk/modulesrc/spatialdb/spatialdb.pyxe.src
===================================================================
--- cs/spatialdata-0.1/trunk/modulesrc/spatialdb/spatialdb.pyxe.src 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/modulesrc/spatialdb/spatialdb.pyxe.src 2008-01-03 22:44:59 UTC (rev 8983)
@@ -619,4 +619,19 @@
return
+ def squash(self, flag, limit):
+ """
+ Set squashed topography/bathymetry flag and minimum elevation of
+ squashing.
+ """
+ # create shim for method 'squash'
+ #embed{ void SCECCVMH_squash(void* pObj, int flag, double limit)
+ assert(0 != pObj);
+ ((spatialdata::spatialdb::SCECCVMH*) pObj)->squash(flag, limit);
+ #}embed
+
+ SCECCVMH_squash(self.thisptr, flag, limit)
+ return
+
+
# End of file
Modified: cs/spatialdata-0.1/trunk/spatialdata/spatialdb/SCECCVMH.py
===================================================================
--- cs/spatialdata-0.1/trunk/spatialdata/spatialdb/SCECCVMH.py 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/spatialdata/spatialdb/SCECCVMH.py 2008-01-03 22:44:59 UTC (rev 8983)
@@ -39,6 +39,8 @@
##
## \b Properties
## @li \b dataDir Directory containing SCEC CVM-H data files.
+ ## @li \b squash Squash topography/bathymetry to sea level.
+ ## @li \b squash_limit Elevation above which topography is squashed.
##
## \b Facilities
## @li none
@@ -48,7 +50,15 @@
dataDir = pyre.inventory.str("data_dir", default=".")
dataDir.meta['tip'] = "Directory containing SCEC CVM-H data files."
+ squash = pyre.inventory.bool("squash", default=False)
+ squash.meta['tip'] = "Squash topography/bathymetry to sea level."
+ from pyre.units.length import km
+ squashLimit = pyre.inventory.dimensional("squash_limit",
+ default=-2.0*km)
+ squashLimit.meta['tip'] = "Elevation above which topography is squashed."
+
+
# PUBLIC METHODS /////////////////////////////////////////////////////
def __init__(self, name="sceccvmh"):
@@ -68,6 +78,7 @@
"""
SpatialDB.initialize(self)
self.cppHandle.dataDir(self.dataDir)
+ self.cppHandle.squash(self.squash, self.squashLimit.value)
return
@@ -79,6 +90,8 @@
"""
SpatialDB._configure(self)
self.dataDir = self.inventory.dataDir
+ self.squash = self.inventory.squash
+ self.squashLimit = self.inventory.squashLimit
return
Modified: cs/spatialdata-0.1/trunk/tests/libtests/spatialdb/TestSCECCVMH.cc
===================================================================
--- cs/spatialdata-0.1/trunk/tests/libtests/spatialdb/TestSCECCVMH.cc 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/tests/libtests/spatialdb/TestSCECCVMH.cc 2008-01-03 22:44:59 UTC (rev 8983)
@@ -57,6 +57,28 @@
} // testDataDir
// ----------------------------------------------------------------------
+// Test squashed().
+void
+spatialdata::spatialdb::TestSCECCVMH::testSquashed(void)
+{ // testSquashed
+ SCECCVMH db;
+
+ const double limitDefault = -2000.0;
+ const double limit = -4000.0;
+
+ CPPUNIT_ASSERT_EQUAL(false, db._squashTopo);
+ CPPUNIT_ASSERT_EQUAL(limitDefault, db._squashLimit);
+
+ db.squash(true);
+ CPPUNIT_ASSERT_EQUAL(true, db._squashTopo);
+ CPPUNIT_ASSERT_EQUAL(limitDefault, db._squashLimit);
+
+ db.squash(true, limit);
+ CPPUNIT_ASSERT_EQUAL(true, db._squashTopo);
+ CPPUNIT_ASSERT_EQUAL(limit, db._squashLimit);
+} // testDataDir
+
+// ----------------------------------------------------------------------
// Test queryVals().
void
spatialdata::spatialdb::TestSCECCVMH::testQueryVals(void)
@@ -125,7 +147,7 @@
cs.datumVert("mean sea level");
cs.initialize();
- const int numLocs = 15;
+ const int numLocs = 17;
const int spaceDim = 3;
const double lonlatelev[] = {
-118.560000, 32.550000, -2450.00,
@@ -143,6 +165,9 @@
-125.000000, 35.000000, -4000.0,
-125.000000, 35.000000, 0.0,
+ -115.000000, 30.000000, -800000.0, // outside domain
+ -115.000000, 30.000000, -4539.0,
+
-117.989344, 34.034148, 5000.00, // above domain
-117.989344, 34.034148, -500000.0, // below domain
};
@@ -170,20 +195,24 @@
-27857.322266, 106.671280, 0.0, -2361.386963, 6351.574219, 3397.191406, 3707.466394,
-28512.111328, 93.515053, 2.0, -2860.919189, 4932.721191, 2924.240397, 2968.477965,
- -99999.0, -99999.0, -99999.0, -99999.0, 7800.0, 3880.0, 4459.294240,
- -99999.0, -99999.0, -99999.0, -99999.0, 7800.0, 3880.0, 4459.294240,
- -99999.0, -99999.0, -99999.0, -99999.0, 7532.0, 3790.666667, 4282.106907,
- -99999.0, -99999.0, -99999.0, -99999.0, 7000.0, 3613.333333, 3998.100000,
- -99999.0, -99999.0, -99999.0, -99999.0, 5924.4444444, 3254.814814, 3514.068087,
- -99999.0, -99999.0, -99999.0, -99999.0, 5000.0, 2946.666667, 3011.300000,
+ -26312.804688, -50.465275, 1.0, -2875.716553, 7819.771973, 3886.590658, 4473.594442,
+ -26312.804688, -50.465275, 1.0, -2875.716553, 7800.935059, 3880.311686, 4459.966373,
+ -26312.804688, -50.465275, 0.0, -2875.716553, 6925.244140, 3588.414714, 3962.891997,
+ -26312.804688, -50.465275, 0.0, -2875.716553, 6585.769531, 3475.256510, 3810.172109,
+ -26312.804688, -50.465275, 0.0, -2875.716553, 6585.769531, 3475.256510, 3810.172109,
+ -26312.804688, -50.465275, 0.0, -2875.716553, 6585.769531, 3475.256510, 3810.172109,
- -28512.111328, 93.515053, -99999.0, -2860.919189, -99999.0, -99999.0, -99999.0,
- -28512.111328, 93.515053, -99999.0, -2860.919189, 7800.0, 3880.0, 4459.294240,
+ -26906.470703, 3.366069, 1.0, -973.515442, 7800.0, 3880.0, 4459.29424,
+ -26906.470703, 3.366069, 2.0, -973.515442, 6642.105469, 3494.035156, 3834.988779,
+
+ -28512.111328, 93.515053, 0.0, -2860.919189, 6649.152832, 3496.384277, 3838.102107,
+ -28512.111328, 93.515053, 0.0, -2860.919189, 8159.479492, 3999.826497, 4751.5155415,
};
const int flags[] = {
0, 0, 0, 0, 0, 0, 0,
- 1, 1, 1, 1, 1, 1,
- 1, 1
+ 0, 0, 0, 0, 0, 0,
+ 0, 0,
+ 0, 0
};
db.queryVals(queryNames, querySize);
@@ -218,19 +247,22 @@
2673.790690, 3.0,
3397.191406, 0.0,
2924.240397, 2.0,
- 3880.000000, -99999.00,
- 3880.000000, -99999.00,
- 3790.666667, -99999.00,
- 3613.333333, -99999.00,
- 3254.814814, -99999.00,
- 2946.666667, -99999.00,
- -99999.0, -99999.00,
- 3880.000000, -99999.00,
+ 3886.590657, 1.0,
+ 3880.311686, 1.0,
+ 3588.414714, 0.0,
+ 3475.256510, 0.0,
+ 3475.256510, 0.0,
+ 3475.256510, 0.0,
+ 3880.000000, 1.0,
+ 3494.035156, 2.0,
+ 3496.384277, 0.0,
+ 3999.826497, 0.0,
};
const int flags[] = {
0, 0, 0, 0, 0, 0, 0,
- 1, 1, 1, 1, 1, 1,
- 1, 1
+ 0, 0, 0, 0, 0, 0,
+ 0, 0,
+ 0, 0
};
db.queryVals(queryNames, querySize);
@@ -264,19 +296,25 @@
2432.217302,
3707.466394,
2968.477965,
- 4459.29424,
- 4459.29424,
- 4282.106907,
- 3998.100000,
- 3514.068087,
- 3011.300000,
- -99999.0,
- 4459.29424,
+
+ 4473.594442,
+ 4459.966373,
+ 3962.891997,
+ 3810.172109,
+ 3810.172109,
+ 3810.172109,
+
+ 4459.29424,
+ 3834.988779,
+
+ 3838.102107,
+ 4751.5155415,
};
const int flags[] = {
0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0,
- 1, 0
+ 0, 0,
+ 0, 0
};
db.queryVals(queryNames, querySize);
@@ -299,6 +337,110 @@
db.close();
} // testQuery
+
+// ----------------------------------------------------------------------
+// Test query() with squashing.
+void
+spatialdata::spatialdb::TestSCECCVMH::testQuerySquashed(void)
+{ // testQuerySquashed
+ SCECCVMH db;
+ db.dataDir(SCECCVMH_DATADIR);
+
+ spatialdata::geocoords::CSGeo cs;
+ cs.ellipsoid("clrk66");
+ cs.datumHoriz("NAD27");
+ cs.datumVert("mean sea level");
+ cs.initialize();
+
+ const int numLocs = 17;
+ const int spaceDim = 3;
+ const double lonlatelev[] = {
+ -118.560000, 32.550000, -2450.00,
+ -118.513208, 33.884888, -1400.00+56.893856,
+ -118.520000, 34.120000, -1400.00-489.975189,
+ -116.400000, 32.340000, -1000.00-801.209961,
+ -118.337765, 34.095691, -1770.00-106.671280,
+ -118.337765, 34.095691, -17700.00,
+ -117.989344, 34.034148, -3000.00,
+
+ -125.000000, 35.000000, -40000.0, // outside domain
+ -125.000000, 35.000000, -35000.0,
+ -125.000000, 35.000000, -28300.0,
+ -125.000000, 35.000000, -15000.0,
+ -125.000000, 35.000000, -4000.0,
+ -125.000000, 35.000000, 0.0-50.465275,
+
+ -115.000000, 30.000000, -800000.0, // outside domain
+ -115.000000, 30.000000, -4539.0,
+
+ -117.989344, 34.034148, 5000.00, // above domain
+ -117.989344, 34.034148, -500000.0, // below domain
+ };
+ const double tolerance = 1.0e-06;
+
+ db.open();
+ { // all values
+ const int querySize = 7;
+ const char* queryNames[] = {
+ "moho-depth",
+ "topo-elev",
+ "vp-tag",
+ "Basement-Depth",
+ "vP",
+ "DenSity",
+ "VS"
+ };
+ const double values[] = {
+ -27991.001953, -1115.499268, 2.0, -1326.415405, 5560.209473, 3133.403158, 3333.334131,
+ //-27991.00, -1115.50, 2.00, -1326.42, 5540.545410, 3126.848633, 3323.004705,
+ -27165.316406, -56.893856, 2.0, -1431.710449, 4384.126953, 2741.375651, 2584.759474,
+ -31178.105469, 489.975189, 2.0, 26.471289, 4142.553223, 2660.851074, 2398.752778,
+ -34526.414062, 801.209961, 2.0, 631.073059, 5226.578125, 3022.192708, 3148.777736,
+ -27857.322266, 106.671280, 3.0, -2361.386963, 4181.372070, 2673.790690, 2432.217302,
+ -27857.322266, 106.671280, 0.0, -2361.386963, 6351.574219, 3397.191406, 3707.466394,
+ -28512.111328, 93.515053, 2.0, -2860.919189, 4932.721191, 2924.240397, 2968.477965,
+
+ -26312.804688, -50.465275, 1.0, -2875.716553, 7819.771973, 3886.590658, 4473.594442,
+ -26312.804688, -50.465275, 1.0, -2875.716553, 7800.935059, 3880.311686, 4459.966373,
+ -26312.804688, -50.465275, 0.0, -2875.716553, 6925.244140, 3588.414714, 3962.891997,
+ -26312.804688, -50.465275, 0.0, -2875.716553, 6585.769531, 3475.256510, 3810.172109,
+ -26312.804688, -50.465275, 0.0, -2875.716553, 6585.769531, 3475.256510, 3810.172109,
+ -26312.804688, -50.465275, 0.0, -2875.716553, 6585.769531, 3475.256510, 3810.172109,
+
+ -26906.470703, 3.36606908, 1.0, -973.515442, 7800.0, 3880.0, 4459.29424,
+ -26906.470703, 3.36606908, 2.0, -973.515442, 6642.105469, 3494.035156, 3834.988779,
+
+ -28512.111328, 93.515053, 0.0, -2860.919189, 6649.152832, 3496.384277, 3838.102107,
+ -28512.111328, 93.515053, 0.0, -2860.919189, 8159.479492, 3999.826497, 4751.5155415,
+ };
+ const int flags[] = {
+ 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0,
+ 0, 0
+ };
+
+ db.squash(true);
+ db.queryVals(queryNames, querySize);
+
+ double data[querySize];
+
+ for (int iLoc=0; iLoc < numLocs; ++iLoc) {
+ int err = db.query(data, querySize, &lonlatelev[iLoc*spaceDim], spaceDim, &cs);
+ CPPUNIT_ASSERT_EQUAL(flags[iLoc], err);
+
+ for (int iVal=0; iVal < querySize; ++iVal) {
+ const double dataE = values[iLoc*querySize+iVal];
+ if (fabs(dataE) > tolerance)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, data[iVal]/dataE, tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(dataE, data[iVal], tolerance);
+ } // for
+ } // for
+ } // all values
+
+ db.close();
+} // testQuerySquashed
#endif
// ----------------------------------------------------------------------
Modified: cs/spatialdata-0.1/trunk/tests/libtests/spatialdb/TestSCECCVMH.hh
===================================================================
--- cs/spatialdata-0.1/trunk/tests/libtests/spatialdb/TestSCECCVMH.hh 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/tests/libtests/spatialdb/TestSCECCVMH.hh 2008-01-03 22:44:59 UTC (rev 8983)
@@ -40,11 +40,13 @@
CPPUNIT_TEST( testConstructor );
CPPUNIT_TEST( testLabel );
CPPUNIT_TEST( testDataDir );
+ CPPUNIT_TEST( testSquashed );
CPPUNIT_TEST( testQueryVals );
CPPUNIT_TEST( testCalcDensity );
CPPUNIT_TEST( testCalcVs );
#if defined(SCECCVMH_DATADIR)
CPPUNIT_TEST( testQuery );
+ CPPUNIT_TEST( testQuerySquashed );
#endif
CPPUNIT_TEST_SUITE_END();
@@ -61,12 +63,18 @@
/// Test dataDir()
void testDataDir(void);
+ /// Test squashed()
+ void testSquashed(void);
+
/// Test queryVals()
void testQueryVals(void);
/// Test query()
void testQuery(void);
+ /// Test querySquashed()
+ void testQuerySquashed(void);
+
/// Test calcDensity()
void testCalcDensity(void);
Added: cs/spatialdata-0.1/trunk/tests/pytests/spatialdb/TestSCECCVMH.py
===================================================================
--- cs/spatialdata-0.1/trunk/tests/pytests/spatialdb/TestSCECCVMH.py 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/tests/pytests/spatialdb/TestSCECCVMH.py 2008-01-03 22:44:59 UTC (rev 8983)
@@ -0,0 +1,110 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+import unittest
+
+import numpy
+
+class TestSCECCVMH(unittest.TestCase):
+
+ def setUp(self):
+ from spatialdata.spatialdb.SCECCVMH import SCECCVMH
+ db = SCECCVMH()
+ db._configure()
+ db.dataDir = "/Users/brad/data/sceccvm-h/vx52/bin"
+ db.initialize()
+ self._db = db
+ return
+
+
+ def test_query(self):
+ locs = numpy.array( [[-118.520000, 34.120000, -1400.00],
+ [-116.400000, 32.340000, -1000.00]],
+ numpy.float64)
+ from spatialdata.geocoords.CSGeo import CSGeo
+ cs = CSGeo()
+ cs._configure()
+ cs.datumHoriz = "NAD27"
+ cs.datumVert = "mean sea level"
+ cs.ellipsoid = "clrk66"
+ cs.initialize()
+ queryVals = ["topo-elev", "moho-depth", "density"]
+ dataE = numpy.array([[489.975189, -31178.105469, 2660.851074],
+ [801.209961, -34526.414062, 3022.192708]],
+ numpy.float64)
+ errE = numpy.array( [0]*2, numpy.int32)
+
+ self._db.open()
+ self._db.queryVals(queryVals)
+ (data, err) = self._db.query(locs, cs, 3)
+ data = numpy.array(data)
+ err = numpy.array(err)
+
+ self.assertEqual(len(errE.shape), len(err.shape))
+ for dE, d in zip(errE.shape, err.shape):
+ self.assertEqual(dE, d)
+ for vE, v in zip(numpy.reshape(errE, -1), numpy.reshape(err, -1)):
+ self.assertEqual(vE, v)
+
+ self.assertEqual(len(dataE.shape), len(data.shape))
+ for dE, d in zip(dataE.shape, data.shape):
+ self.assertEqual(dE, d)
+ for vE, v in zip(numpy.reshape(dataE, -1), numpy.reshape(data, -1)):
+ self.assertAlmostEqual(vE, v, 6)
+
+ self._db.close()
+ return
+
+
+ def test_querySquash(self):
+ locs = numpy.array( [[-118.520000, 34.120000, -1400.00],
+ [-116.400000, 32.340000, -1000.00]],
+ numpy.float64)
+ from spatialdata.geocoords.CSGeo import CSGeo
+ cs = CSGeo()
+ cs._configure()
+ cs.datumHoriz = "NAD27"
+ cs.datumVert = "mean sea level"
+ cs.ellipsoid = "clrk66"
+ cs.initialize()
+ queryVals = ["topo-elev", "moho-depth", "density"]
+ dataE = numpy.array([[489.975189, -31178.105469, 2660.851074],
+ [801.209961, -34526.414062, 3022.192708]],
+ numpy.float64)
+ errE = numpy.array( [0]*2, numpy.int32)
+
+ self._db.open()
+ self._db.queryVals(queryVals)
+ self._db.squash = True
+ from pyre.units.length import km
+ self._db.squashLimit = -2.5*km
+ (data, err) = self._db.query(locs, cs, 3)
+ data = numpy.array(data)
+ err = numpy.array(err)
+
+ self.assertEqual(len(errE.shape), len(err.shape))
+ for dE, d in zip(errE.shape, err.shape):
+ self.assertEqual(dE, d)
+ for vE, v in zip(numpy.reshape(errE, -1), numpy.reshape(err, -1)):
+ self.assertEqual(vE, v)
+
+ self.assertEqual(len(dataE.shape), len(data.shape))
+ for dE, d in zip(dataE.shape, data.shape):
+ self.assertEqual(dE, d)
+ for vE, v in zip(numpy.reshape(dataE, -1), numpy.reshape(data, -1)):
+ self.assertAlmostEqual(vE, v, 6)
+
+ self._db.close()
+ return
+
+
+# End of file
Modified: cs/spatialdata-0.1/trunk/tests/pytests/spatialdb/TestSCECCVMH.py.in
===================================================================
--- cs/spatialdata-0.1/trunk/tests/pytests/spatialdb/TestSCECCVMH.py.in 2008-01-03 22:08:37 UTC (rev 8982)
+++ cs/spatialdata-0.1/trunk/tests/pytests/spatialdb/TestSCECCVMH.py.in 2008-01-03 22:44:59 UTC (rev 8983)
@@ -65,4 +65,46 @@
return
+ def test_querySquash(self):
+ locs = numpy.array( [[-118.520000, 34.120000, -1400.00],
+ [-116.400000, 32.340000, -1000.00]],
+ numpy.float64)
+ from spatialdata.geocoords.CSGeo import CSGeo
+ cs = CSGeo()
+ cs._configure()
+ cs.datumHoriz = "NAD27"
+ cs.datumVert = "mean sea level"
+ cs.ellipsoid = "clrk66"
+ cs.initialize()
+ queryVals = ["topo-elev", "moho-depth", "density"]
+ dataE = numpy.array([[489.975189, -31178.105469, 2660.851074],
+ [801.209961, -34526.414062, 3022.192708]],
+ numpy.float64)
+ errE = numpy.array( [0]*2, numpy.int32)
+
+ self._db.open()
+ self._db.queryVals(queryVals)
+ self._db.squash = True
+ from pyre.units.length import km
+ self._db.squashLimit = -2.5*km
+ (data, err) = self._db.query(locs, cs, 3)
+ data = numpy.array(data)
+ err = numpy.array(err)
+
+ self.assertEqual(len(errE.shape), len(err.shape))
+ for dE, d in zip(errE.shape, err.shape):
+ self.assertEqual(dE, d)
+ for vE, v in zip(numpy.reshape(errE, -1), numpy.reshape(err, -1)):
+ self.assertEqual(vE, v)
+
+ self.assertEqual(len(dataE.shape), len(data.shape))
+ for dE, d in zip(dataE.shape, data.shape):
+ self.assertEqual(dE, d)
+ for vE, v in zip(numpy.reshape(dataE, -1), numpy.reshape(data, -1)):
+ self.assertAlmostEqual(vE, v, 6)
+
+ self._db.close()
+ return
+
+
# End of file
More information about the cig-commits
mailing list