1 // Boost.Geometry - gis-projections (based on PROJ4) 2 3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. 4 5 // This file was modified by Oracle on 2017, 2018, 2019. 6 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates. 7 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle. 8 9 // Use, modification and distribution is subject to the Boost Software License, 10 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at 11 // http://www.boost.org/LICENSE_1_0.txt) 12 13 // This file is converted from PROJ4, http://trac.osgeo.org/proj 14 // PROJ4 is originally written by Gerald Evenden (then of the USGS) 15 // PROJ4 is maintained by Frank Warmerdam 16 // PROJ4 is converted to Boost.Geometry by Barend Gehrels 17 18 // Last updated version of proj: 5.0.0 19 20 // Original copyright notice: 21 22 // Purpose: Implementation of the HEALPix and rHEALPix projections. 23 // For background see <http://code.scenzgrid.org/index.php/p/scenzgrid-py/source/tree/master/docs/rhealpix_dggs.pdf>. 24 // Authors: Alex Raichev (raichev@cs.auckland.ac.nz) 25 // Michael Speth (spethm@landcareresearch.co.nz) 26 // Notes: Raichev implemented these projections in Python and 27 // Speth translated them into C here. 28 29 // Copyright (c) 2001, Thomas Flemming, tf@ttqv.com 30 31 // Permission is hereby granted, free of charge, to any person obtaining a 32 // copy of this software and associated documentation files (the "Software"), 33 // to deal in the Software without restriction, including without limitation 34 // the rights to use, copy, modify, merge, publish, distribute, sublicense, 35 // and/or sell copies of the Software, and to permit persons to whom the 36 // Software is furnished to do so, subject to the following conditions: 37 38 // The above copyright notice and this permission notice shall be included 39 // in all copies or substantial portions of the Software. 40 41 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 42 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 43 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 44 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 45 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 46 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 47 // DEALINGS IN THE SOFTWARE. 48 49 #ifndef BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP 50 #define BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP 51 52 #include <boost/geometry/srs/projections/impl/base_static.hpp> 53 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> 54 #include <boost/geometry/srs/projections/impl/factory_entry.hpp> 55 #include <boost/geometry/srs/projections/impl/pj_auth.hpp> 56 #include <boost/geometry/srs/projections/impl/pj_param.hpp> 57 #include <boost/geometry/srs/projections/impl/pj_qsfn.hpp> 58 #include <boost/geometry/srs/projections/impl/projects.hpp> 59 60 #include <boost/geometry/util/math.hpp> 61 62 namespace boost { namespace geometry 63 { 64 65 namespace projections 66 { 67 #ifndef DOXYGEN_NO_DETAIL 68 namespace detail { namespace healpix 69 { 70 71 /* Fuzz to handle rounding errors: */ 72 static const double epsilon = 1e-15; 73 74 template <typename T> 75 struct par_healpix 76 { 77 T qp; 78 detail::apa<T> apa; 79 int north_square; 80 int south_square; 81 }; 82 83 template <typename T> 84 struct cap_map 85 { 86 T x, y; /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */ 87 int cn; /* An integer 0--3 indicating the position of the polar cap. */ 88 enum region_type {north, south, equatorial} region; 89 }; 90 template <typename T> 91 struct point_xy 92 { 93 T x, y; 94 }; 95 96 /* IDENT, R1, R2, R3, R1 inverse, R2 inverse, R3 inverse:*/ 97 static double rot[7][2][2] = { 98 /* Identity matrix */ 99 {{1, 0},{0, 1}}, 100 /* Matrix for counterclockwise rotation by pi/2: */ 101 {{ 0,-1},{ 1, 0}}, 102 /* Matrix for counterclockwise rotation by pi: */ 103 {{-1, 0},{ 0,-1}}, 104 /* Matrix for counterclockwise rotation by 3*pi/2: */ 105 {{ 0, 1},{-1, 0}}, 106 {{ 0, 1},{-1, 0}}, // 3*pi/2 107 {{-1, 0},{ 0,-1}}, // pi 108 {{ 0,-1},{ 1, 0}} // pi/2 109 }; 110 111 /** 112 * Returns the sign of the double. 113 * @param v the parameter whose sign is returned. 114 * @return 1 for positive number, -1 for negative, and 0 for zero. 115 **/ 116 template <typename T> pj_sign(T const & v)117 inline T pj_sign (T const& v) 118 { 119 return v > 0 ? 1 : (v < 0 ? -1 : 0); 120 } 121 /** 122 * 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}}}. 123 * @param index ranges from -3 to 3. 124 */ get_rotate_index(int index)125 inline int get_rotate_index(int index) 126 { 127 switch(index) { 128 case 0: 129 return 0; 130 case 1: 131 return 1; 132 case 2: 133 return 2; 134 case 3: 135 return 3; 136 case -1: 137 return 4; 138 case -2: 139 return 5; 140 case -3: 141 return 6; 142 } 143 return 0; 144 } 145 /** 146 * Return 1 if point (testx, testy) lies in the interior of the polygon 147 * determined by the vertices in vert, and return 0 otherwise. 148 * See http://paulbourke.net/geometry/polygonmesh/ for more details. 149 * @param nvert the number of vertices in the polygon. 150 * @param vert the (x, y)-coordinates of the polygon's vertices 151 **/ 152 template <typename T> pnpoly(int nvert,T vert[][2],T const & testx,T const & testy)153 inline int pnpoly(int nvert, T vert[][2], T const& testx, T const& testy) 154 { 155 int i; 156 int counter = 0; 157 T xinters; 158 point_xy<T> p1, p2; 159 160 /* Check for boundrary cases */ 161 for (i = 0; i < nvert; i++) { 162 if (testx == vert[i][0] && testy == vert[i][1]) { 163 return 1; 164 } 165 } 166 167 p1.x = vert[0][0]; 168 p1.y = vert[0][1]; 169 170 for (i = 1; i < nvert; i++) { 171 p2.x = vert[i % nvert][0]; 172 p2.y = vert[i % nvert][1]; 173 if (testy > (std::min)(p1.y, p2.y) && 174 testy <= (std::max)(p1.y, p2.y) && 175 testx <= (std::max)(p1.x, p2.x) && 176 p1.y != p2.y) 177 { 178 xinters = (testy-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x; 179 if (p1.x == p2.x || testx <= xinters) 180 counter++; 181 } 182 p1 = p2; 183 } 184 185 if (counter % 2 == 0) { 186 return 0; 187 } else { 188 return 1; 189 } 190 } 191 /** 192 * Return 1 if (x, y) lies in (the interior or boundary of) the image of the 193 * HEALPix projection (in case proj=0) or in the image the rHEALPix projection 194 * (in case proj=1), and return 0 otherwise. 195 * @param north_square the position of the north polar square (rHEALPix only) 196 * @param south_square the position of the south polar square (rHEALPix only) 197 **/ 198 template <typename T> in_image(T const & x,T const & y,int proj,int north_square,int south_square)199 inline int in_image(T const& x, T const& y, int proj, int north_square, int south_square) 200 { 201 static const T pi = detail::pi<T>(); 202 static const T half_pi = detail::half_pi<T>(); 203 static const T fourth_pi = detail::fourth_pi<T>(); 204 205 if (proj == 0) { 206 T healpixVertsJit[][2] = { 207 {-pi - epsilon, fourth_pi}, 208 {-3.0*fourth_pi, half_pi + epsilon}, 209 {-half_pi, fourth_pi + epsilon}, 210 {-fourth_pi, half_pi + epsilon}, 211 {0.0, fourth_pi + epsilon}, 212 {fourth_pi, half_pi + epsilon}, 213 {half_pi, fourth_pi + epsilon}, 214 {3.0*fourth_pi, half_pi + epsilon}, 215 {pi + epsilon, fourth_pi}, 216 {pi + epsilon, -fourth_pi}, 217 {3.0*fourth_pi, -half_pi - epsilon}, 218 {half_pi, -fourth_pi - epsilon}, 219 {fourth_pi, -half_pi - epsilon}, 220 {0.0, -fourth_pi - epsilon}, 221 {-fourth_pi, -half_pi - epsilon}, 222 {-half_pi, -fourth_pi - epsilon}, 223 {-3.0*fourth_pi, -half_pi - epsilon}, 224 {-pi - epsilon, -fourth_pi} 225 }; 226 return pnpoly((int)sizeof(healpixVertsJit)/ 227 sizeof(healpixVertsJit[0]), healpixVertsJit, x, y); 228 } else { 229 T rhealpixVertsJit[][2] = { 230 {-pi - epsilon, fourth_pi + epsilon}, 231 {-pi + north_square*half_pi - epsilon, fourth_pi + epsilon}, 232 {-pi + north_square*half_pi - epsilon, 3.0*fourth_pi + epsilon}, 233 {-pi + (north_square + 1.0)*half_pi + epsilon, 3.0*fourth_pi + epsilon}, 234 {-pi + (north_square + 1.0)*half_pi + epsilon, fourth_pi + epsilon}, 235 {pi + epsilon, fourth_pi + epsilon}, 236 {pi + epsilon, -fourth_pi - epsilon}, 237 {-pi + (south_square + 1.0)*half_pi + epsilon, -fourth_pi - epsilon}, 238 {-pi + (south_square + 1.0)*half_pi + epsilon, -3.0*fourth_pi - epsilon}, 239 {-pi + south_square*half_pi - epsilon, -3.0*fourth_pi - epsilon}, 240 {-pi + south_square*half_pi - epsilon, -fourth_pi - epsilon}, 241 {-pi - epsilon, -fourth_pi - epsilon} 242 }; 243 244 return pnpoly((int)sizeof(rhealpixVertsJit)/ 245 sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y); 246 } 247 } 248 /** 249 * Return the authalic latitude of latitude alpha (if inverse=0) or 250 * return the approximate latitude of authalic latitude alpha (if inverse=1). 251 * P contains the relavent ellipsoid parameters. 252 **/ 253 template <typename Parameters, typename T> auth_lat(const Parameters & par,const par_healpix<T> & proj_parm,T const & alpha,int inverse)254 inline T auth_lat(const Parameters& par, const par_healpix<T>& proj_parm, T const& alpha, int inverse) 255 { 256 if (inverse == 0) { 257 /* Authalic latitude. */ 258 T q = pj_qsfn(sin(alpha), par.e, 1.0 - par.es); 259 T qp = proj_parm.qp; 260 T ratio = q/qp; 261 262 if (math::abs(ratio) > 1) { 263 /* Rounding error. */ 264 ratio = pj_sign(ratio); 265 } 266 267 return asin(ratio); 268 } else { 269 /* Approximation to inverse authalic latitude. */ 270 return pj_authlat(alpha, proj_parm.apa); 271 } 272 } 273 /** 274 * Return the HEALPix projection of the longitude-latitude point lp on 275 * the unit sphere. 276 **/ 277 template <typename T> healpix_sphere(T const & lp_lam,T const & lp_phi,T & xy_x,T & xy_y)278 inline void healpix_sphere(T const& lp_lam, T const& lp_phi, T& xy_x, T& xy_y) 279 { 280 static const T pi = detail::pi<T>(); 281 static const T half_pi = detail::half_pi<T>(); 282 static const T fourth_pi = detail::fourth_pi<T>(); 283 284 T lam = lp_lam; 285 T phi = lp_phi; 286 T phi0 = asin(T(2.0)/T(3.0)); 287 288 /* equatorial region */ 289 if ( fabsl(phi) <= phi0) { 290 xy_x = lam; 291 xy_y = 3.0*pi/8.0*sin(phi); 292 } else { 293 T lamc; 294 T sigma = sqrt(3.0*(1 - math::abs(sin(phi)))); 295 T cn = floor(2*lam / pi + 2); 296 if (cn >= 4) { 297 cn = 3; 298 } 299 lamc = -3*fourth_pi + half_pi*cn; 300 xy_x = lamc + (lam - lamc)*sigma; 301 xy_y = pj_sign(phi)*fourth_pi*(2 - sigma); 302 } 303 return; 304 } 305 /** 306 * Return the inverse of healpix_sphere(). 307 **/ 308 template <typename T> healpix_sphere_inverse(T const & xy_x,T const & xy_y,T & lp_lam,T & lp_phi)309 inline void healpix_sphere_inverse(T const& xy_x, T const& xy_y, T& lp_lam, T& lp_phi) 310 { 311 static const T pi = detail::pi<T>(); 312 static const T half_pi = detail::half_pi<T>(); 313 static const T fourth_pi = detail::fourth_pi<T>(); 314 315 T x = xy_x; 316 T y = xy_y; 317 T y0 = fourth_pi; 318 319 /* Equatorial region. */ 320 if (math::abs(y) <= y0) { 321 lp_lam = x; 322 lp_phi = asin(8.0*y/(3.0*pi)); 323 } else if (fabsl(y) < half_pi) { 324 T cn = floor(2.0*x/pi + 2.0); 325 T xc, tau; 326 if (cn >= 4) { 327 cn = 3; 328 } 329 xc = -3.0*fourth_pi + (half_pi)*cn; 330 tau = 2.0 - 4.0*fabsl(y)/pi; 331 lp_lam = xc + (x - xc)/tau; 332 lp_phi = pj_sign(y)*asin(1.0 - math::pow(tau, 2)/3.0); 333 } else { 334 lp_lam = -1.0*pi; 335 lp_phi = pj_sign(y)*half_pi; 336 } 337 return; 338 } 339 /** 340 * Return the vector sum a + b, where a and b are 2-dimensional vectors. 341 * @param ret holds a + b. 342 **/ 343 template <typename T> vector_add(const T a[2],const T b[2],T ret[2])344 inline void vector_add(const T a[2], const T b[2], T ret[2]) 345 { 346 int i; 347 for(i = 0; i < 2; i++) { 348 ret[i] = a[i] + b[i]; 349 } 350 } 351 /** 352 * Return the vector difference a - b, where a and b are 2-dimensional vectors. 353 * @param ret holds a - b. 354 **/ 355 template <typename T> vector_sub(const T a[2],const T b[2],T ret[2])356 inline void vector_sub(const T a[2], const T b[2], T ret[2]) 357 { 358 int i; 359 for(i = 0; i < 2; i++) { 360 ret[i] = a[i] - b[i]; 361 } 362 } 363 /** 364 * Return the 2 x 1 matrix product a*b, where a is a 2 x 2 matrix and 365 * b is a 2 x 1 matrix. 366 * @param ret holds a*b. 367 **/ 368 template <typename T1, typename T2> dot_product(const T1 a[2][2],const T2 b[2],T2 ret[2])369 inline void dot_product(const T1 a[2][2], const T2 b[2], T2 ret[2]) 370 { 371 int i, j; 372 int length = 2; 373 for(i = 0; i < length; i++) { 374 ret[i] = 0; 375 for(j = 0; j < length; j++) { 376 ret[i] += a[i][j]*b[j]; 377 } 378 } 379 } 380 /** 381 * Return the number of the polar cap, the pole point coordinates, and 382 * the region that (x, y) lies in. 383 * If inverse=0, then assume (x,y) lies in the image of the HEALPix 384 * projection of the unit sphere. 385 * If inverse=1, then assume (x,y) lies in the image of the 386 * (north_square, south_square)-rHEALPix projection of the unit sphere. 387 **/ 388 template <typename T> get_cap(T x,T const & y,int north_square,int south_square,int inverse)389 inline cap_map<T> get_cap(T x, T const& y, int north_square, int south_square, 390 int inverse) 391 { 392 static const T pi = detail::pi<T>(); 393 static const T half_pi = detail::half_pi<T>(); 394 static const T fourth_pi = detail::fourth_pi<T>(); 395 396 cap_map<T> capmap; 397 T c; 398 capmap.x = x; 399 capmap.y = y; 400 if (inverse == 0) { 401 if (y > fourth_pi) { 402 capmap.region = cap_map<T>::north; 403 c = half_pi; 404 } else if (y < -fourth_pi) { 405 capmap.region = cap_map<T>::south; 406 c = -half_pi; 407 } else { 408 capmap.region = cap_map<T>::equatorial; 409 capmap.cn = 0; 410 return capmap; 411 } 412 /* polar region */ 413 if (x < -half_pi) { 414 capmap.cn = 0; 415 capmap.x = (-3.0*fourth_pi); 416 capmap.y = c; 417 } else if (x >= -half_pi && x < 0) { 418 capmap.cn = 1; 419 capmap.x = -fourth_pi; 420 capmap.y = c; 421 } else if (x >= 0 && x < half_pi) { 422 capmap.cn = 2; 423 capmap.x = fourth_pi; 424 capmap.y = c; 425 } else { 426 capmap.cn = 3; 427 capmap.x = 3.0*fourth_pi; 428 capmap.y = c; 429 } 430 } else { 431 if (y > fourth_pi) { 432 capmap.region = cap_map<T>::north; 433 capmap.x = (-3.0*fourth_pi + north_square*half_pi); 434 capmap.y = half_pi; 435 x = x - north_square*half_pi; 436 } else if (y < -fourth_pi) { 437 capmap.region = cap_map<T>::south; 438 capmap.x = (-3.0*fourth_pi + south_square*pi/2); 439 capmap.y = -half_pi; 440 x = x - south_square*half_pi; 441 } else { 442 capmap.region = cap_map<T>::equatorial; 443 capmap.cn = 0; 444 return capmap; 445 } 446 /* Polar Region, find the HEALPix polar cap number that 447 x, y moves to when rHEALPix polar square is disassembled. */ 448 if (capmap.region == cap_map<T>::north) { 449 if (y >= -x - fourth_pi - epsilon && y < x + 5.0*fourth_pi - epsilon) { 450 capmap.cn = (north_square + 1) % 4; 451 } else if (y > -x -fourth_pi + epsilon && y >= x + 5.0*fourth_pi - epsilon) { 452 capmap.cn = (north_square + 2) % 4; 453 } else if (y <= -x -fourth_pi + epsilon && y > x + 5.0*fourth_pi + epsilon) { 454 capmap.cn = (north_square + 3) % 4; 455 } else { 456 capmap.cn = north_square; 457 } 458 } else if (capmap.region == cap_map<T>::south) { 459 if (y <= x + fourth_pi + epsilon && y > -x - 5.0*fourth_pi + epsilon) { 460 capmap.cn = (south_square + 1) % 4; 461 } else if (y < x + fourth_pi - epsilon && y <= -x - 5.0*fourth_pi + epsilon) { 462 capmap.cn = (south_square + 2) % 4; 463 } else if (y >= x + fourth_pi - epsilon && y < -x - 5.0*fourth_pi - epsilon) { 464 capmap.cn = (south_square + 3) % 4; 465 } else { 466 capmap.cn = south_square; 467 } 468 } 469 } 470 return capmap; 471 } 472 /** 473 * Rearrange point (x, y) in the HEALPix projection by 474 * combining the polar caps into two polar squares. 475 * Put the north polar square in position north_square and 476 * the south polar square in position south_square. 477 * If inverse=1, then uncombine the polar caps. 478 * @param north_square integer between 0 and 3. 479 * @param south_square integer between 0 and 3. 480 **/ 481 template <typename T> combine_caps(T & xy_x,T & xy_y,int north_square,int south_square,int inverse)482 inline void combine_caps(T& xy_x, T& xy_y, int north_square, int south_square, 483 int inverse) 484 { 485 static const T half_pi = detail::half_pi<T>(); 486 static const T fourth_pi = detail::fourth_pi<T>(); 487 488 T v[2]; 489 T c[2]; 490 T vector[2]; 491 T v_min_c[2]; 492 T ret_dot[2]; 493 const double (*tmpRot)[2]; 494 int pole = 0; 495 496 cap_map<T> capmap = get_cap(xy_x, xy_y, north_square, south_square, inverse); 497 if (capmap.region == cap_map<T>::equatorial) { 498 xy_x = capmap.x; 499 xy_y = capmap.y; 500 return; 501 } 502 503 v[0] = xy_x; v[1] = xy_y; 504 c[0] = capmap.x; c[1] = capmap.y; 505 506 if (inverse == 0) { 507 /* Rotate (xy_x, xy_y) about its polar cap tip and then translate it to 508 north_square or south_square. */ 509 510 if (capmap.region == cap_map<T>::north) { 511 pole = north_square; 512 tmpRot = rot[get_rotate_index(capmap.cn - pole)]; 513 } else { 514 pole = south_square; 515 tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))]; 516 } 517 } else { 518 /* Inverse function. 519 Unrotate (xy_x, xy_y) and then translate it back. */ 520 521 /* disassemble */ 522 if (capmap.region == cap_map<T>::north) { 523 pole = north_square; 524 tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))]; 525 } else { 526 pole = south_square; 527 tmpRot = rot[get_rotate_index(capmap.cn - pole)]; 528 } 529 } 530 531 vector_sub(v, c, v_min_c); 532 dot_product(tmpRot, v_min_c, ret_dot); 533 534 { 535 T a[2]; 536 /* Workaround cppcheck git issue */ 537 T* pa = a; 538 // TODO: in proj4 5.0.0 this line is used instead 539 //pa[0] = -3.0*fourth_pi + ((inverse == 0) ? 0 : capmap.cn) *half_pi; 540 pa[0] = -3.0*fourth_pi + ((inverse == 0) ? pole : capmap.cn) *half_pi; 541 pa[1] = half_pi; 542 vector_add(ret_dot, a, vector); 543 } 544 545 xy_x = vector[0]; 546 xy_y = vector[1]; 547 } 548 549 template <typename T, typename Parameters> 550 struct base_healpix_ellipsoid 551 { 552 par_healpix<T> m_proj_parm; 553 554 // FORWARD(e_healpix_forward) ellipsoid 555 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::healpix::base_healpix_ellipsoid556 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const 557 { 558 lp_lat = auth_lat(par, m_proj_parm, lp_lat, 0); 559 return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); 560 } 561 562 // INVERSE(e_healpix_inverse) ellipsoid 563 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::healpix::base_healpix_ellipsoid564 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const 565 { 566 /* Check whether (x, y) lies in the HEALPix image. */ 567 if (in_image(xy_x, xy_y, 0, 0, 0) == 0) { 568 lp_lon = HUGE_VAL; 569 lp_lat = HUGE_VAL; 570 BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); 571 } 572 healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); 573 lp_lat = auth_lat(par, m_proj_parm, lp_lat, 1); 574 } 575 get_nameboost::geometry::projections::detail::healpix::base_healpix_ellipsoid576 static inline std::string get_name() 577 { 578 return "healpix_ellipsoid"; 579 } 580 581 }; 582 583 template <typename T, typename Parameters> 584 struct base_healpix_spheroid 585 { 586 par_healpix<T> m_proj_parm; 587 588 // FORWARD(s_healpix_forward) sphere 589 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::healpix::base_healpix_spheroid590 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 591 { 592 return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); 593 } 594 595 // INVERSE(s_healpix_inverse) sphere 596 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::healpix::base_healpix_spheroid597 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const 598 { 599 /* Check whether (x, y) lies in the HEALPix image */ 600 if (in_image(xy_x, xy_y, 0, 0, 0) == 0) { 601 lp_lon = HUGE_VAL; 602 lp_lat = HUGE_VAL; 603 BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); 604 } 605 return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); 606 } 607 get_nameboost::geometry::projections::detail::healpix::base_healpix_spheroid608 static inline std::string get_name() 609 { 610 return "healpix_spheroid"; 611 } 612 613 }; 614 615 template <typename T, typename Parameters> 616 struct base_rhealpix_ellipsoid 617 { 618 par_healpix<T> m_proj_parm; 619 620 // FORWARD(e_rhealpix_forward) ellipsoid 621 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::healpix::base_rhealpix_ellipsoid622 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const 623 { 624 lp_lat = auth_lat(par, m_proj_parm, lp_lat, 0); 625 healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); 626 combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0); 627 } 628 629 // INVERSE(e_rhealpix_inverse) ellipsoid 630 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::healpix::base_rhealpix_ellipsoid631 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const 632 { 633 /* Check whether (x, y) lies in the rHEALPix image. */ 634 if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) { 635 lp_lon = HUGE_VAL; 636 lp_lat = HUGE_VAL; 637 BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); 638 } 639 combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1); 640 healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); 641 lp_lat = auth_lat(par, m_proj_parm, lp_lat, 1); 642 } 643 get_nameboost::geometry::projections::detail::healpix::base_rhealpix_ellipsoid644 static inline std::string get_name() 645 { 646 return "rhealpix_ellipsoid"; 647 } 648 649 }; 650 651 template <typename T, typename Parameters> 652 struct base_rhealpix_spheroid 653 { 654 par_healpix<T> m_proj_parm; 655 656 // FORWARD(s_rhealpix_forward) sphere 657 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::healpix::base_rhealpix_spheroid658 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 659 { 660 healpix_sphere(lp_lon, lp_lat, xy_x, xy_y); 661 combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0); 662 } 663 664 // INVERSE(s_rhealpix_inverse) sphere 665 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::healpix::base_rhealpix_spheroid666 inline void inv(Parameters const& , T xy_x, T xy_y, T& lp_lon, T& lp_lat) const 667 { 668 /* Check whether (x, y) lies in the rHEALPix image. */ 669 if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) { 670 lp_lon = HUGE_VAL; 671 lp_lat = HUGE_VAL; 672 BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) ); 673 } 674 combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1); 675 return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat); 676 } 677 get_nameboost::geometry::projections::detail::healpix::base_rhealpix_spheroid678 static inline std::string get_name() 679 { 680 return "rhealpix_spheroid"; 681 } 682 683 }; 684 685 // HEALPix 686 template <typename Parameters, typename T> setup_healpix(Parameters & par,par_healpix<T> & proj_parm)687 inline void setup_healpix(Parameters& par, par_healpix<T>& proj_parm) 688 { 689 if (par.es != 0.0) { 690 proj_parm.apa = pj_authset<T>(par.es); /* For auth_lat(). */ 691 proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */ 692 par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */ 693 pj_calc_ellipsoid_params(par, par.a, par.es); /* Ensure we have a consistent parameter set */ 694 } else { 695 } 696 } 697 698 // rHEALPix 699 template <typename Params, typename Parameters, typename T> setup_rhealpix(Params const & params,Parameters & par,par_healpix<T> & proj_parm)700 inline void setup_rhealpix(Params const& params, Parameters& par, par_healpix<T>& proj_parm) 701 { 702 proj_parm.north_square = pj_get_param_i<srs::spar::north_square>(params, "north_square", srs::dpar::north_square); 703 proj_parm.south_square = pj_get_param_i<srs::spar::south_square>(params, "south_square", srs::dpar::south_square); 704 /* Check for valid north_square and south_square inputs. */ 705 if ((proj_parm.north_square < 0) || (proj_parm.north_square > 3)) { 706 BOOST_THROW_EXCEPTION( projection_exception(error_axis) ); 707 } 708 if ((proj_parm.south_square < 0) || (proj_parm.south_square > 3)) { 709 BOOST_THROW_EXCEPTION( projection_exception(error_axis) ); 710 } 711 if (par.es != 0.0) { 712 proj_parm.apa = pj_authset<T>(par.es); /* For auth_lat(). */ 713 proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */ 714 par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */ 715 // TODO: why not the same as in healpix? 716 //pj_calc_ellipsoid_params(par, par.a, par.es); 717 par.ra = 1.0/par.a; 718 } else { 719 } 720 } 721 722 }} // namespace detail::healpix 723 #endif // doxygen 724 725 /*! 726 \brief HEALPix projection 727 \ingroup projections 728 \tparam Geographic latlong point type 729 \tparam Cartesian xy point type 730 \tparam Parameters parameter type 731 \par Projection characteristics 732 - Spheroid 733 - Ellipsoid 734 \par Example 735 \image html ex_healpix.gif 736 */ 737 template <typename T, typename Parameters> 738 struct healpix_ellipsoid : public detail::healpix::base_healpix_ellipsoid<T, Parameters> 739 { 740 template <typename Params> healpix_ellipsoidboost::geometry::projections::healpix_ellipsoid741 inline healpix_ellipsoid(Params const& , Parameters & par) 742 { 743 detail::healpix::setup_healpix(par, this->m_proj_parm); 744 } 745 }; 746 747 /*! 748 \brief HEALPix projection 749 \ingroup projections 750 \tparam Geographic latlong point type 751 \tparam Cartesian xy point type 752 \tparam Parameters parameter type 753 \par Projection characteristics 754 - Spheroid 755 - Ellipsoid 756 \par Example 757 \image html ex_healpix.gif 758 */ 759 template <typename T, typename Parameters> 760 struct healpix_spheroid : public detail::healpix::base_healpix_spheroid<T, Parameters> 761 { 762 template <typename Params> healpix_spheroidboost::geometry::projections::healpix_spheroid763 inline healpix_spheroid(Params const& , Parameters & par) 764 { 765 detail::healpix::setup_healpix(par, this->m_proj_parm); 766 } 767 }; 768 769 /*! 770 \brief rHEALPix projection 771 \ingroup projections 772 \tparam Geographic latlong point type 773 \tparam Cartesian xy point type 774 \tparam Parameters parameter type 775 \par Projection characteristics 776 - Spheroid 777 - Ellipsoid 778 \par Projection parameters 779 - north_square (integer) 780 - south_square (integer) 781 \par Example 782 \image html ex_rhealpix.gif 783 */ 784 template <typename T, typename Parameters> 785 struct rhealpix_ellipsoid : public detail::healpix::base_rhealpix_ellipsoid<T, Parameters> 786 { 787 template <typename Params> rhealpix_ellipsoidboost::geometry::projections::rhealpix_ellipsoid788 inline rhealpix_ellipsoid(Params const& params, Parameters & par) 789 { 790 detail::healpix::setup_rhealpix(params, par, this->m_proj_parm); 791 } 792 }; 793 794 /*! 795 \brief rHEALPix projection 796 \ingroup projections 797 \tparam Geographic latlong point type 798 \tparam Cartesian xy point type 799 \tparam Parameters parameter type 800 \par Projection characteristics 801 - Spheroid 802 - Ellipsoid 803 \par Projection parameters 804 - north_square (integer) 805 - south_square (integer) 806 \par Example 807 \image html ex_rhealpix.gif 808 */ 809 template <typename T, typename Parameters> 810 struct rhealpix_spheroid : public detail::healpix::base_rhealpix_spheroid<T, Parameters> 811 { 812 template <typename Params> rhealpix_spheroidboost::geometry::projections::rhealpix_spheroid813 inline rhealpix_spheroid(Params const& params, Parameters & par) 814 { 815 detail::healpix::setup_rhealpix(params, par, this->m_proj_parm); 816 } 817 }; 818 819 #ifndef DOXYGEN_NO_DETAIL 820 namespace detail 821 { 822 823 // Static projection BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_healpix,healpix_spheroid,healpix_ellipsoid)824 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_healpix, healpix_spheroid, healpix_ellipsoid) 825 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_rhealpix, rhealpix_spheroid, rhealpix_ellipsoid) 826 827 // Factory entry(s) 828 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(healpix_entry, healpix_spheroid, healpix_ellipsoid) 829 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(rhealpix_entry, rhealpix_spheroid, rhealpix_ellipsoid) 830 831 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(healpix_init) 832 { 833 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(healpix, healpix_entry) 834 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(rhealpix, rhealpix_entry) 835 } 836 837 } // namespace detail 838 #endif // doxygen 839 840 } // namespace projections 841 842 }} // namespace boost::geometry 843 844 #endif // BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP 845 846