diff --git a/cxx/isce3/core/Projections.cpp b/cxx/isce3/core/Projections.cpp index dcd539253..8769460d1 100644 --- a/cxx/isce3/core/Projections.cpp +++ b/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)); + } } } diff --git a/cxx/isce3/core/Projections.h b/cxx/isce3/core/Projections.h index cf5461b6d..187c98ec6 100644 --- a/cxx/isce3/core/Projections.h +++ b/cxx/isce3/core/Projections.h @@ -3,6 +3,8 @@ #include #include +#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); diff --git a/tests/cxx/isce3/Sources.cmake b/tests/cxx/isce3/Sources.cmake index 80c5f37cb..4f75d4b7f 100644 --- a/tests/cxx/isce3/Sources.cmake +++ b/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 diff --git a/tests/cxx/isce3/core/projections/genericepsg.cpp b/tests/cxx/isce3/core/projections/genericepsg.cpp new file mode 100644 index 000000000..a24ae9b89 --- /dev/null +++ b/tests/cxx/isce3/core/projections/genericepsg.cpp @@ -0,0 +1,81 @@ +#include +#include +#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(); +} diff --git a/tests/cxx/isce3/core/projections/projtest.h b/tests/cxx/isce3/core/projections/projtest.h index e18fe4ae1..d89776c26 100644 --- a/tests/cxx/isce3/core/projections/projtest.h +++ b/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