Merge pull request #95 from baagaard-usgs/fix-coordinate-system-type-detection

Fix coordinate system type detection in CSGeo for computing surface normal

Author: baagaard-usgs

CHANGES.md

@@ -1,3 +1,9 @@
+## Main branch
+
+* **Fixed**
+  * Detect Proj string coordinate system type when computing surface normal. Default to (0, 0, 1) if coordinate system type is not recognized.
+  * Add missing include directive for cstddef to CoordSys.hh.
+
 ## Version 3.1.0 (2023/12/13)
 
 * Add `AnalyticDB` for a spatial database composed of analytic functions.

libsrc/spatialdata/geocoords/CSGeo.cc

@@ -121,6 +121,7 @@ spatialdata::geocoords::CSGeo::computeSurfaceNormal(double* dir,
         case PJ_TYPE_GEOGRAPHIC_3D_CRS:
         case PJ_TYPE_GEODETIC_CRS:
         case PJ_TYPE_PROJECTED_CRS:
+        case PJ_TYPE_OTHER_COORDINATE_OPERATION:
             for (size_t i = 0; i < numLocs; ++i) {
                 dir[i*numDims+0] = +0.0;
                 dir[i*numDims+1] = +0.0;
@@ -146,9 +147,14 @@ spatialdata::geocoords::CSGeo::computeSurfaceNormal(double* dir,
             break;
         } // PJ_TYPE_GEOCENTRIC_CRS
         default: {
-            std::ostringstream msg;
-            msg << "Unknown coordinate system type (" << projType << ") for coordinate system '" << _string << "'.";
-            throw std::logic_error(msg.str());
+            std::cout << "Internal error: Coordinate system type (" << projType
+                      << ") not recognized for coordinate system '" << _string << "' "
+                      << "when computing normal of ground surface. Using default value of (0, 0, +1).";
+            for (size_t i = 0; i < numLocs; ++i) {
+                dir[i*numDims+0] = +0.0;
+                dir[i*numDims+1] = +0.0;
+                dir[i*numDims+2] = +1.0;
+            } // for
         } // default
 
         } // switch

tests/libtests/geocoords/TestCSGeo.cc

@@ -152,6 +152,28 @@ spatialdata::geocoords::TestCSGeo::testComputeSurfaceNormal(void) {
         delete[] dirs;dirs = 0;
     } // UTM zone 10
 
+    { // Transverse mercator projection
+        cs.setString("+proj=tmerc +datum=WGS84 +lon_0=-100.0 +lat_0=+40.0 +k=0.9996 +units=km");
+        const size_t numLocs = 4;
+        const size_t numDims = 3;
+        const size_t size = numLocs * numDims;
+        const double coords[size] = {
+            57.0, 41.0, -1.0,
+            40.0, 40.0,  2.0,
+            30.0, 42.0,  3.0,
+            55.0, 41.0,  4.0,
+        };
+        double* dirs = new double[size];
+        cs.setSpaceDim(numDims);
+        cs.computeSurfaceNormal(dirs, coords, numLocs, numDims);
+        for (size_t iLoc = 0, i = 0; iLoc < numLocs; ++iLoc) {
+            CHECK(0.0 == dirs[i++]);
+            CHECK(0.0 == dirs[i++]);
+            CHECK(1.0 == dirs[i++]);
+        } // for
+        delete[] dirs;dirs = 0;
+    } // Transverse mercator projection
+
     { // Geocentric ECEF
         cs.setString("EPSG:4978");