Adding support for Generic EPSG codes via GDAL

cxx/isce3/core/Projections.cpp

@@ -364,6 +364,87 @@ int CEA::inverse(const Vec3& enu, Vec3& llh) const
     return 0;
 }
 
+GenericEPSG::GenericEPSG(int epsg) : ProjectionBase(epsg)
+{
+    // Initiate SRS with given EPSG code
+    if (_srs.importFromEPSG(epsg)) {
+        std::string errstr =
+                "In GenericEPSG::GenericEPSG - Invalid EPSG Code. Received ";
+        errstr += std::to_string(epsg);
+        errstr += ".";
+        throw std::invalid_argument(errstr);
+    }
+    _srs.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+
+    // Initialize 4979 for latlong on WGS84
+    if (_srs_llh.importFromEPSG(4979)) {
+        std::string errstr = 
+            "In Generic EPSG::GenericEPSG - Failed to initialize EPSG:4979.";
+        throw std::invalid_argument(errstr);
+    }
+    _srs_llh.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+
+    // Setup forward transform
+    _poct_fwd = OGRCreateCoordinateTransformation(&_srs_llh, &_srs);
+    if (_poct_fwd == nullptr)
+    {
+        std::string errstr = 
+            "In GenericEPSG::GeneticEPSG. Failed to initialize fwd transform. EPSG: ";
+        errstr += std::to_string(epsg);
+        errstr += ".";
+        throw std::invalid_argument(errstr);
+    }
+
+    // Setup inverse transform
+    _poct_inv = OGRCreateCoordinateTransformation(&_srs, &_srs_llh);
+    if (_poct_inv == nullptr)
+    {
+        std::string errstr = 
+            "In GenericEPSG::GeneticEPSG. Failed to initialize inv transform. EPSG: ";
+        errstr += std::to_string(epsg);
+        errstr += ".";
+        throw std::invalid_argument(errstr);
+    }
+}
+
+int GenericEPSG::forward(const Vec3& llh, Vec3& enu) const
+{
+    double x, y, z;
+    x = 180.0 * llh[0] / M_PI;
+    y = 180.0 * llh[1] / M_PI;
+    z = llh[2];
+
+    if (_poct_fwd->Transform(1, &x, &y, &z))
+    {
+        enu[0] = x;
+        enu[1] = y;
+        enu[2] = z;
+
+        return 0;
+    }
+    return 1;
+}
+
+int GenericEPSG::inverse(const Vec3& enu, Vec3& llh) const
+{
+    double x, y, z;
+    x = enu[0];
+    y = enu[1];
+    z = enu[2];
+
+    if (_poct_inv->Transform(1, &x, &y, &z))
+    {
+        llh[0] = x * M_PI / 180.0;
+        llh[1] = y * M_PI / 180.0;
+        llh[2] = z;
+
+        return 0;
+    }
+    return 1;
+}
+
+
+
 ProjectionBase* createProj(int epsgcode)
 {
     // Check for Lat/Lon
@@ -386,9 +467,13 @@ ProjectionBase* createProj(int epsgcode)
     else if (epsgcode == 6933) {
         return new CEA;
     } else {
-        throw isce3::except::RuntimeError(ISCE_SRCINFO(),
-                                          "Unknown EPSG code (in factory): " +
-                                                  std::to_string(epsgcode));
+        try {
+                return new GenericEPSG {epsgcode};
+        } catch (std::invalid_argument) {
+            throw isce3::except::RuntimeError(ISCE_SRCINFO(), 
+                    "Unknown EPSG code (in factory): " +
+                        std::to_string(epsgcode));
+        }
     }
 }
 

cxx/isce3/core/Projections.h

@@ -3,6 +3,8 @@
 #include <iostream>
 #include <memory>
 
+#include "ogr_spatialref.h"
+
 #include "Constants.h"
 #include "Ellipsoid.h"
 #include "Vector.h"
@@ -256,6 +258,47 @@ inline void CEA::print() const
               << "EPSG: " << code() << std::endl;
 }
 
+/**
+ * Generic EPSG code based extension of ProjBase
+ *
+ * Interpretation of EPSG code is handed off to GDAL.
+ */
+class GenericEPSG : public ProjectionBase {
+    // wkt corresponding to the EPSG code
+    OGRSpatialReference _srs;
+    OGRSpatialReference _srs_llh;
+    OGRCoordinateTransformation * _poct_fwd = nullptr;
+    OGRCoordinateTransformation * _poct_inv = nullptr;
+
+public:
+    GenericEPSG(int);
+
+    /** @copydoc ProjectionBase::print() */
+    void print() const override;
+
+    /** Transform from llh (rad) to xyz - units depend on EPSG */
+    int forward(const Vec3& llh, Vec3& xyz) const override;
+
+    /** Transform from xyz in EPSG-dependent units to llh (rad) */
+    int inverse(const Vec3& xyz, Vec3& llh) const override;
+
+    /** Destructor **/
+    ~GenericEPSG() {
+        if (_poct_fwd) {
+            delete _poct_fwd;
+        }
+        if (_poct_inv) {
+            delete _poct_inv;
+        }
+    }
+};
+
+inline void GenericEPSG::print() const
+{
+    std::cout << "Generic EPSG: " <<  code() << std::endl
+              << "WKT: " << _srs.exportToWkt() << std::endl;
+}
+
 // This is to create a projection system from the EPSG code
 ProjectionBase* createProj(int epsg);
 

tests/cxx/isce3/Sources.cmake

@@ -20,6 +20,7 @@ core/projections/cea.cpp
 core/projections/geocent.cpp
 core/projections/polar.cpp
 core/projections/utm.cpp
+core/projections/genericepsg.cpp
 core/serialization/serializeAttitude.cpp
 core/serialization/serializeDoppler.cpp
 core/serialization/serializeOrbit.cpp

tests/cxx/isce3/core/projections/genericepsg.cpp

@@ -0,0 +1,81 @@
+#include <gtest/gtest.h>
+#include <isce3/core/Projections.h>
+#include "projtest.h"
+
+isce3::core::GenericEPSG EU(3035);
+isce3::core::GenericEPSG Aus(3577);
+isce3::core::GenericEPSG Kor(5186);
+
+struct GenericEPSGTest : public ::testing::Test {
+    unsigned int fails;
+    virtual void SetUp() { fails = 0; }
+    virtual void TearDown() {
+        if (fails > 0) {
+            std::cerr << "GenericEPSG::TearDown sees failures" << std::endl;
+        }
+    }
+};
+
+#define epsgTest(...) PROJ_TEST(GenericEPSGTest, __VA_ARGS__)
+#define makeprojTest(...) MAKE_PROJ_TEST(GenericEPSGTest, __VA_ARGS__)
+
+epsgTest(EU, EUOrigin,
+        {10.0*M_PI/180.0, 52.0*M_PI/180.0, 0.},
+        {4321000., 3210000., 0.});
+epsgTest(Aus, AusOrigin,
+        {132.0*M_PI/180.0, 0., 0.},
+        {0.,0.,0.});
+epsgTest(Kor, KorOrigin,
+        {127.0*M_PI/180.0, 38.0*M_PI/180.0, 0.},
+        {200000.0, 600000.0, 0.0});
+
+//EU tests - GDAL transform only appears to be good to mm
+//Probably due to large x0 and y0?
+epsgTest(EU, EU1,
+        {0.07716832314413966, 1.0295281464986623, 0.},
+        {4.0e+06, 4.0e+06, 0.}, 5.0e-4);
+
+epsgTest(EU, EU2,
+        {-0.053932484569879118874, 1.052410267404324750729, 0.},
+        {3.6e+06, 4.2e+06, 0.}, 5.0e-4);
+
+epsgTest(EU, EU3,
+        {0.19684437483694725, 0.9843154958421356, 0.},
+        {4.4e+06, 3.7e+06, 0.}, 5.0e-4);
+
+//Australia tests
+epsgTest(Aus, Aus1,
+        {2.326601641674536, -0.59163508189832, 0.},
+        {1.2e+05, -3.7e+06, 0.});
+
+epsgTest(Aus, Aus2,
+        {2.313562879599467, -0.27900251672315923, 0.},
+        {6.0e+04, -1.7e+06, 0.});
+
+epsgTest(Aus, Aus3,
+        {2.344362841755795, -0.7993597057678191, 0.},
+        {1.9e+05, -5e+06, 0.});
+
+//South Korea tests - This only appears to be good to 1cm
+epsgTest(Kor, Kor1,
+        {2.2412392328140913, 0.6001667563213443, 0.},
+        {3.3e+05, 2.0e+05, 0.});
+
+epsgTest(Kor, Kor2,
+        {2.135176403804947, 0.6930318564897303, 0.},
+        {-2.0e+05, 8.0e+05, 0.});
+
+epsgTest(Kor, Kor3,
+        {2.103517202529827, 0.6836580237834484, 0.},
+        {-3.6e+05, 7.5e+05, 0.});
+
+// Test creating projections
+makeprojTest(3035, makeEU);
+makeprojTest(3577, makeAustralia);
+makeprojTest(5176, makeSouthKorea);
+
+int main(int argc, char **argv) {
+
+    ::testing::InitGoogleTest(&argc, argv);
+    return RUN_ALL_TESTS();
+}

tests/cxx/isce3/core/projections/projtest.h

@@ -10,17 +10,18 @@
 using isce3::core::ProjectionBase;
 using isce3::core::Vec3;
 
-auto projTest(const ProjectionBase& p, const Vec3& ref_llh, const Vec3& ref_xyz)
+auto projTest(const ProjectionBase& p, const Vec3& ref_llh, const Vec3& ref_xyz, double tol = 1.0e-6)
 {
     Vec3 xyz = p.forward(ref_llh);
-    EXPECT_NEAR(xyz[0], ref_xyz[0], 1e-6);
-    EXPECT_NEAR(xyz[1], ref_xyz[1], 1e-6);
-    EXPECT_NEAR(xyz[2], ref_xyz[2], 1e-6);
+    EXPECT_NEAR(xyz[0], ref_xyz[0], tol);
+    EXPECT_NEAR(xyz[1], ref_xyz[1], tol);
+    EXPECT_NEAR(xyz[2], ref_xyz[2], tol);
 
     Vec3 llh = p.inverse(ref_xyz);
     EXPECT_NEAR(llh[0], ref_llh[0], 1e-9);
     EXPECT_NEAR(llh[1], ref_llh[1], 1e-9);
     EXPECT_NEAR(llh[2], ref_llh[2], 1e-6);
+
 }
 
 #define PROJ_TEST(testclass, proj, name, ...)                                  \
@@ -30,3 +31,18 @@ auto projTest(const ProjectionBase& p, const Vec3& ref_llh, const Vec3& ref_xyz)
         fails += ::testing::Test::HasFailure();                                \
     }                                                                          \
     struct consume_semicolon
+
+auto makeProjTest(const int code)
+{
+    ProjectionBase *p = isce3::core::createProj(code);
+    ASSERT_EQ(p->code(), code);
+    delete p;
+}
+
+#define MAKE_PROJ_TEST(testclass, code, name)                                 \
+    TEST_F(testclass, name)                                                   \
+    {                                                                         \
+        makeProjTest(code);                                                   \
+        fails += ::testing::Test::HasFailure();                               \
+    }                                                                         \
+    struct consume_semicolon