[cig-commits] r8302 - cs/spatialdata-0.1/trunk/libsrc/spatialdb

[brad at geodynamics.org]

Author: brad
Date: 2007-11-18 20:15:16 -0800 (Sun, 18 Nov 2007)
New Revision: 8302

Modified:
   cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc
   cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc
Log:
More work on SCEC CVM-H implementation. Fixed coordinate transformation bug and geometry bug. Still need to implement queryVals() and add unit tests.

Modified: cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc	2007-11-18 23:59:14 UTC (rev 8301)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/GocadVoxet.cc	2007-11-19 04:15:16 UTC (rev 8302)
@@ -23,8 +23,6 @@
 #include <sstream> // USES std::ostringsgream
 #include <assert.h> // USES assert()
 
-#include <iostream> // TEMPORARY
-
 #if defined(WORDS_BIGENDIAN)
 #define NATIVE_BIG_ENDIAN
 #else
@@ -84,24 +82,6 @@
       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;
@@ -244,9 +224,9 @@
   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;
+  _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
 
 // ----------------------------------------------------------------------

Modified: cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc
===================================================================
--- cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc	2007-11-18 23:59:14 UTC (rev 8301)
+++ cs/spatialdata-0.1/trunk/libsrc/spatialdb/SCECCVMH.cc	2007-11-19 04:15:16 UTC (rev 8302)
@@ -49,8 +49,9 @@
   proj.units("m");
   proj.projOptions("+zone=11");
 
-  _csUTM->datumHoriz("WGS84");
+  _csUTM->datumHoriz("NAD27");
   _csUTM->datumVert("mean sea level");
+  _csUTM->ellipsoid("clrk66");
   _csUTM->projector(proj);
   _csUTM->initialize();
 } // constructor
@@ -170,16 +171,16 @@
   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) {
@@ -196,17 +197,21 @@
       density = 1000.0;
   } else {
     outsideVoxet = _crustMantleVp->query(&vp, _xyzUTM);
-    if (!outsideVoxet)
+    if (!outsideVoxet) {
       outsideVoxet = _crustMantleVs->query(&vs, _xyzUTM);
-    density = vp / 3.0 + 1280.0;
+      density = vp / 3.0 + 1280.0;
+    } else {
+      vs = vp; // No data value
+      density = vp; // No data value
+    } // else
   } // 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;
+  std::cout << "elev: " << elev
+	    << ", basement depth: " << baseDepth
+	    << ", moho depth: " << mohoDepth
+	    << ", vp: " << vp
+	    << ", vs: " << vs
+	    << ", density: " << density << std::endl;
 } // query