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 // This code was entirely written by Nathan Wagner 23 // and is in the public domain. 24 25 // Permission is hereby granted, free of charge, to any person obtaining a 26 // copy of this software and associated documentation files (the "Software"), 27 // to deal in the Software without restriction, including without limitation 28 // the rights to use, copy, modify, merge, publish, distribute, sublicense, 29 // and/or sell copies of the Software, and to permit persons to whom the 30 // Software is furnished to do so, subject to the following conditions: 31 32 // The above copyright notice and this permission notice shall be included 33 // in all copies or substantial portions of the Software. 34 35 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 36 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 37 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 38 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 39 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 40 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 41 // DEALINGS IN THE SOFTWARE. 42 43 #ifndef BOOST_GEOMETRY_PROJECTIONS_ISEA_HPP 44 #define BOOST_GEOMETRY_PROJECTIONS_ISEA_HPP 45 46 #include <sstream> 47 48 #include <boost/core/ignore_unused.hpp> 49 50 #include <boost/geometry/core/assert.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_param.hpp> 56 #include <boost/geometry/srs/projections/impl/projects.hpp> 57 58 #include <boost/geometry/util/math.hpp> 59 60 namespace boost { namespace geometry 61 { 62 63 namespace projections 64 { 65 #ifndef DOXYGEN_NO_DETAIL 66 namespace detail { namespace isea 67 { 68 static const double epsilon = std::numeric_limits<double>::epsilon(); 69 70 /* sqrt(5)/M_PI */ 71 static const double isea_scale = 0.8301572857837594396028083; 72 /* 26.565051177 degrees */ 73 static const double v_lat = 0.46364760899944494524; 74 /* 52.62263186 */ 75 static const double e_rad = 0.91843818702186776133; 76 /* 10.81231696 */ 77 static const double f_rad = 0.18871053072122403508; 78 /* R tan(g) sin(60) */ 79 static const double table_g = 0.6615845383; 80 /* H = 0.25 R tan g = */ 81 static const double table_h = 0.1909830056; 82 //static const double RPRIME = 0.91038328153090290025; 83 static const double isea_std_lat = 1.01722196792335072101; 84 static const double isea_std_lon = .19634954084936207740; 85 86 template <typename T> deg30_rad()87 inline T deg30_rad() { return T(30) * geometry::math::d2r<T>(); } 88 template <typename T> deg120_rad()89 inline T deg120_rad() { return T(120) * geometry::math::d2r<T>(); } 90 template <typename T> deg72_rad()91 inline T deg72_rad() { return T(72) * geometry::math::d2r<T>(); } 92 template <typename T> deg90_rad()93 inline T deg90_rad() { return geometry::math::half_pi<T>(); } 94 template <typename T> deg144_rad()95 inline T deg144_rad() { return T(144) * geometry::math::d2r<T>(); } 96 template <typename T> deg36_rad()97 inline T deg36_rad() { return T(36) * geometry::math::d2r<T>(); } 98 template <typename T> deg108_rad()99 inline T deg108_rad() { return T(108) * geometry::math::d2r<T>(); } 100 template <typename T> deg180_rad()101 inline T deg180_rad() { return geometry::math::pi<T>(); } 102 downtri(int tri)103 inline bool downtri(int tri) { return (((tri - 1) / 5) % 2 == 1); } 104 105 /* 106 * Proj 4 provides its own entry points into 107 * the code, so none of the library functions 108 * need to be global 109 */ 110 111 struct hex { 112 int iso; 113 int x, y, z; 114 }; 115 116 /* y *must* be positive down as the xy /iso conversion assumes this */ 117 inline hex_xy(struct hex * h)118 int hex_xy(struct hex *h) { 119 if (!h->iso) return 1; 120 if (h->x >= 0) { 121 h->y = -h->y - (h->x+1)/2; 122 } else { 123 /* need to round toward -inf, not toward zero, so x-1 */ 124 h->y = -h->y - h->x/2; 125 } 126 h->iso = 0; 127 128 return 1; 129 } 130 131 inline hex_iso(struct hex * h)132 int hex_iso(struct hex *h) { 133 if (h->iso) return 1; 134 135 if (h->x >= 0) { 136 h->y = (-h->y - (h->x+1)/2); 137 } else { 138 /* need to round toward -inf, not toward zero, so x-1 */ 139 h->y = (-h->y - (h->x)/2); 140 } 141 142 h->z = -h->x - h->y; 143 h->iso = 1; 144 return 1; 145 } 146 147 template <typename T> 148 inline hexbin2(T const & width,T x,T y,int * i,int * j)149 int hexbin2(T const& width, T x, T y, int *i, int *j) 150 { 151 T z, rx, ry, rz; 152 T abs_dx, abs_dy, abs_dz; 153 int ix, iy, iz, s; 154 struct hex h; 155 156 static const T cos_deg30 = cos(deg30_rad<T>()); 157 158 x = x / cos_deg30; /* rotated X coord */ 159 y = y - x / 2.0; /* adjustment for rotated X */ 160 161 /* adjust for actual hexwidth */ 162 x /= width; 163 y /= width; 164 165 z = -x - y; 166 167 rx = floor(x + 0.5); 168 ix = (int)rx; 169 ry = floor(y + 0.5); 170 iy = (int)ry; 171 rz = floor(z + 0.5); 172 iz = (int)rz; 173 174 s = ix + iy + iz; 175 176 if (s) { 177 abs_dx = fabs(rx - x); 178 abs_dy = fabs(ry - y); 179 abs_dz = fabs(rz - z); 180 181 if (abs_dx >= abs_dy && abs_dx >= abs_dz) { 182 ix -= s; 183 } else if (abs_dy >= abs_dx && abs_dy >= abs_dz) { 184 iy -= s; 185 } else { 186 iz -= s; 187 } 188 } 189 h.x = ix; 190 h.y = iy; 191 h.z = iz; 192 h.iso = 1; 193 194 hex_xy(&h); 195 *i = h.x; 196 *j = h.y; 197 return ix * 100 + iy; 198 } 199 200 //enum isea_poly { isea_none = 0, isea_icosahedron = 20 }; 201 //enum isea_topology { isea_hexagon=6, isea_triangle=3, isea_diamond=4 }; 202 enum isea_address_form { 203 /*isea_addr_geo,*/ isea_addr_q2di, isea_addr_seqnum, 204 /*isea_addr_interleave,*/ isea_addr_plane, isea_addr_q2dd, 205 isea_addr_projtri, isea_addr_vertex2dd, isea_addr_hex 206 }; 207 208 template <typename T> 209 struct isea_dgg { 210 T o_lat, o_lon, o_az; /* orientation, radians */ 211 T radius; /* radius of the earth in meters, ignored 1.0 */ 212 unsigned long serial; 213 //int pole; /* true if standard snyder */ 214 int aperture; /* valid values depend on partitioning method */ 215 int resolution; 216 int triangle; /* triangle of last transformed point */ 217 int quad; /* quad of last transformed point */ 218 //isea_poly polyhedron; /* ignored, icosahedron */ 219 //isea_topology topology; /* ignored, hexagon */ 220 isea_address_form output; /* an isea_address_form */ 221 }; 222 223 template <typename T> 224 struct isea_pt { 225 T x, y; 226 }; 227 228 template <typename T> 229 struct isea_geo { 230 T lon, lat; 231 }; 232 233 template <typename T> 234 struct isea_address { 235 T x,y; /* or i,j or lon,lat depending on type */ 236 int type; /* enum isea_address_form */ 237 int number; 238 }; 239 240 /* ENDINC */ 241 242 enum snyder_polyhedron { 243 snyder_poly_hexagon = 0, snyder_poly_pentagon = 1, 244 snyder_poly_tetrahedron = 2, snyder_poly_cube = 3, 245 snyder_poly_octahedron = 4, snyder_poly_dodecahedron = 5, 246 snyder_poly_icosahedron = 6 247 }; 248 249 template <typename T> 250 struct snyder_constants { 251 T g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b; 252 }; 253 254 template <typename T> constants()255 inline const snyder_constants<T> * constants() 256 { 257 /* TODO put these in radians to avoid a later conversion */ 258 static snyder_constants<T> result[] = { 259 {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0}, 260 {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027}, 261 {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 262 {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 263 {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 264 {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 265 {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0} 266 }; 267 return result; 268 } 269 270 template <typename T> vertex()271 inline const isea_geo<T> * vertex() 272 { 273 static isea_geo<T> result[] = { 274 { 0.0, deg90_rad<T>()}, 275 { deg180_rad<T>(), v_lat}, 276 {-deg108_rad<T>(), v_lat}, 277 {-deg36_rad<T>(), v_lat}, 278 { deg36_rad<T>(), v_lat}, 279 { deg108_rad<T>(), v_lat}, 280 {-deg144_rad<T>(), -v_lat}, 281 {-deg72_rad<T>(), -v_lat}, 282 { 0.0, -v_lat}, 283 { deg72_rad<T>(), -v_lat}, 284 { deg144_rad<T>(), -v_lat}, 285 { 0.0, -deg90_rad<T>()} 286 }; 287 return result; 288 } 289 290 /* TODO make an isea_pt array of the vertices as well */ 291 292 static int tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11}; 293 294 /* triangle Centers */ 295 template <typename T> icostriangles()296 inline const isea_geo<T> * icostriangles() 297 { 298 static isea_geo<T> result[] = { 299 { 0.0, 0.0}, 300 {-deg144_rad<T>(), e_rad}, 301 {-deg72_rad<T>(), e_rad}, 302 { 0.0, e_rad}, 303 { deg72_rad<T>(), e_rad}, 304 { deg144_rad<T>(), e_rad}, 305 {-deg144_rad<T>(), f_rad}, 306 {-deg72_rad<T>(), f_rad}, 307 { 0.0, f_rad}, 308 { deg72_rad<T>(), f_rad}, 309 { deg144_rad<T>(), f_rad}, 310 {-deg108_rad<T>(), -f_rad}, 311 {-deg36_rad<T>(), -f_rad}, 312 { deg36_rad<T>(), -f_rad}, 313 { deg108_rad<T>(), -f_rad}, 314 { deg180_rad<T>(), -f_rad}, 315 {-deg108_rad<T>(), -e_rad}, 316 {-deg36_rad<T>(), -e_rad}, 317 { deg36_rad<T>(), -e_rad}, 318 { deg108_rad<T>(), -e_rad}, 319 { deg180_rad<T>(), -e_rad}, 320 }; 321 return result; 322 } 323 324 template <typename T> az_adjustment(int triangle)325 inline T az_adjustment(int triangle) 326 { 327 T adj; 328 329 isea_geo<T> v; 330 isea_geo<T> c; 331 332 v = vertex<T>()[tri_v1[triangle]]; 333 c = icostriangles<T>()[triangle]; 334 335 /* TODO looks like the adjustment is always either 0 or 180 */ 336 /* at least if you pick your vertex carefully */ 337 adj = atan2(cos(v.lat) * sin(v.lon - c.lon), 338 cos(c.lat) * sin(v.lat) 339 - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon)); 340 return adj; 341 } 342 343 template <typename T> isea_triangle_xy(int triangle)344 inline isea_pt<T> isea_triangle_xy(int triangle) 345 { 346 isea_pt<T> c; 347 T Rprime = 0.91038328153090290025; 348 349 triangle = (triangle - 1) % 20; 350 351 c.x = table_g * ((triangle % 5) - 2) * 2.0; 352 if (triangle > 9) { 353 c.x += table_g; 354 } 355 switch (triangle / 5) { 356 case 0: 357 c.y = 5.0 * table_h; 358 break; 359 case 1: 360 c.y = table_h; 361 break; 362 case 2: 363 c.y = -table_h; 364 break; 365 case 3: 366 c.y = -5.0 * table_h; 367 break; 368 default: 369 /* should be impossible */ 370 BOOST_THROW_EXCEPTION( projection_exception() ); 371 }; 372 c.x *= Rprime; 373 c.y *= Rprime; 374 375 return c; 376 } 377 378 /* snyder eq 14 */ 379 template <typename T> sph_azimuth(T const & f_lon,T const & f_lat,T const & t_lon,T const & t_lat)380 inline T sph_azimuth(T const& f_lon, T const& f_lat, T const& t_lon, T const& t_lat) 381 { 382 T az; 383 384 az = atan2(cos(t_lat) * sin(t_lon - f_lon), 385 cos(f_lat) * sin(t_lat) 386 - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon) 387 ); 388 return az; 389 } 390 391 /* coord needs to be in radians */ 392 template <typename T> isea_snyder_forward(isea_geo<T> * ll,isea_pt<T> * out)393 inline int isea_snyder_forward(isea_geo<T> * ll, isea_pt<T> * out) 394 { 395 static T const two_pi = detail::two_pi<T>(); 396 static T const d2r = geometry::math::d2r<T>(); 397 398 int i; 399 400 /* 401 * spherical distance from center of polygon face to any of its 402 * vertexes on the globe 403 */ 404 T g; 405 406 /* 407 * spherical angle between radius vector to center and adjacent edge 408 * of spherical polygon on the globe 409 */ 410 T G; 411 412 /* 413 * plane angle between radius vector to center and adjacent edge of 414 * plane polygon 415 */ 416 T theta; 417 418 /* additional variables from snyder */ 419 T q, Rprime, H, Ag, Azprime, Az, dprime, f, rho, 420 x, y; 421 422 /* variables used to store intermediate results */ 423 T cot_theta, tan_g, az_offset; 424 425 /* how many multiples of 60 degrees we adjust the azimuth */ 426 int Az_adjust_multiples; 427 428 snyder_constants<T> c; 429 430 /* 431 * TODO by locality of reference, start by trying the same triangle 432 * as last time 433 */ 434 435 /* TODO put these constants in as radians to begin with */ 436 c = constants<T>()[snyder_poly_icosahedron]; 437 theta = c.theta * d2r; 438 g = c.g * d2r; 439 G = c.G * d2r; 440 441 for (i = 1; i <= 20; i++) { 442 T z; 443 isea_geo<T> center; 444 445 center = icostriangles<T>()[i]; 446 447 /* step 1 */ 448 z = acos(sin(center.lat) * sin(ll->lat) 449 + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon)); 450 451 /* not on this triangle */ 452 if (z > g + 0.000005) { /* TODO DBL_EPSILON */ 453 continue; 454 } 455 456 Az = sph_azimuth(center.lon, center.lat, ll->lon, ll->lat); 457 458 /* step 2 */ 459 460 /* This calculates "some" vertex coordinate */ 461 az_offset = az_adjustment<T>(i); 462 463 Az -= az_offset; 464 465 /* TODO I don't know why we do this. It's not in snyder */ 466 /* maybe because we should have picked a better vertex */ 467 if (Az < 0.0) { 468 Az += two_pi; 469 } 470 /* 471 * adjust Az for the point to fall within the range of 0 to 472 * 2(90 - theta) or 60 degrees for the hexagon, by 473 * and therefore 120 degrees for the triangle 474 * of the icosahedron 475 * subtracting or adding multiples of 60 degrees to Az and 476 * recording the amount of adjustment 477 */ 478 479 Az_adjust_multiples = 0; 480 while (Az < 0.0) { 481 Az += deg120_rad<T>(); 482 Az_adjust_multiples--; 483 } 484 while (Az > deg120_rad<T>() + epsilon) { 485 Az -= deg120_rad<T>(); 486 Az_adjust_multiples++; 487 } 488 489 /* step 3 */ 490 cot_theta = 1.0 / tan(theta); 491 tan_g = tan(g); /* TODO this is a constant */ 492 493 /* Calculate q from eq 9. */ 494 /* TODO cot_theta is cot(30) */ 495 q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta); 496 497 /* not in this triangle */ 498 if (z > q + 0.000005) { 499 continue; 500 } 501 /* step 4 */ 502 503 /* Apply equations 5-8 and 10-12 in order */ 504 505 /* eq 5 */ 506 /* Rprime = 0.9449322893 * R; */ 507 /* R' in the paper is for the truncated */ 508 Rprime = 0.91038328153090290025; 509 510 /* eq 6 */ 511 H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G)); 512 513 /* eq 7 */ 514 /* Ag = (Az + G + H - deg180_rad) * M_PI * R * R / deg180_rad; */ 515 Ag = Az + G + H - deg180_rad<T>(); 516 517 /* eq 8 */ 518 Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta); 519 520 /* eq 10 */ 521 /* cot(theta) = 1.73205080756887729355 */ 522 dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta); 523 524 /* eq 11 */ 525 f = dprime / (2.0 * Rprime * sin(q / 2.0)); 526 527 /* eq 12 */ 528 rho = 2.0 * Rprime * f * sin(z / 2.0); 529 530 /* 531 * add back the same 60 degree multiple adjustment from step 532 * 2 to Azprime 533 */ 534 535 Azprime += deg120_rad<T>() * Az_adjust_multiples; 536 537 /* calculate rectangular coordinates */ 538 539 x = rho * sin(Azprime); 540 y = rho * cos(Azprime); 541 542 /* 543 * TODO 544 * translate coordinates to the origin for the particular 545 * hexagon on the flattened polyhedral map plot 546 */ 547 548 out->x = x; 549 out->y = y; 550 551 return i; 552 } 553 554 /* 555 * should be impossible, this implies that the coordinate is not on 556 * any triangle 557 */ 558 559 //fprintf(stderr, "impossible transform: %f %f is not on any triangle\n", 560 // ll->lon * geometry::math::r2d<double>(), ll->lat * geometry::math::r2d<double>()); 561 std::stringstream ss; 562 ss << "impossible transform: " << ll->lon * geometry::math::r2d<T>() 563 << " " << ll->lat * geometry::math::r2d<T>() << " is not on any triangle."; 564 565 BOOST_THROW_EXCEPTION( projection_exception(ss.str()) ); 566 567 /* not reached */ 568 return 0; /* supresses a warning */ 569 } 570 571 /* 572 * return the new coordinates of any point in orginal coordinate system. 573 * Define a point (newNPold) in orginal coordinate system as the North Pole in 574 * new coordinate system, and the great circle connect the original and new 575 * North Pole as the lon0 longitude in new coordinate system, given any point 576 * in orginal coordinate system, this function return the new coordinates. 577 */ 578 579 580 /* formula from Snyder, Map Projections: A working manual, p31 */ 581 /* 582 * old north pole at np in new coordinates 583 * could be simplified a bit with fewer intermediates 584 * 585 * TODO take a result pointer 586 */ 587 template <typename T> snyder_ctran(isea_geo<T> * np,isea_geo<T> * pt)588 inline isea_geo<T> snyder_ctran(isea_geo<T> * np, isea_geo<T> * pt) 589 { 590 static T const pi = detail::pi<T>(); 591 static T const two_pi = detail::two_pi<T>(); 592 593 isea_geo<T> npt; 594 T alpha, phi, lambda, lambda0, beta, lambdap, phip; 595 T sin_phip; 596 T lp_b; /* lambda prime minus beta */ 597 T cos_p, sin_a; 598 599 phi = pt->lat; 600 lambda = pt->lon; 601 alpha = np->lat; 602 beta = np->lon; 603 lambda0 = beta; 604 605 cos_p = cos(phi); 606 sin_a = sin(alpha); 607 608 /* mpawm 5-7 */ 609 sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0); 610 611 /* mpawm 5-8b */ 612 613 /* use the two argument form so we end up in the right quadrant */ 614 lp_b = atan2(cos_p * sin(lambda - lambda0), 615 (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi))); 616 617 lambdap = lp_b + beta; 618 619 /* normalize longitude */ 620 /* TODO can we just do a modulus ? */ 621 lambdap = fmod(lambdap, two_pi); 622 while (lambdap > pi) 623 lambdap -= two_pi; 624 while (lambdap < -pi) 625 lambdap += two_pi; 626 627 phip = asin(sin_phip); 628 629 npt.lat = phip; 630 npt.lon = lambdap; 631 632 return npt; 633 } 634 635 template <typename T> isea_ctran(isea_geo<T> * np,isea_geo<T> * pt,T const & lon0)636 inline isea_geo<T> isea_ctran(isea_geo<T> * np, isea_geo<T> * pt, T const& lon0) 637 { 638 static T const pi = detail::pi<T>(); 639 static T const two_pi = detail::two_pi<T>(); 640 641 isea_geo<T> npt; 642 643 np->lon += pi; 644 npt = snyder_ctran(np, pt); 645 np->lon -= pi; 646 647 npt.lon -= (pi - lon0 + np->lon); 648 649 /* 650 * snyder is down tri 3, isea is along side of tri1 from vertex 0 to 651 * vertex 1 these are 180 degrees apart 652 */ 653 npt.lon += pi; 654 /* normalize longitude */ 655 npt.lon = fmod(npt.lon, two_pi); 656 while (npt.lon > pi) 657 npt.lon -= two_pi; 658 while (npt.lon < -pi) 659 npt.lon += two_pi; 660 661 return npt; 662 } 663 664 /* in radians */ 665 666 /* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */ 667 668 template <typename T> isea_grid_init(isea_dgg<T> * g)669 inline int isea_grid_init(isea_dgg<T> * g) 670 { 671 if (!g) 672 return 0; 673 674 //g->polyhedron = isea_icosahedron; 675 g->o_lat = isea_std_lat; 676 g->o_lon = isea_std_lon; 677 g->o_az = 0.0; 678 g->aperture = 4; 679 g->resolution = 6; 680 g->radius = 1.0; 681 //g->topology = isea_hexagon; 682 683 return 1; 684 } 685 686 template <typename T> isea_orient_isea(isea_dgg<T> * g)687 inline int isea_orient_isea(isea_dgg<T> * g) 688 { 689 if (!g) 690 return 0; 691 g->o_lat = isea_std_lat; 692 g->o_lon = isea_std_lon; 693 g->o_az = 0.0; 694 return 1; 695 } 696 697 template <typename T> isea_orient_pole(isea_dgg<T> * g)698 inline int isea_orient_pole(isea_dgg<T> * g) 699 { 700 static T const half_pi = detail::half_pi<T>(); 701 702 if (!g) 703 return 0; 704 g->o_lat = half_pi; 705 g->o_lon = 0.0; 706 g->o_az = 0; 707 return 1; 708 } 709 710 template <typename T> isea_transform(isea_dgg<T> * g,isea_geo<T> * in,isea_pt<T> * out)711 inline int isea_transform(isea_dgg<T> * g, isea_geo<T> * in, 712 isea_pt<T> * out) 713 { 714 isea_geo<T> i, pole; 715 int tri; 716 717 pole.lat = g->o_lat; 718 pole.lon = g->o_lon; 719 720 i = isea_ctran(&pole, in, g->o_az); 721 722 tri = isea_snyder_forward(&i, out); 723 out->x *= g->radius; 724 out->y *= g->radius; 725 g->triangle = tri; 726 727 return tri; 728 } 729 730 731 template <typename T> isea_rotate(isea_pt<T> * pt,T const & degrees)732 inline void isea_rotate(isea_pt<T> * pt, T const& degrees) 733 { 734 static T const d2r = geometry::math::d2r<T>(); 735 static T const two_pi = detail::two_pi<T>(); 736 737 T rad; 738 739 T x, y; 740 741 rad = -degrees * d2r; 742 while (rad >= two_pi) rad -= two_pi; 743 while (rad <= -two_pi) rad += two_pi; 744 745 x = pt->x * cos(rad) + pt->y * sin(rad); 746 y = -pt->x * sin(rad) + pt->y * cos(rad); 747 748 pt->x = x; 749 pt->y = y; 750 } 751 752 template <typename T> isea_tri_plane(int tri,isea_pt<T> * pt,T const & radius)753 inline int isea_tri_plane(int tri, isea_pt<T> *pt, T const& radius) 754 { 755 isea_pt<T> tc; /* center of triangle */ 756 757 if (downtri(tri)) { 758 isea_rotate(pt, 180.0); 759 } 760 tc = isea_triangle_xy<T>(tri); 761 tc.x *= radius; 762 tc.y *= radius; 763 pt->x += tc.x; 764 pt->y += tc.y; 765 766 return tri; 767 } 768 769 /* convert projected triangle coords to quad xy coords, return quad number */ 770 template <typename T> isea_ptdd(int tri,isea_pt<T> * pt)771 inline int isea_ptdd(int tri, isea_pt<T> *pt) 772 { 773 int downtri, quad; 774 775 downtri = (((tri - 1) / 5) % 2 == 1); 776 quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1; 777 778 isea_rotate(pt, downtri ? 240.0 : 60.0); 779 if (downtri) { 780 pt->x += 0.5; 781 /* pt->y += cos(30.0 * M_PI / 180.0); */ 782 pt->y += .86602540378443864672; 783 } 784 return quad; 785 } 786 787 template <typename T> isea_dddi_ap3odd(isea_dgg<T> * g,int quad,isea_pt<T> * pt,isea_pt<T> * di)788 inline int isea_dddi_ap3odd(isea_dgg<T> *g, int quad, isea_pt<T> *pt, isea_pt<T> *di) 789 { 790 static T const pi = detail::pi<T>(); 791 792 isea_pt<T> v; 793 T hexwidth; 794 T sidelength; /* in hexes */ 795 int d, i; 796 int maxcoord; 797 hex h; 798 799 /* This is the number of hexes from apex to base of a triangle */ 800 sidelength = (math::pow(T(2), g->resolution) + T(1)) / T(2); 801 802 /* apex to base is cos(30deg) */ 803 hexwidth = cos(pi / 6.0) / sidelength; 804 805 /* TODO I think sidelength is always x.5, so 806 * (int)sidelength * 2 + 1 might be just as good 807 */ 808 maxcoord = (int) (sidelength * 2.0 + 0.5); 809 810 v = *pt; 811 hexbin2(hexwidth, v.x, v.y, &h.x, &h.y); 812 h.iso = 0; 813 hex_iso(&h); 814 815 d = h.x - h.z; 816 i = h.x + h.y + h.y; 817 818 /* 819 * you want to test for max coords for the next quad in the same 820 * "row" first to get the case where both are max 821 */ 822 if (quad <= 5) { 823 if (d == 0 && i == maxcoord) { 824 /* north pole */ 825 quad = 0; 826 d = 0; 827 i = 0; 828 } else if (i == maxcoord) { 829 /* upper right in next quad */ 830 quad += 1; 831 if (quad == 6) 832 quad = 1; 833 i = maxcoord - d; 834 d = 0; 835 } else if (d == maxcoord) { 836 /* lower right in quad to lower right */ 837 quad += 5; 838 d = 0; 839 } 840 } else if (quad >= 6) { 841 if (i == 0 && d == maxcoord) { 842 /* south pole */ 843 quad = 11; 844 d = 0; 845 i = 0; 846 } else if (d == maxcoord) { 847 /* lower right in next quad */ 848 quad += 1; 849 if (quad == 11) 850 quad = 6; 851 d = maxcoord - i; 852 i = 0; 853 } else if (i == maxcoord) { 854 /* upper right in quad to upper right */ 855 quad = (quad - 4) % 5; 856 i = 0; 857 } 858 } 859 860 di->x = d; 861 di->y = i; 862 863 g->quad = quad; 864 return quad; 865 } 866 867 template <typename T> isea_dddi(isea_dgg<T> * g,int quad,isea_pt<T> * pt,isea_pt<T> * di)868 inline int isea_dddi(isea_dgg<T> *g, int quad, isea_pt<T> *pt, isea_pt<T> *di) 869 { 870 isea_pt<T> v; 871 T hexwidth; 872 int sidelength; /* in hexes */ 873 hex h; 874 875 if (g->aperture == 3 && g->resolution % 2 != 0) { 876 return isea_dddi_ap3odd(g, quad, pt, di); 877 } 878 /* todo might want to do this as an iterated loop */ 879 if (g->aperture >0) { 880 sidelength = (int) (math::pow(T(g->aperture), T(g->resolution / T(2))) + T(0.5)); 881 } else { 882 sidelength = g->resolution; 883 } 884 885 hexwidth = 1.0 / sidelength; 886 887 v = *pt; 888 isea_rotate(&v, -30.0); 889 hexbin2(hexwidth, v.x, v.y, &h.x, &h.y); 890 h.iso = 0; 891 hex_iso(&h); 892 893 /* we may actually be on another quad */ 894 if (quad <= 5) { 895 if (h.x == 0 && h.z == -sidelength) { 896 /* north pole */ 897 quad = 0; 898 h.z = 0; 899 h.y = 0; 900 h.x = 0; 901 } else if (h.z == -sidelength) { 902 quad = quad + 1; 903 if (quad == 6) 904 quad = 1; 905 h.y = sidelength - h.x; 906 h.z = h.x - sidelength; 907 h.x = 0; 908 } else if (h.x == sidelength) { 909 quad += 5; 910 h.y = -h.z; 911 h.x = 0; 912 } 913 } else if (quad >= 6) { 914 if (h.z == 0 && h.x == sidelength) { 915 /* south pole */ 916 quad = 11; 917 h.x = 0; 918 h.y = 0; 919 h.z = 0; 920 } else if (h.x == sidelength) { 921 quad = quad + 1; 922 if (quad == 11) 923 quad = 6; 924 h.x = h.y + sidelength; 925 h.y = 0; 926 h.z = -h.x; 927 } else if (h.y == -sidelength) { 928 quad -= 4; 929 h.y = 0; 930 h.z = -h.x; 931 } 932 } 933 di->x = h.x; 934 di->y = -h.z; 935 936 g->quad = quad; 937 return quad; 938 } 939 940 template <typename T> isea_ptdi(isea_dgg<T> * g,int tri,isea_pt<T> * pt,isea_pt<T> * di)941 inline int isea_ptdi(isea_dgg<T> *g, int tri, isea_pt<T> *pt, 942 isea_pt<T> *di) 943 { 944 isea_pt<T> v; 945 int quad; 946 947 v = *pt; 948 quad = isea_ptdd(tri, &v); 949 quad = isea_dddi(g, quad, &v, di); 950 return quad; 951 } 952 953 /* q2di to seqnum */ 954 template <typename T> isea_disn(isea_dgg<T> * g,int quad,isea_pt<T> * di)955 inline int isea_disn(isea_dgg<T> *g, int quad, isea_pt<T> *di) 956 { 957 int sidelength; 958 int sn, height; 959 int hexes; 960 961 if (quad == 0) { 962 g->serial = 1; 963 return g->serial; 964 } 965 /* hexes in a quad */ 966 hexes = (int) (math::pow(T(g->aperture), T(g->resolution)) + T(0.5)); 967 if (quad == 11) { 968 g->serial = 1 + 10 * hexes + 1; 969 return g->serial; 970 } 971 if (g->aperture == 3 && g->resolution % 2 == 1) { 972 height = (int) (math::pow(T(g->aperture), T((g->resolution - 1) / T(2)))); 973 sn = ((int) di->x) * height; 974 sn += ((int) di->y) / height; 975 sn += (quad - 1) * hexes; 976 sn += 2; 977 } else { 978 sidelength = (int) (math::pow(T(g->aperture), T(g->resolution / T(2))) + T(0.5)); 979 sn = (int) ((quad - 1) * hexes + sidelength * di->x + di->y + 2); 980 } 981 982 g->serial = sn; 983 return sn; 984 } 985 986 /* TODO just encode the quad in the d or i coordinate 987 * quad is 0-11, which can be four bits. 988 * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf 989 */ 990 /* convert a q2di to global hex coord */ 991 template <typename T> isea_hex(isea_dgg<T> * g,int tri,isea_pt<T> * pt,isea_pt<T> * hex)992 inline int isea_hex(isea_dgg<T> *g, int tri, isea_pt<T> *pt, 993 isea_pt<T> *hex) 994 { 995 isea_pt<T> v; 996 #ifdef BOOST_GEOMETRY_PROJECTIONS_FIXME 997 int sidelength; 998 int d, i, x, y; 999 #endif // BOOST_GEOMETRY_PROJECTIONS_FIXME 1000 int quad; 1001 1002 quad = isea_ptdi(g, tri, pt, &v); 1003 1004 hex->x = ((int)v.x << 4) + quad; 1005 hex->y = v.y; 1006 1007 return 1; 1008 #ifdef BOOST_GEOMETRY_PROJECTIONS_FIXME 1009 d = (int)v.x; 1010 i = (int)v.y; 1011 1012 /* Aperture 3 odd resolutions */ 1013 if (g->aperture == 3 && g->resolution % 2 != 0) { 1014 int offset = (int)(pow(T(3.0), T(g->resolution - 1)) + 0.5); 1015 1016 d += offset * ((g->quad-1) % 5); 1017 i += offset * ((g->quad-1) % 5); 1018 1019 if (quad == 0) { 1020 d = 0; 1021 i = offset; 1022 } else if (quad == 11) { 1023 d = 2 * offset; 1024 i = 0; 1025 } else if (quad > 5) { 1026 d += offset; 1027 } 1028 1029 x = (2*d - i) /3; 1030 y = (2*i - d) /3; 1031 1032 hex->x = x + offset / 3; 1033 hex->y = y + 2 * offset / 3; 1034 return 1; 1035 } 1036 1037 /* aperture 3 even resolutions and aperture 4 */ 1038 sidelength = (int) (pow(T(g->aperture), T(g->resolution / 2.0)) + 0.5); 1039 if (g->quad == 0) { 1040 hex->x = 0; 1041 hex->y = sidelength; 1042 } else if (g->quad == 11) { 1043 hex->x = sidelength * 2; 1044 hex->y = 0; 1045 } else { 1046 hex->x = d + sidelength * ((g->quad-1) % 5); 1047 if (g->quad > 5) hex->x += sidelength; 1048 hex->y = i + sidelength * ((g->quad-1) % 5); 1049 } 1050 1051 return 1; 1052 #endif // BOOST_GEOMETRY_PROJECTIONS_FIXME 1053 } 1054 1055 template <typename T> isea_forward(isea_dgg<T> * g,isea_geo<T> * in)1056 inline isea_pt<T> isea_forward(isea_dgg<T> *g, isea_geo<T> *in) 1057 { 1058 int tri; 1059 isea_pt<T> out, coord; 1060 1061 tri = isea_transform(g, in, &out); 1062 1063 if (g->output == isea_addr_plane) { 1064 isea_tri_plane(tri, &out, g->radius); 1065 return out; 1066 } 1067 1068 /* convert to isea standard triangle size */ 1069 out.x = out.x / g->radius * isea_scale; 1070 out.y = out.y / g->radius * isea_scale; 1071 out.x += 0.5; 1072 out.y += 2.0 * .14433756729740644112; 1073 1074 switch (g->output) { 1075 case isea_addr_projtri: 1076 /* nothing to do, already in projected triangle */ 1077 break; 1078 case isea_addr_vertex2dd: 1079 g->quad = isea_ptdd(tri, &out); 1080 break; 1081 case isea_addr_q2dd: 1082 /* Same as above, we just don't print as much */ 1083 g->quad = isea_ptdd(tri, &out); 1084 break; 1085 case isea_addr_q2di: 1086 g->quad = isea_ptdi(g, tri, &out, &coord); 1087 return coord; 1088 break; 1089 case isea_addr_seqnum: 1090 isea_ptdi(g, tri, &out, &coord); 1091 /* disn will set g->serial */ 1092 isea_disn(g, g->quad, &coord); 1093 return coord; 1094 break; 1095 case isea_addr_hex: 1096 isea_hex(g, tri, &out, &coord); 1097 return coord; 1098 break; 1099 default: 1100 // isea_addr_plane handled above 1101 BOOST_GEOMETRY_ASSERT(false); 1102 break; 1103 } 1104 1105 return out; 1106 } 1107 /* 1108 * Proj 4 integration code follows 1109 */ 1110 1111 template <typename T> 1112 struct par_isea 1113 { 1114 isea_dgg<T> dgg; 1115 }; 1116 1117 template <typename T, typename Parameters> 1118 struct base_isea_spheroid 1119 { 1120 par_isea<T> m_proj_parm; 1121 1122 // FORWARD(s_forward) 1123 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::isea::base_isea_spheroid1124 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 1125 { 1126 isea_pt<T> out; 1127 isea_geo<T> in; 1128 1129 in.lon = lp_lon; 1130 in.lat = lp_lat; 1131 1132 isea_dgg<T> copy = this->m_proj_parm.dgg; 1133 out = isea_forward(©, &in); 1134 1135 xy_x = out.x; 1136 xy_y = out.y; 1137 } 1138 get_nameboost::geometry::projections::detail::isea::base_isea_spheroid1139 static inline std::string get_name() 1140 { 1141 return "isea_spheroid"; 1142 } 1143 1144 }; 1145 1146 template <typename T> isea_orient_init(srs::detail::proj4_parameters const & params,par_isea<T> & proj_parm)1147 inline void isea_orient_init(srs::detail::proj4_parameters const& params, 1148 par_isea<T>& proj_parm) 1149 { 1150 std::string opt = pj_get_param_s(params, "orient"); 1151 if (! opt.empty()) { 1152 if (opt == std::string("isea")) { 1153 isea_orient_isea(&proj_parm.dgg); 1154 } else if (opt == std::string("pole")) { 1155 isea_orient_pole(&proj_parm.dgg); 1156 } else { 1157 BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) ); 1158 } 1159 } 1160 } 1161 1162 template <typename T> isea_orient_init(srs::dpar::parameters<T> const & params,par_isea<T> & proj_parm)1163 inline void isea_orient_init(srs::dpar::parameters<T> const& params, 1164 par_isea<T>& proj_parm) 1165 { 1166 typename srs::dpar::parameters<T>::const_iterator 1167 it = pj_param_find(params, srs::dpar::orient); 1168 if (it != params.end()) { 1169 srs::dpar::value_orient o = static_cast<srs::dpar::value_orient>(it->template get_value<int>()); 1170 if (o == srs::dpar::orient_isea) { 1171 isea_orient_isea(&proj_parm.dgg); 1172 } else if (o == srs::dpar::orient_pole) { 1173 isea_orient_pole(&proj_parm.dgg); 1174 } else { 1175 BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) ); 1176 } 1177 } 1178 } 1179 1180 template <typename T> isea_mode_init(srs::detail::proj4_parameters const & params,par_isea<T> & proj_parm)1181 inline void isea_mode_init(srs::detail::proj4_parameters const& params, 1182 par_isea<T>& proj_parm) 1183 { 1184 std::string opt = pj_get_param_s(params, "mode"); 1185 if (! opt.empty()) { 1186 if (opt == std::string("plane")) { 1187 proj_parm.dgg.output = isea_addr_plane; 1188 } else if (opt == std::string("di")) { 1189 proj_parm.dgg.output = isea_addr_q2di; 1190 } else if (opt == std::string("dd")) { 1191 proj_parm.dgg.output = isea_addr_q2dd; 1192 } else if (opt == std::string("hex")) { 1193 proj_parm.dgg.output = isea_addr_hex; 1194 } else { 1195 BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) ); 1196 } 1197 } 1198 } 1199 1200 template <typename T> isea_mode_init(srs::dpar::parameters<T> const & params,par_isea<T> & proj_parm)1201 inline void isea_mode_init(srs::dpar::parameters<T> const& params, 1202 par_isea<T>& proj_parm) 1203 { 1204 typename srs::dpar::parameters<T>::const_iterator 1205 it = pj_param_find(params, srs::dpar::mode); 1206 if (it != params.end()) { 1207 srs::dpar::value_mode m = static_cast<srs::dpar::value_mode>(it->template get_value<int>()); 1208 if (m == srs::dpar::mode_plane) { 1209 proj_parm.dgg.output = isea_addr_plane; 1210 } else if (m == srs::dpar::mode_di) { 1211 proj_parm.dgg.output = isea_addr_q2di; 1212 } else if (m == srs::dpar::mode_dd) { 1213 proj_parm.dgg.output = isea_addr_q2dd; 1214 } else if (m == srs::dpar::mode_hex) { 1215 proj_parm.dgg.output = isea_addr_hex; 1216 } else { 1217 BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) ); 1218 } 1219 } 1220 } 1221 1222 // Icosahedral Snyder Equal Area 1223 template <typename Params, typename T> setup_isea(Params const & params,par_isea<T> & proj_parm)1224 inline void setup_isea(Params const& params, par_isea<T>& proj_parm) 1225 { 1226 std::string opt; 1227 1228 isea_grid_init(&proj_parm.dgg); 1229 1230 proj_parm.dgg.output = isea_addr_plane; 1231 /* proj_parm.dgg.radius = par.a; / * otherwise defaults to 1 */ 1232 /* calling library will scale, I think */ 1233 1234 isea_orient_init(params, proj_parm); 1235 1236 pj_param_r<srs::spar::azi>(params, "azi", srs::dpar::azi, proj_parm.dgg.o_az); 1237 pj_param_r<srs::spar::lon_0>(params, "lon_0", srs::dpar::lon_0, proj_parm.dgg.o_lon); 1238 pj_param_r<srs::spar::lat_0>(params, "lat_0", srs::dpar::lat_0, proj_parm.dgg.o_lat); 1239 // TODO: this parameter is set below second time 1240 pj_param_i<srs::spar::aperture>(params, "aperture", srs::dpar::aperture, proj_parm.dgg.aperture); 1241 // TODO: this parameter is set below second time 1242 pj_param_i<srs::spar::resolution>(params, "resolution", srs::dpar::resolution, proj_parm.dgg.resolution); 1243 1244 isea_mode_init(params, proj_parm); 1245 1246 // TODO: pj_param_exists -> pj_get_param_b ? 1247 if (pj_param_exists<srs::spar::rescale>(params, "rescale", srs::dpar::rescale)) { 1248 proj_parm.dgg.radius = isea_scale; 1249 } 1250 1251 if (pj_param_i<srs::spar::resolution>(params, "resolution", srs::dpar::resolution, proj_parm.dgg.resolution)) { 1252 /* empty */ 1253 } else { 1254 proj_parm.dgg.resolution = 4; 1255 } 1256 1257 if (pj_param_i<srs::spar::aperture>(params, "aperture", srs::dpar::aperture, proj_parm.dgg.aperture)) { 1258 /* empty */ 1259 } else { 1260 proj_parm.dgg.aperture = 3; 1261 } 1262 } 1263 1264 }} // namespace detail::isea 1265 #endif // doxygen 1266 1267 /*! 1268 \brief Icosahedral Snyder Equal Area projection 1269 \ingroup projections 1270 \tparam Geographic latlong point type 1271 \tparam Cartesian xy point type 1272 \tparam Parameters parameter type 1273 \par Projection characteristics 1274 - Spheroid 1275 \par Projection parameters 1276 - orient (string) 1277 - azi: Azimuth (or Gamma) (degrees) 1278 - lon_0: Central meridian (degrees) 1279 - lat_0: Latitude of origin (degrees) 1280 - aperture (integer) 1281 - resolution (integer) 1282 - mode (string) 1283 - rescale 1284 \par Example 1285 \image html ex_isea.gif 1286 */ 1287 template <typename T, typename Parameters> 1288 struct isea_spheroid : public detail::isea::base_isea_spheroid<T, Parameters> 1289 { 1290 template <typename Params> isea_spheroidboost::geometry::projections::isea_spheroid1291 inline isea_spheroid(Params const& params, Parameters const& ) 1292 { 1293 detail::isea::setup_isea(params, this->m_proj_parm); 1294 } 1295 }; 1296 1297 #ifndef DOXYGEN_NO_DETAIL 1298 namespace detail 1299 { 1300 1301 // Static projection BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_F(srs::spar::proj_isea,isea_spheroid)1302 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_F(srs::spar::proj_isea, isea_spheroid) 1303 1304 // Factory entry(s) 1305 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_F(isea_entry, isea_spheroid) 1306 1307 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(isea_init) 1308 { 1309 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(isea, isea_entry) 1310 } 1311 1312 } // namespace detail 1313 #endif // doxygen 1314 1315 } // namespace projections 1316 1317 }} // namespace boost::geometry 1318 1319 #endif // BOOST_GEOMETRY_PROJECTIONS_ISEA_HPP 1320 1321