boost_mirror/geometry
1
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2025-05-11 13:34:10 +00:00
Code Issues Packages Projects Releases Wiki Activity

[srs] Support axis orientation using the +axis proj4 argument

Browse Source
This commit is contained in:
Vissarion Fisikopoulos 2022-05-20 17:21:06 +03:00
parent 442d03cef5
commit 38d47dd7cb
9 changed files with 221 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
include/boost/geometry/srs/projections/dpar.hpp
View File

@ -487,6 +487,11 @@ enum name_units
vunits
};
enum name_axis
{
axis = 86 // 3 element list of numbers
};
template <typename T>
struct parameter
{
@ -660,6 +665,17 @@ struct parameter
}
#endif
parameter(name_axis id, std::initializer_list<int> 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 <typename V>
V const& get_value() const

8
include/boost/geometry/srs/projections/impl/pj_datum_set.hpp
View File

@ -103,7 +103,7 @@ inline const pj_datums_type<T>* pj_datum_find_datum(srs::dpar::parameters<T> con
{
typename srs::dpar::parameters<T>::const_iterator
it = pj_param_find(params, srs::dpar::datum);
if (it != params.end())
{
const pj_datums_type<T>* pj_datums = pj_get_datums<T>().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<T> const& params,
{
out = it->template get_value<srs::detail::nadgrids>();
}
return ! out.empty();
}
@ -284,7 +284,7 @@ inline bool pj_datum_find_towgs84(srs::dpar::parameters<T> const& params,
{
typename srs::dpar::parameters<T>::const_iterator
it = pj_param_find(params, srs::dpar::towgs84);
if (it != params.end())
{
srs::detail::towgs84<T> const&

10
include/boost/geometry/srs/projections/impl/pj_fwd.hpp
View File

@ -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

56
include/boost/geometry/srs/projections/impl/pj_init.hpp
View File

@ -427,6 +427,59 @@ inline void pj_init_pm(srs::spar::parameters<Ps...> const& params, T& val)
>::apply(params, val);
}
/************************************************************************/
/* pj_init_axis() */
/************************************************************************/
template <typename Params, typename T>
inline void pj_init_axis(Params const& params, parameters<T> & 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<T> pj_init(Params const& params)
/* prime meridian */
pj_init_pm(params, pin.from_greenwich);
/* set axis orientation */
pj_init_axis(params, pin);
return pin;
}

18
include/boost/geometry/srs/projections/impl/pj_inv.hpp
View File

@ -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 */

6
include/boost/geometry/srs/projections/impl/projects.hpp
View File

@ -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)
{}

91
include/boost/geometry/srs/projections/par_data.hpp
View File

@ -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 <typename It>
axis(It first, It last)
{
assign(first, last);
}
axis(std::initializer_list<int> 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 <typename It>
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<int> 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

11
include/boost/geometry/srs/projections/spar.hpp
View File

@ -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<int> l) : base_t(l) {}
};
template <typename Units>
struct vunits
{

16
test/srs/projections.cpp
View File

@ -413,7 +413,21 @@ void test_srs()
// EPSG:9833 Hyperbolic Cassini-Soldner, SRS EPSG:3139
test_both<P>("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<P>("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<P>("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<P>("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<P>("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* [])