diff options
Diffstat (limited to 'boost/geometry/srs/projections/proj/healpix.hpp')
-rw-r--r-- | boost/geometry/srs/projections/proj/healpix.hpp | 897 |
1 files changed, 897 insertions, 0 deletions
diff --git a/boost/geometry/srs/projections/proj/healpix.hpp b/boost/geometry/srs/projections/proj/healpix.hpp new file mode 100644 index 0000000000..d2f5a8c081 --- /dev/null +++ b/boost/geometry/srs/projections/proj/healpix.hpp @@ -0,0 +1,897 @@ +#ifndef BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP +#define BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP + +// Boost.Geometry - extensions-gis-projections (based on PROJ4) +// This file is automatically generated. DO NOT EDIT. + +// Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017, 2018. +// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +// This file is converted from PROJ4, http://trac.osgeo.org/proj +// PROJ4 is originally written by Gerald Evenden (then of the USGS) +// PROJ4 is maintained by Frank Warmerdam +// PROJ4 is converted to Boost.Geometry by Barend Gehrels + +// Last updated version of proj: 4.9.1 + +// Original copyright notice: + +// Purpose: Implementation of the HEALPix and rHEALPix projections. +// For background see <http://code.scenzgrid.org/index.php/p/scenzgrid-py/source/tree/master/docs/rhealpix_dggs.pdf>. +// Authors: Alex Raichev (raichev@cs.auckland.ac.nz) +// Michael Speth (spethm@landcareresearch.co.nz) +// Notes: Raichev implemented these projections in Python and +// Speth translated them into C here. +// Copyright (c) 2001, Thomas Flemming, tf@ttqv.com + +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: + +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. + +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. + +#include <boost/geometry/util/math.hpp> + +#include <boost/geometry/srs/projections/impl/base_static.hpp> +#include <boost/geometry/srs/projections/impl/base_dynamic.hpp> +#include <boost/geometry/srs/projections/impl/projects.hpp> +#include <boost/geometry/srs/projections/impl/factory_entry.hpp> +#include <boost/geometry/srs/projections/impl/pj_auth.hpp> +#include <boost/geometry/srs/projections/impl/pj_qsfn.hpp> + +namespace boost { namespace geometry +{ + +namespace srs { namespace par4 +{ + struct healpix {}; + struct rhealpix {}; + +}} //namespace srs::par4 + +namespace projections +{ + #ifndef DOXYGEN_NO_DETAIL + namespace detail { namespace healpix + { + + static const double EPS = 1e-15; + + template <typename T> + struct par_healpix + { + int north_square; + int south_square; + T qp; + T apa[APA_SIZE]; + }; + + /* Matrix for counterclockwise rotation by pi/2: */ + /* Matrix for counterclockwise rotation by pi: */ + /* Matrix for counterclockwise rotation by 3*pi/2: */ + /* Identity matrix */ + /* IDENT, R1, R2, R3, R1 inverse, R2 inverse, R3 inverse:*/ + /* Fuzz to handle rounding errors: */ + template <typename T> + struct CapMap + { + int cn; /* An integer 0--3 indicating the position of the polar cap. */ + T x, y; /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */ + enum Region {north, south, equatorial} region; + }; + template <typename T> + struct Point + { + T x, y; + }; + static double rot[7][2][2] = {{{1, 0},{0, 1}}, {{ 0,-1},{ 1, 0}}, {{-1, 0},{ 0,-1}}, {{ 0, 1},{-1, 0}}, {{ 0, 1},{-1, 0}}, {{-1, 0},{ 0,-1}}, {{ 0,-1},{ 1, 0}}}; + + /** + * Returns the sign of the double. + * @param v the parameter whose sign is returned. + * @return 1 for positive number, -1 for negative, and 0 for zero. + **/ + template <typename T> + inline T pj_sign (T const& v) + { + return v > 0 ? 1 : (v < 0 ? -1 : 0); + } + /** + * Return the index of the matrix in {{{1, 0},{0, 1}}, {{ 0,-1},{ 1, 0}}, {{-1, 0},{ 0,-1}}, {{ 0, 1},{-1, 0}}, {{ 0, 1},{-1, 0}}, {{-1, 0},{ 0,-1}}, {{ 0,-1},{ 1, 0}}}. + * @param index ranges from -3 to 3. + */ + inline int get_rotate_index(int index) + { + switch(index) { + case 0: + return 0; + case 1: + return 1; + case 2: + return 2; + case 3: + return 3; + case -1: + return 4; + case -2: + return 5; + case -3: + return 6; + } + return 0; + } + /** + * Return 1 if point (testx, testy) lies in the interior of the polygon + * determined by the vertices in vert, and return 0 otherwise. + * See http://paulbourke.net/geometry/polygonmesh/ for more details. + * @param nvert the number of vertices in the polygon. + * @param vert the (x, y)-coordinates of the polygon's vertices + **/ + template <typename T> + inline int pnpoly(int nvert, T vert[][2], T const& testx, T const& testy) + { + int i, c = 0; + int counter = 0; + T xinters; + Point<T> p1, p2; + /* Check for boundrary cases */ + for (i = 0; i < nvert; i++) { + if (testx == vert[i][0] && testy == vert[i][1]) { + return 1; + } + } + p1.x = vert[0][0]; + p1.y = vert[0][1]; + for (i = 1; i < nvert; i++) { + p2.x = vert[i % nvert][0]; + p2.y = vert[i % nvert][1]; + if (testy > (std::min)(p1.y, p2.y)) { + if (testy <= (std::max)(p1.y, p2.y)) { + if (testx <= (std::max)(p1.x, p2.x)) { + if (p1.y != p2.y) { + xinters = (testy-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x; + if (p1.x == p2.x || testx <= xinters) { + counter++; + } + } + } + } + } + p1 = p2; + } + if (counter % 2 == 0) { + return 0; + } else { + return 1; + } + return c; + } + /** + * Return 1 if (x, y) lies in (the interior or boundary of) the image of the + * HEALPix projection (in case proj=0) or in the image the rHEALPix projection + * (in case proj=1), and return 0 otherwise. + * @param north_square the position of the north polar square (rHEALPix only) + * @param south_square the position of the south polar square (rHEALPix only) + **/ + template <typename T> + inline int in_image(T const& x, T const& y, int proj, int north_square, int south_square) + { + static const T ONEPI = detail::ONEPI<T>(); + + if (proj == 0) { + T healpixVertsJit[][2] = { + {-1.0*ONEPI- EPS, ONEPI/4.0}, + {-3.0*ONEPI/4.0, ONEPI/2.0 + EPS}, + {-1.0*ONEPI/2.0, ONEPI/4.0 + EPS}, + {-1.0*ONEPI/4.0, ONEPI/2.0 + EPS}, + {0.0, ONEPI/4.0 + EPS}, + {ONEPI/4.0, ONEPI/2.0 + EPS}, + {ONEPI/2.0, ONEPI/4.0 + EPS}, + {3.0*ONEPI/4.0, ONEPI/2.0 + EPS}, + {ONEPI+ EPS, ONEPI/4.0}, + {ONEPI+ EPS, -1.0*ONEPI/4.0}, + {3.0*ONEPI/4.0, -1.0*ONEPI/2.0 - EPS}, + {ONEPI/2.0, -1.0*ONEPI/4.0 - EPS}, + {ONEPI/4.0, -1.0*ONEPI/2.0 - EPS}, + {0.0, -1.0*ONEPI/4.0 - EPS}, + {-1.0*ONEPI/4.0, -1.0*ONEPI/2.0 - EPS}, + {-1.0*ONEPI/2.0, -1.0*ONEPI/4.0 - EPS}, + {-3.0*ONEPI/4.0, -1.0*ONEPI/2.0 - EPS}, + {-1.0*ONEPI - EPS, -1.0*ONEPI/4.0} + }; + return pnpoly((int)sizeof(healpixVertsJit)/ + sizeof(healpixVertsJit[0]), healpixVertsJit, x, y); + } else { + T rhealpixVertsJit[][2] = { + {-1.0*ONEPI - EPS, ONEPI/4.0 + EPS}, + {-1.0*ONEPI + north_square*ONEPI/2.0- EPS, ONEPI/4.0 + EPS}, + {-1.0*ONEPI + north_square*ONEPI/2.0- EPS, 3*ONEPI/4.0 + EPS}, + {-1.0*ONEPI + (north_square + 1.0)*ONEPI/2.0 + EPS, 3*ONEPI/4.0 + EPS}, + {-1.0*ONEPI + (north_square + 1.0)*ONEPI/2.0 + EPS, ONEPI/4.0 + EPS}, + {ONEPI + EPS, ONEPI/4.0 + EPS}, + {ONEPI + EPS, -1.0*ONEPI/4.0 - EPS}, + {-1.0*ONEPI + (south_square + 1.0)*ONEPI/2.0 + EPS, -1.0*ONEPI/4.0 - EPS}, + {-1.0*ONEPI + (south_square + 1.0)*ONEPI/2.0 + EPS, -3.0*ONEPI/4.0 - EPS}, + {-1.0*ONEPI + south_square*ONEPI/2.0 - EPS, -3.0*ONEPI/4.0 - EPS}, + {-1.0*ONEPI + south_square*ONEPI/2.0 - EPS, -1.0*ONEPI/4.0 - EPS}, + {-1.0*ONEPI - EPS, -1.0*ONEPI/4.0 - EPS}}; + return pnpoly((int)sizeof(rhealpixVertsJit)/ + sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y); + } + } + /** + * Return the authalic latitude of latitude alpha (if inverse=0) or + * return the approximate latitude of authalic latitude alpha (if inverse=1). + * P contains the relavent ellipsoid parameters. + **/ + template <typename Parameters, typename T> + inline T auth_lat(const Parameters& par, const par_healpix<T>& proj_parm, T const& alpha, int inverse) + { + if (inverse == 0) { + /* Authalic latitude. */ + T q = pj_qsfn(sin(alpha), par.e, 1.0 - par.es); + T qp = proj_parm.qp; + T ratio = q/qp; + if (fabsl(ratio) > 1) { + /* Rounding error. */ + ratio = pj_sign(ratio); + } + return asin(ratio); + } else { + /* Approximation to inverse authalic latitude. */ + return pj_authlat(alpha, proj_parm.apa); + } + } + /** + * Return the HEALPix projection of the longitude-latitude point lp on + * the unit sphere. + **/ + template <typename T> + inline void healpix_sphere(T const& lp_lam, T const& lp_phi, T& xy_x, T& xy_y) + { + static const T ONEPI = detail::ONEPI<T>(); + + T lam = lp_lam; + T phi = lp_phi; + T phi0 = asin(T(2.0)/T(3.0)); + + /* equatorial region */ + if ( fabsl(phi) <= phi0) { + xy_x = lam; + xy_y = 3.0*ONEPI/8.0*sin(phi); + } else { + T lamc; + T sigma = sqrt(3.0*(1 - fabsl(sin(phi)))); + T cn = floor(2*lam / ONEPI + 2); + if (cn >= 4) { + cn = 3; + } + lamc = -3*ONEPI/4 + (ONEPI/2)*cn; + xy_x = lamc + (lam - lamc)*sigma; + xy_y = pj_sign(phi)*ONEPI/4*(2 - sigma); + } + return; + } + /** + * Return the inverse of healpix_sphere(). + **/ + template <typename T> + inline void healpix_sphere_inverse(T const& xy_x, T const& xy_y, T& lp_lam, T& lp_phi) + { + static const T ONEPI = detail::ONEPI<T>(); + + T x = xy_x; + T y = xy_y; + T y0 = ONEPI/4.0; + /* Equatorial region. */ + if (fabsl(y) <= y0) { + lp_lam = x; + lp_phi = asin(8.0*y/(3.0*ONEPI)); + } else if (fabsl(y) < ONEPI/2.0) { + T cn = floor(2.0*x/ONEPI + 2.0); + T xc, tau; + if (cn >= 4) { + cn = 3; + } + xc = -3.0*ONEPI/4.0 + (ONEPI/2.0)*cn; + tau = 2.0 - 4.0*fabsl(y)/ONEPI; + lp_lam = xc + (x - xc)/tau; + lp_phi = pj_sign(y)*asin(1.0 - pow(tau , 2.0)/3.0); + } else { + lp_lam = -1.0*ONEPI; + lp_phi = pj_sign(y)*ONEPI/2.0; + } + return; + } + /** + * Return the vector sum a + b, where a and b are 2-dimensional vectors. + * @param ret holds a + b. + **/ + template <typename T> + inline void vector_add(T a[2], T b[2], T *ret) + { + int i; + for(i = 0; i < 2; i++) { + ret[i] = a[i] + b[i]; + } + } + /** + * Return the vector difference a - b, where a and b are 2-dimensional vectors. + * @param ret holds a - b. + **/ + template <typename T> + inline void vector_sub(T a[2], T b[2], T*ret) + { + int i; + for(i = 0; i < 2; i++) { + ret[i] = a[i] - b[i]; + } + } + /** + * Return the 2 x 1 matrix product a*b, where a is a 2 x 2 matrix and + * b is a 2 x 1 matrix. + * @param ret holds a*b. + **/ + template <typename T> + inline void dot_product(T a[2][2], T b[2], T *ret) + { + int i, j; + int length = 2; + for(i = 0; i < length; i++) { + ret[i] = 0; + for(j = 0; j < length; j++) { + ret[i] += a[i][j]*b[j]; + } + } + } + /** + * Return the number of the polar cap, the pole point coordinates, and + * the region that (x, y) lies in. + * If inverse=0, then assume (x,y) lies in the image of the HEALPix + * projection of the unit sphere. + * If inverse=1, then assume (x,y) lies in the image of the + * (north_square, south_square)-rHEALPix projection of the unit sphere. + **/ + template <typename T> + inline CapMap<T> get_cap(T x, T const& y, int north_square, int south_square, + int inverse) + { + static const T ONEPI = detail::ONEPI<T>(); + + CapMap<T> capmap; + T c; + capmap.x = x; + capmap.y = y; + if (inverse == 0) { + if (y > ONEPI/4.0) { + capmap.region = CapMap<T>::north; + c = ONEPI/2.0; + } else if (y < -1*ONEPI/4.0) { + capmap.region = CapMap<T>::south; + c = -1*ONEPI/2.0; + } else { + capmap.region = CapMap<T>::equatorial; + capmap.cn = 0; + return capmap; + } + /* polar region */ + if (x < -1*ONEPI/2.0) { + capmap.cn = 0; + capmap.x = (-1*3.0*ONEPI/4.0); + capmap.y = c; + } else if (x >= -1*ONEPI/2.0 && x < 0) { + capmap.cn = 1; + capmap.x = -1*ONEPI/4.0; + capmap.y = c; + } else if (x >= 0 && x < ONEPI/2.0) { + capmap.cn = 2; + capmap.x = ONEPI/4.0; + capmap.y = c; + } else { + capmap.cn = 3; + capmap.x = 3.0*ONEPI/4.0; + capmap.y = c; + } + return capmap; + } else { + T eps; + if (y > ONEPI/4.0) { + capmap.region = CapMap<T>::north; + capmap.x = (-3.0*ONEPI/4.0 + north_square*ONEPI/2.0); + capmap.y = ONEPI/2.0; + x = x - north_square*ONEPI/2.0; + } else if (y < -1*ONEPI/4.0) { + capmap.region = CapMap<T>::south; + capmap.x = (-3.0*ONEPI/4.0 + south_square*ONEPI/2); + capmap.y = -1*ONEPI/2.0; + x = x - south_square*ONEPI/2.0; + } else { + capmap.region = CapMap<T>::equatorial; + capmap.cn = 0; + return capmap; + } + /* Polar Region, find the HEALPix polar cap number that + x, y moves to when rHEALPix polar square is disassembled. */ + eps = 1e-15; /* Kludge. Fuzz to avoid some rounding errors. */ + if (capmap.region == CapMap<T>::north) { + if (y >= -1*x - ONEPI/4.0 - eps && y < x + 5.0*ONEPI/4.0 - eps) { + capmap.cn = (north_square + 1) % 4; + } else if (y > -1*x -1*ONEPI/4.0 + eps && y >= x + 5.0*ONEPI/4.0 - eps) { + capmap.cn = (north_square + 2) % 4; + } else if (y <= -1*x -1*ONEPI/4.0 + eps && y > x + 5.0*ONEPI/4.0 + eps) { + capmap.cn = (north_square + 3) % 4; + } else { + capmap.cn = north_square; + } + } else if (capmap.region == CapMap<T>::south) { + if (y <= x + ONEPI/4.0 + eps && y > -1*x - 5.0*ONEPI/4 + eps) { + capmap.cn = (south_square + 1) % 4; + } else if (y < x + ONEPI/4.0 - eps && y <= -1*x - 5.0*ONEPI/4.0 + eps) { + capmap.cn = (south_square + 2) % 4; + } else if (y >= x + ONEPI/4.0 - eps && y < -1*x - 5.0*ONEPI/4.0 - eps) { + capmap.cn = (south_square + 3) % 4; + } else { + capmap.cn = south_square; + } + } + return capmap; + } + } + /** + * Rearrange point (x, y) in the HEALPix projection by + * combining the polar caps into two polar squares. + * Put the north polar square in position north_square and + * the south polar square in position south_square. + * If inverse=1, then uncombine the polar caps. + * @param north_square integer between 0 and 3. + * @param south_square integer between 0 and 3. + **/ + template <typename T> + inline void combine_caps(T& xy_x, T& xy_y, int north_square, int south_square, + int inverse) + { + static const T ONEPI = detail::ONEPI<T>(); + + T v[2]; + T a[2]; + T vector[2]; + T v_min_c[2]; + T ret_dot[2]; + CapMap<T> capmap = get_cap(xy_x, xy_y, north_square, south_square, inverse); + if (capmap.region == CapMap<T>::equatorial) { + xy_x = capmap.x; + xy_y = capmap.y; + return; + } + v[0] = xy_x; + v[1] = xy_y; + if (inverse == 0) { + /* Rotate (xy_x, xy_y) about its polar cap tip and then translate it to + north_square or south_square. */ + int pole = 0; + T (*tmpRot)[2]; + T c[2] = {capmap.x, capmap.y}; + if (capmap.region == CapMap<T>::north) { + pole = north_square; + a[0] = (-3.0*ONEPI/4.0 + pole*ONEPI/2); + a[1] = (ONEPI/2.0 + pole*0); + tmpRot = rot[get_rotate_index(capmap.cn - pole)]; + vector_sub(v, c, v_min_c); + dot_product(tmpRot, v_min_c, ret_dot); + vector_add(ret_dot, a, vector); + } else { + pole = south_square; + a[0] = (-3.0*ONEPI/4.0 + pole*ONEPI/2); + a[1] = (ONEPI/-2.0 + pole*0); + tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))]; + vector_sub(v, c, v_min_c); + dot_product(tmpRot, v_min_c, ret_dot); + vector_add(ret_dot, a, vector); + } + xy_x = vector[0]; + xy_y = vector[1]; + return; + } else { + /* Inverse function. + Unrotate (xy_x, xy_y) and then translate it back. */ + int pole = 0; + T (*tmpRot)[2]; + T c[2] = {capmap.x, capmap.y}; + /* disassemble */ + if (capmap.region == CapMap<T>::north) { + pole = north_square; + a[0] = (-3.0*ONEPI/4.0 + capmap.cn*ONEPI/2); + a[1] = (ONEPI/2.0 + capmap.cn*0); + tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))]; + vector_sub(v, c, v_min_c); + dot_product(tmpRot, v_min_c, ret_dot); + vector_add(ret_dot, a, vector); + } else { + pole = south_square; + a[0] = (-3.0*ONEPI/4.0 + capmap.cn*ONEPI/2); + a[1] = (ONEPI/-2.0 + capmap.cn*0); + tmpRot = rot[get_rotate_index(capmap.cn - pole)]; + vector_sub(v, c, v_min_c); + dot_product(tmpRot, v_min_c, ret_dot); + vector_add(ret_dot, a, vector); + } + xy_x = vector[0]; + xy_y = vector[1]; + return; + } + } + + // template class, using CRTP to implement forward/inverse + template <typename CalculationType, typename Parameters> + struct base_healpix_ellipsoid : public base_t_fi<base_healpix_ellipsoid<CalculationType, Parameters>, + CalculationType, Parameters> + { + + typedef CalculationType geographic_type; + typedef CalculationType cartesian_type; + + par_healpix<CalculationType> m_proj_parm; + + inline base_healpix_ellipsoid(const Parameters& par) + : base_t_fi<base_healpix_ellipsoid<CalculationType, Parameters>, + CalculationType, Parameters>(*this, par) {} + + // FORWARD(e_healpix_forward) ellipsoid + // Project coordinates from geographic (lon, lat) to cartesian (x, y) + inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const + { + lp_lat = auth_lat(this->params(), m_proj_parm, lp_lat, 0); + return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); + } + + // INVERSE(e_healpix_inverse) ellipsoid + // Project coordinates from cartesian (x, y) to geographic (lon, lat) + inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const + { + /* Check whether (x, y) lies in the HEALPix image. */ + if (in_image(xy_x, xy_y, 0, 0, 0) == 0) { + lp_lon = HUGE_VAL; + lp_lat = HUGE_VAL; + BOOST_THROW_EXCEPTION( projection_exception(-15) ); + } + healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); + lp_lat = auth_lat(this->params(), m_proj_parm, lp_lat, 1); + } + + static inline std::string get_name() + { + return "healpix_ellipsoid"; + } + + }; + + // template class, using CRTP to implement forward/inverse + template <typename CalculationType, typename Parameters> + struct base_healpix_spheroid : public base_t_fi<base_healpix_spheroid<CalculationType, Parameters>, + CalculationType, Parameters> + { + + typedef CalculationType geographic_type; + typedef CalculationType cartesian_type; + + par_healpix<CalculationType> m_proj_parm; + + inline base_healpix_spheroid(const Parameters& par) + : base_t_fi<base_healpix_spheroid<CalculationType, Parameters>, + CalculationType, Parameters>(*this, par) {} + + // FORWARD(s_healpix_forward) sphere + // Project coordinates from geographic (lon, lat) to cartesian (x, y) + inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const + { + return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); + } + + // INVERSE(s_healpix_inverse) sphere + // Project coordinates from cartesian (x, y) to geographic (lon, lat) + inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const + { + /* Check whether (x, y) lies in the HEALPix image */ + if (in_image(xy_x, xy_y, 0, 0, 0) == 0) { + lp_lon = HUGE_VAL; + lp_lat = HUGE_VAL; + BOOST_THROW_EXCEPTION( projection_exception(-15) ); + } + return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); + } + + static inline std::string get_name() + { + return "healpix_spheroid"; + } + + }; + + // template class, using CRTP to implement forward/inverse + template <typename CalculationType, typename Parameters> + struct base_rhealpix_ellipsoid : public base_t_fi<base_rhealpix_ellipsoid<CalculationType, Parameters>, + CalculationType, Parameters> + { + + typedef CalculationType geographic_type; + typedef CalculationType cartesian_type; + + par_healpix<CalculationType> m_proj_parm; + + inline base_rhealpix_ellipsoid(const Parameters& par) + : base_t_fi<base_rhealpix_ellipsoid<CalculationType, Parameters>, + CalculationType, Parameters>(*this, par) {} + + // FORWARD(e_rhealpix_forward) ellipsoid + // Project coordinates from geographic (lon, lat) to cartesian (x, y) + inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const + { + lp_lat = auth_lat(this->params(), m_proj_parm, lp_lat, 0); + healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); + combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0); + } + + // INVERSE(e_rhealpix_inverse) ellipsoid + // Project coordinates from cartesian (x, y) to geographic (lon, lat) + inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const + { + /* Check whether (x, y) lies in the rHEALPix image. */ + if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) { + lp_lon = HUGE_VAL; + lp_lat = HUGE_VAL; + BOOST_THROW_EXCEPTION( projection_exception(-15) ); + } + combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1); + healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); + lp_lat = auth_lat(this->params(), m_proj_parm, lp_lat, 1); + } + + static inline std::string get_name() + { + return "rhealpix_ellipsoid"; + } + + }; + + // template class, using CRTP to implement forward/inverse + template <typename CalculationType, typename Parameters> + struct base_rhealpix_spheroid : public base_t_fi<base_rhealpix_spheroid<CalculationType, Parameters>, + CalculationType, Parameters> + { + + typedef CalculationType geographic_type; + typedef CalculationType cartesian_type; + + par_healpix<CalculationType> m_proj_parm; + + inline base_rhealpix_spheroid(const Parameters& par) + : base_t_fi<base_rhealpix_spheroid<CalculationType, Parameters>, + CalculationType, Parameters>(*this, par) {} + + // FORWARD(s_rhealpix_forward) sphere + // Project coordinates from geographic (lon, lat) to cartesian (x, y) + inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const + { + healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); + combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0); + } + + // INVERSE(s_rhealpix_inverse) sphere + // Project coordinates from cartesian (x, y) to geographic (lon, lat) + inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const + { + /* Check whether (x, y) lies in the rHEALPix image. */ + if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) { + lp_lon = HUGE_VAL; + lp_lat = HUGE_VAL; + BOOST_THROW_EXCEPTION( projection_exception(-15) ); + } + combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1); + return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); + } + + static inline std::string get_name() + { + return "rhealpix_spheroid"; + } + + }; + + // HEALPix + template <typename Parameters, typename T> + inline void setup_healpix(Parameters& par, par_healpix<T>& proj_parm) + { + if (par.es) { + pj_authset(par.es, proj_parm.apa); /* For auth_lat(). */ + proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */ + par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */ + par.ra = 1.0/par.a; + } else { + } + } + + // rHEALPix + template <typename Parameters, typename T> + inline void setup_rhealpix(Parameters& par, par_healpix<T>& proj_parm) + { + proj_parm.north_square = pj_param(par.params,"inorth_square").i; + proj_parm.south_square = pj_param(par.params,"isouth_square").i; + /* Check for valid north_square and south_square inputs. */ + if (proj_parm.north_square < 0 || proj_parm.north_square > 3) { + BOOST_THROW_EXCEPTION( projection_exception(-47) ); + } + if (proj_parm.south_square < 0 || proj_parm.south_square > 3) { + BOOST_THROW_EXCEPTION( projection_exception(-47) ); + } + if (par.es) { + pj_authset(par.es, proj_parm.apa); /* For auth_lat(). */ + proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */ + par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */ + par.ra = 1.0/par.a; + } else { + } + } + + }} // namespace detail::healpix + #endif // doxygen + + /*! + \brief HEALPix projection + \ingroup projections + \tparam Geographic latlong point type + \tparam Cartesian xy point type + \tparam Parameters parameter type + \par Projection characteristics + - Spheroid + - Ellipsoid + \par Example + \image html ex_healpix.gif + */ + template <typename CalculationType, typename Parameters> + struct healpix_ellipsoid : public detail::healpix::base_healpix_ellipsoid<CalculationType, Parameters> + { + inline healpix_ellipsoid(const Parameters& par) : detail::healpix::base_healpix_ellipsoid<CalculationType, Parameters>(par) + { + detail::healpix::setup_healpix(this->m_par, this->m_proj_parm); + } + }; + + /*! + \brief HEALPix projection + \ingroup projections + \tparam Geographic latlong point type + \tparam Cartesian xy point type + \tparam Parameters parameter type + \par Projection characteristics + - Spheroid + - Ellipsoid + \par Example + \image html ex_healpix.gif + */ + template <typename CalculationType, typename Parameters> + struct healpix_spheroid : public detail::healpix::base_healpix_spheroid<CalculationType, Parameters> + { + inline healpix_spheroid(const Parameters& par) : detail::healpix::base_healpix_spheroid<CalculationType, Parameters>(par) + { + detail::healpix::setup_healpix(this->m_par, this->m_proj_parm); + } + }; + + /*! + \brief rHEALPix projection + \ingroup projections + \tparam Geographic latlong point type + \tparam Cartesian xy point type + \tparam Parameters parameter type + \par Projection characteristics + - Spheroid + - Ellipsoid + \par Projection parameters + - north_square (integer) + - south_square (integer) + \par Example + \image html ex_rhealpix.gif + */ + template <typename CalculationType, typename Parameters> + struct rhealpix_ellipsoid : public detail::healpix::base_rhealpix_ellipsoid<CalculationType, Parameters> + { + inline rhealpix_ellipsoid(const Parameters& par) : detail::healpix::base_rhealpix_ellipsoid<CalculationType, Parameters>(par) + { + detail::healpix::setup_rhealpix(this->m_par, this->m_proj_parm); + } + }; + + /*! + \brief rHEALPix projection + \ingroup projections + \tparam Geographic latlong point type + \tparam Cartesian xy point type + \tparam Parameters parameter type + \par Projection characteristics + - Spheroid + - Ellipsoid + \par Projection parameters + - north_square (integer) + - south_square (integer) + \par Example + \image html ex_rhealpix.gif + */ + template <typename CalculationType, typename Parameters> + struct rhealpix_spheroid : public detail::healpix::base_rhealpix_spheroid<CalculationType, Parameters> + { + inline rhealpix_spheroid(const Parameters& par) : detail::healpix::base_rhealpix_spheroid<CalculationType, Parameters>(par) + { + detail::healpix::setup_rhealpix(this->m_par, this->m_proj_parm); + } + }; + + #ifndef DOXYGEN_NO_DETAIL + namespace detail + { + + // Static projection + BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::healpix, healpix_spheroid, healpix_ellipsoid) + BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::rhealpix, rhealpix_spheroid, rhealpix_ellipsoid) + + // Factory entry(s) + template <typename CalculationType, typename Parameters> + class healpix_entry : public detail::factory_entry<CalculationType, Parameters> + { + public : + virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const + { + if (par.es) + return new base_v_fi<healpix_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par); + else + return new base_v_fi<healpix_spheroid<CalculationType, Parameters>, CalculationType, Parameters>(par); + } + }; + + template <typename CalculationType, typename Parameters> + class rhealpix_entry : public detail::factory_entry<CalculationType, Parameters> + { + public : + virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const + { + if (par.es) + return new base_v_fi<rhealpix_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par); + else + return new base_v_fi<rhealpix_spheroid<CalculationType, Parameters>, CalculationType, Parameters>(par); + } + }; + + template <typename CalculationType, typename Parameters> + inline void healpix_init(detail::base_factory<CalculationType, Parameters>& factory) + { + factory.add_to_factory("healpix", new healpix_entry<CalculationType, Parameters>); + factory.add_to_factory("rhealpix", new rhealpix_entry<CalculationType, Parameters>); + } + + } // namespace detail + #endif // doxygen + +} // namespace projections + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP + |