From 38d47dd7cbb1bb793b6b1a285bfa511490ceddd1 Mon Sep 17 00:00:00 2001 From: Vissarion Fisikopoulos Date: Fri, 20 May 2022 17:21:06 +0300 Subject: [PATCH] [srs] Support axis orientation using the +axis proj4 argument --- .../boost/geometry/srs/projections/dpar.hpp | 17 ++++ .../srs/projections/impl/pj_datum_set.hpp | 8 +- .../geometry/srs/projections/impl/pj_fwd.hpp | 10 +- .../geometry/srs/projections/impl/pj_init.hpp | 56 ++++++++++++ .../geometry/srs/projections/impl/pj_inv.hpp | 18 +++- .../srs/projections/impl/projects.hpp | 6 +- .../geometry/srs/projections/par_data.hpp | 91 +++++++++++++++++++ .../boost/geometry/srs/projections/spar.hpp | 11 +++ test/srs/projections.cpp | 16 +++- 9 files changed, 221 insertions(+), 12 deletions(-) diff --git a/include/boost/geometry/srs/projections/dpar.hpp b/include/boost/geometry/srs/projections/dpar.hpp index fd799088d..3720669f6 100644 --- a/include/boost/geometry/srs/projections/dpar.hpp +++ b/include/boost/geometry/srs/projections/dpar.hpp @@ -487,6 +487,11 @@ enum name_units vunits }; +enum name_axis +{ + axis = 86 // 3 element list of numbers +}; + template struct parameter { @@ -660,6 +665,17 @@ struct parameter } #endif + parameter(name_axis id, std::initializer_list v) + : m_id(id) + , m_value(srs::detail::axis(v)) + { + std::size_t n = v.size(); + if (n != 3) + { + BOOST_THROW_EXCEPTION( projection_exception("Invalid number of axis elements. Should be 3.") ); + } + } + parameter(name_units id, value_units v) : m_id(id), m_value(int(v)) {} @@ -694,6 +710,7 @@ public: bool is_id_equal(name_sweep const& id) const { return m_id == int(id); } bool is_id_equal(name_towgs84 const& id) const { return m_id == int(id); } bool is_id_equal(name_units const& id) const { return m_id == int(id); } + bool is_id_equal(name_axis const& id) const { return m_id == int(id); } template V const& get_value() const diff --git a/include/boost/geometry/srs/projections/impl/pj_datum_set.hpp b/include/boost/geometry/srs/projections/impl/pj_datum_set.hpp index e44491e25..3f7e7ac82 100644 --- a/include/boost/geometry/srs/projections/impl/pj_datum_set.hpp +++ b/include/boost/geometry/srs/projections/impl/pj_datum_set.hpp @@ -103,7 +103,7 @@ inline const pj_datums_type* pj_datum_find_datum(srs::dpar::parameters con { typename srs::dpar::parameters::const_iterator it = pj_param_find(params, srs::dpar::datum); - + if (it != params.end()) { const pj_datums_type* pj_datums = pj_get_datums().first; @@ -184,7 +184,7 @@ inline bool pj_datum_find_nadgrids(srs::detail::proj4_parameters const& params, { std::string::size_type end = snadgrids.find(',', i); std::string name = snadgrids.substr(i, end - i); - + i = end; if (end != std::string::npos) ++i; @@ -207,7 +207,7 @@ inline bool pj_datum_find_nadgrids(srs::dpar::parameters const& params, { out = it->template get_value(); } - + return ! out.empty(); } @@ -284,7 +284,7 @@ inline bool pj_datum_find_towgs84(srs::dpar::parameters const& params, { typename srs::dpar::parameters::const_iterator it = pj_param_find(params, srs::dpar::towgs84); - + if (it != params.end()) { srs::detail::towgs84 const& diff --git a/include/boost/geometry/srs/projections/impl/pj_fwd.hpp b/include/boost/geometry/srs/projections/impl/pj_fwd.hpp index 01b7b97a2..786be4bc0 100644 --- a/include/boost/geometry/srs/projections/impl/pj_fwd.hpp +++ b/include/boost/geometry/srs/projections/impl/pj_fwd.hpp @@ -92,8 +92,14 @@ inline void pj_fwd(Prj const& prj, P const& par, LL const& ll, XY& xy) prj.fwd(par, lp_lon, lp_lat, x, y); - geometry::set<0>(xy, par.fr_meter * (par.a * x + par.x0)); - geometry::set<1>(xy, par.fr_meter * (par.a * y + par.y0)); + if (par.axis[0] == 0) + { + geometry::set<0>(xy, par.sign[0] * par.fr_meter * (par.a * x + par.x0)); + geometry::set<1>(xy, par.sign[1] * par.fr_meter * (par.a * y + par.y0)); + } else { + geometry::set<1>(xy, par.sign[1] * par.fr_meter * (par.a * x + par.x0)); + geometry::set<0>(xy, par.sign[0] * par.fr_meter * (par.a * y + par.y0)); + } } } // namespace detail diff --git a/include/boost/geometry/srs/projections/impl/pj_init.hpp b/include/boost/geometry/srs/projections/impl/pj_init.hpp index 2e36d639c..a17f3a22d 100644 --- a/include/boost/geometry/srs/projections/impl/pj_init.hpp +++ b/include/boost/geometry/srs/projections/impl/pj_init.hpp @@ -427,6 +427,59 @@ inline void pj_init_pm(srs::spar::parameters const& params, T& val) >::apply(params, val); } +/************************************************************************/ +/* pj_init_axis() */ +/************************************************************************/ + +template +inline void pj_init_axis(Params const& params, parameters & projdef) +{ + std::string axis = pj_get_param_s(params, "axis"); + if(! axis.empty()) + { + //projdef.axis.assign(axis.begin(), axis.end()); + for (std::size_t i = 0; i < axis.length(); ++i) + { + switch(axis[i]) + { + case 'w': + projdef.sign[i] = -1; + projdef.axis[i] = 0; + break; + case 'e': + projdef.sign[i] = 1; + projdef.axis[i] = 0; + break; + case 's': + projdef.sign[i] = -1; + projdef.axis[i] = 1; + break; + case 'n': + projdef.sign[i] = 1; + projdef.axis[i] = 1; + break; + case 'd': + projdef.sign[i] = -1; + projdef.axis[i] = 2; + break; + case 'u': + projdef.sign[i] = 1; + projdef.axis[i] = 2; + break; + default: + BOOST_THROW_EXCEPTION( projection_exception(error_axis) ); + } + } + // Currently not support elevation + if (projdef.axis[0] + projdef.axis[1] != 1) + { + BOOST_THROW_EXCEPTION( projection_exception(error_axis) ); + } + } + +} + + /************************************************************************/ /* pj_init() */ /* */ @@ -515,6 +568,9 @@ inline parameters pj_init(Params const& params) /* prime meridian */ pj_init_pm(params, pin.from_greenwich); + /* set axis orientation */ + pj_init_axis(params, pin); + return pin; } diff --git a/include/boost/geometry/srs/projections/impl/pj_inv.hpp b/include/boost/geometry/srs/projections/impl/pj_inv.hpp index aa65af974..0a3d55cf0 100644 --- a/include/boost/geometry/srs/projections/impl/pj_inv.hpp +++ b/include/boost/geometry/srs/projections/impl/pj_inv.hpp @@ -62,12 +62,22 @@ inline void pj_inv(PRJ const& prj, PAR const& par, XY const& xy, LL& ll) /* can't do as much preliminary checking as with forward */ /* descale and de-offset */ - calc_t xy_x = (geometry::get<0>(xy) * par.to_meter - par.x0) * par.ra; - calc_t xy_y = (geometry::get<1>(xy) * par.to_meter - par.y0) * par.ra; - calc_t lon = 0, lat = 0; + calc_t lon = 0; + calc_t lat = 0; + calc_t xy_x = 0; + calc_t xy_y = 0; + + if (par.axis[0] == 1) + { + xy_x = (geometry::get<1>(xy) * par.to_meter * par.sign[1] - par.x0) * par.ra; + xy_y = (geometry::get<0>(xy) * par.to_meter * par.sign[0] - par.y0) * par.ra; + } else { + xy_x = (geometry::get<0>(xy) * par.to_meter * par.sign[0] - par.x0) * par.ra; + xy_y = (geometry::get<1>(xy) * par.to_meter * par.sign[1] - par.y0) * par.ra; + } prj.inv(par, xy_x, xy_y, lon, lat); /* inverse project */ - + lon += par.lam0; /* reduce from del lp.lam */ if (!par.over) lon = adjlon(lon); /* adjust longitude to CM */ diff --git a/include/boost/geometry/srs/projections/impl/projects.hpp b/include/boost/geometry/srs/projections/impl/projects.hpp index 6ca757252..a3b8a64fb 100644 --- a/include/boost/geometry/srs/projections/impl/projects.hpp +++ b/include/boost/geometry/srs/projections/impl/projects.hpp @@ -106,7 +106,7 @@ struct pj_consts T to_meter, fr_meter; /* cartesian scaling */ T vto_meter, vfr_meter; /* Vertical scaling. Internal unit [m] */ - // D A T U M S A N D H E I G H T S Y S T E M S + // D A T U M S A N D H E I G H T S Y S T E M S T from_greenwich; /* prime meridian offset (in radians) */ T long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/ @@ -128,6 +128,9 @@ struct pj_consts //enum pj_io_units left; /* Flags for input/output coordinate types */ //enum pj_io_units right; + srs::detail::axis axis; + srs::detail::axis sign; + // Initialize all variables pj_consts() : a(0), ra(0) @@ -140,6 +143,7 @@ struct pj_consts , datum_type(datum_unknown) , is_long_wrap_set(false) , over(false), geoc(false), is_latlong(false), is_geocent(false) + , axis(0,1,2), sign(1,1,1) //the default east, northing, elevation //, need_ellps(true) //, left(PJ_IO_UNITS_ANGULAR), right(PJ_IO_UNITS_CLASSIC) {} diff --git a/include/boost/geometry/srs/projections/par_data.hpp b/include/boost/geometry/srs/projections/par_data.hpp index c86767612..1af67e374 100644 --- a/include/boost/geometry/srs/projections/par_data.hpp +++ b/include/boost/geometry/srs/projections/par_data.hpp @@ -188,6 +188,97 @@ private: T m_data[7]; }; +struct axis +{ + static const std::size_t static_capacity = 3; + + typedef std::size_t size_type; + typedef int value_type; + typedef int* iterator; + typedef const int* const_iterator; + typedef int& reference; + typedef const int& const_reference; + + axis() + : m_size(3) + , m_data{0, 0, 0} + {} + + template + axis(It first, It last) + { + assign(first, last); + } + + axis(std::initializer_list l) + { + assign(l.begin(), l.end()); + } + + axis(int const& v0, int const& v1, int const& v2) + : m_size(3) + { + m_data[0] = v0; + m_data[1] = v1; + m_data[2] = v2; + } + + void push_back(int const& v) + { + BOOST_GEOMETRY_ASSERT(m_size < static_capacity); + m_data[m_size] = v; + ++m_size; + } + + template + void assign(It first, It last) + { + for (m_size = 0 ; first != last && m_size < 3 ; ++first, ++m_size) + m_data[m_size] = *first; + } + + void assign(std::initializer_list l) + { + assign(l.begin(), l.end()); + } + + const_reference operator[](size_type i) const + { + BOOST_GEOMETRY_ASSERT(i < m_size); + return m_data[i]; + } + + reference operator[](size_type i) + { + BOOST_GEOMETRY_ASSERT(i < m_size); + return m_data[i]; + } + + size_type size() const + { + return m_size; + } + + bool empty() const + { + return m_size == 0; + } + + void clear() + { + m_size = 0; + } + + iterator begin() { return m_data; } + iterator end() { return m_data + m_size; } + const_iterator begin() const { return m_data; } + const_iterator end() const { return m_data + m_size; } + +private: + size_type m_size; + int m_data[3]; +}; + } // namespace detail #endif // DOXYGEN_NO_DETAIL diff --git a/include/boost/geometry/srs/projections/spar.hpp b/include/boost/geometry/srs/projections/spar.hpp index ef15e1544..30b7dfb05 100644 --- a/include/boost/geometry/srs/projections/spar.hpp +++ b/include/boost/geometry/srs/projections/spar.hpp @@ -493,6 +493,17 @@ struct towgs84 #endif }; +struct axis + : srs::detail::axis +{ + typedef srs::detail::axis base_t; + + axis(int const& v0, int const& v1, int const& v2) + : base_t(v0, v1, v2) + {} + axis(std::initializer_list l) : base_t(l) {} +}; + template struct vunits { diff --git a/test/srs/projections.cpp b/test/srs/projections.cpp index eb36aebf5..90afa3635 100644 --- a/test/srs/projections.cpp +++ b/test/srs/projections.cpp @@ -413,7 +413,21 @@ void test_srs() // EPSG:9833 Hyperbolic Cassini-Soldner, SRS EPSG:3139 test_both

("cass", 179.99433651, -16.841456514, 322174, 268950, "+proj=cass +hyperbolic +lat_0=-16.25 +lon_0=179.33333332 +x_0=251727.9155424\ - +y_0=334519.953768 +towgs84=51,391,-36,0,0,0,0"); + +y_0=334519.953768 +towgs84=51,391,-36,0,0,0,0 +axis=enu"); + + // test +axis argument + test_both

("cass", 179.99433651, -16.841456514, -322174, 268950, + "+proj=cass +hyperbolic +lat_0=-16.25 +lon_0=179.33333332 +x_0=251727.9155424\ + +y_0=334519.953768 +towgs84=51,391,-36,0,0,0,0 +axis=wnu"); + test_both

("cass", 179.99433651, -16.841456514, 322174, -268950, + "+proj=cass +hyperbolic +lat_0=-16.25 +lon_0=179.33333332 +x_0=251727.9155424\ + +y_0=334519.953768 +towgs84=51,391,-36,0,0,0,0 +axis=esu"); + test_both

("cass", 179.99433651, -16.841456514, -322174, -268950, + "+proj=cass +hyperbolic +lat_0=-16.25 +lon_0=179.33333332 +x_0=251727.9155424\ + +y_0=334519.953768 +towgs84=51,391,-36,0,0,0,0 +axis=wsu"); + test_both

("cass", 179.99433651, -16.841456514, 268950, 322174, + "+proj=cass +hyperbolic +lat_0=-16.25 +lon_0=179.33333332 +x_0=251727.9155424\ + +y_0=334519.953768 +towgs84=51,391,-36,0,0,0,0 +axis=neu"); } int test_main(int, char* [])