5 enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 };
6 enum isea_topology { ISEA_HEXAGON=6, ISEA_TRIANGLE=3, ISEA_DIAMOND=4 };
7 enum isea_address_form { ISEA_GEO, ISEA_Q2DI, ISEA_SEQNUM, ISEA_INTERLEAVE,
8 ISEA_PLANE, ISEA_Q2DD, ISEA_PROJTRI, ISEA_VERTEX2DD, ISEA_HEX
12 int polyhedron; /* ignored, icosahedron */
13 double o_lat, o_lon, o_az; /* orientation, radians */
14 int pole; /* true if standard snyder */
15 int topology; /* ignored, hexagon */
16 int aperture; /* valid values depend on partitioning method */
18 double radius; /* radius of the earth in meters, ignored 1.0 */
19 int output; /* an isea_address_form */
20 int triangle; /* triangle of last transformed point */
21 int quad; /* quad of last transformed point */
34 int type; /* enum isea_address_form */
36 double x,y; /* or i,j or lon,lat depending on type */
41 enum snyder_polyhedron {
42 SNYDER_POLY_HEXAGON, SNYDER_POLY_PENTAGON,
43 SNYDER_POLY_TETRAHEDRON, SNYDER_POLY_CUBE,
44 SNYDER_POLY_OCTAHEDRON, SNYDER_POLY_DODECAHEDRON,
45 SNYDER_POLY_ICOSAHEDRON
48 struct snyder_constants {
49 double g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b;
52 /* TODO put these in radians to avoid a later conversion */
54 struct snyder_constants constants[] = {
55 {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0},
56 {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027},
61 {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0},
67 #define DEG60 1.04719755119659774614
68 #define DEG120 2.09439510239319549229
69 #define DEG72 1.25663706143591729537
70 #define DEG90 1.57079632679489661922
71 #define DEG144 2.51327412287183459075
72 #define DEG36 0.62831853071795864768
73 #define DEG108 1.88495559215387594306
76 #define ISEA_SCALE 0.8301572857837594396028083
78 /* 26.565051177 degrees */
79 #define V_LAT 0.46364760899944494524
81 #define RAD2DEG (180.0/M_PI)
82 #define DEG2RAD (M_PI/180.0)
85 struct isea_geo vertex[] = {
100 /* TODO make an isea_pt array of the vertices as well */
102 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};
105 #define E_RAD 0.91843818702186776133
108 #define F_RAD 0.18871053072122403508
110 /* triangle Centers */
111 struct isea_geo icostriangles[] = {
136 az_adjustment(int triangle)
143 v = vertex[tri_v1[triangle]];
144 c = icostriangles[triangle];
146 /* TODO looks like the adjustment is always either 0 or 180 */
147 /* at least if you pick your vertex carefully */
148 adj = atan2(cos(v.lat) * sin(v.lon - c.lon),
149 cos(c.lat) * sin(v.lat)
150 - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon));
154 /* R tan(g) sin(60) */
155 #define TABLE_G 0.6615845383
157 /* H = 0.25 R tan g = */
158 #define TABLE_H 0.1909830056
160 #define RPRIME 0.91038328153090290025
164 isea_triangle_xy(int triangle)
167 double Rprime = 0.91038328153090290025;
169 triangle = (triangle - 1) % 20;
171 c.x = TABLE_G * ((triangle % 5) - 2) * 2.0;
175 switch (triangle / 5) {
186 c.y = -5.0 * TABLE_H;
189 /* should be impossible */
200 sph_azimuth(double f_lon, double f_lat, double t_lon, double t_lat)
204 az = atan2(cos(t_lat) * sin(t_lon - f_lon),
205 cos(f_lat) * sin(t_lat)
206 - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)
211 /* coord needs to be in radians */
214 isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
219 * spherical distance from center of polygon face to any of its
220 * vertexes on the globe
225 * spherical angle between radius vector to center and adjacent edge
226 * of spherical polygon on the globe
231 * plane angle between radius vector to center and adjacent edge of
236 /* additional variables from snyder */
237 double q, Rprime, H, Ag, Azprime, Az, dprime, f, rho,
240 /* variables used to store intermediate results */
241 double cot_theta, tan_g, az_offset;
243 /* how many multiples of 60 degrees we adjust the azimuth */
244 int Az_adjust_multiples;
246 struct snyder_constants c;
249 * TODO by locality of reference, start by trying the same triangle
253 /* TODO put these constants in as radians to begin with */
254 c = constants[SNYDER_POLY_ICOSAHEDRON];
255 theta = c.theta * DEG2RAD;
259 for (i = 1; i <= 20; i++) {
261 struct isea_geo center;
263 center = icostriangles[i];
267 z = sph_distance(center.lon, center.lat, ll->lon, ll->lat);
269 z = acos(sin(center.lat) * sin(ll->lat)
270 + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));
273 /* not on this triangle */
274 if (z > g + 0.000005) { /* TODO DBL_EPSILON */
277 Az = sph_azimuth(ll->lon, ll->lat, center.lon, center.lat);
279 Az = atan2(cos(ll->lat) * sin(ll->lon - center.lon),
280 cos(center.lat) * sin(ll->lat)
281 - sin(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon)
286 /* This calculates "some" vertex coordinate */
287 az_offset = az_adjustment(i);
291 /* TODO I don't know why we do this. It's not in snyder */
292 /* maybe because we should have picked a better vertex */
297 * adjust Az for the point to fall within the range of 0 to
298 * 2(90 - theta) or 60 degrees for the hexagon, by
299 * and therefore 120 degrees for the triangle
301 * subtracting or adding multiples of 60 degrees to Az and
302 * recording the amount of adjustment
305 Az_adjust_multiples = 0;
308 Az_adjust_multiples--;
310 while (Az > DEG120 + DBL_EPSILON) {
312 Az_adjust_multiples++;
316 cot_theta = 1.0 / tan(theta);
317 tan_g = tan(g); /* TODO this is a constant */
319 /* Calculate q from eq 9. */
320 /* TODO cot_theta is cot(30) */
321 q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);
323 /* not in this triangle */
324 if (z > q + 0.000005) {
329 /* Apply equations 5-8 and 10-12 in order */
332 /* Rprime = 0.9449322893 * R; */
333 /* R' in the paper is for the truncated */
334 Rprime = 0.91038328153090290025;
337 H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));
340 /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */
341 Ag = Az + G + H - DEG180;
344 Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta);
347 /* cot(theta) = 1.73205080756887729355 */
348 dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);
351 f = dprime / (2.0 * Rprime * sin(q / 2.0));
354 rho = 2.0 * Rprime * f * sin(z / 2.0);
357 * add back the same 60 degree multiple adjustment from step
361 Azprime += DEG120 * Az_adjust_multiples;
363 /* calculate rectangular coordinates */
365 x = rho * sin(Azprime);
366 y = rho * cos(Azprime);
370 * translate coordinates to the origin for the particular
371 * hexagon on the flattened polyhedral map plot
381 * should be impossible, this implies that the coordinate is not on
385 fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",
386 ll->lon * RAD2DEG, ll->lat * RAD2DEG);
391 return 0; /* supresses a warning */
395 * return the new coordinates of any point in orginal coordinate system.
396 * Define a point (newNPold) in orginal coordinate system as the North Pole in
397 * new coordinate system, and the great circle connect the original and new
398 * North Pole as the lon0 longitude in new coordinate system, given any point
399 * in orginal coordinate system, this function return the new coordinates.
402 #define PRECISION 0.0000000000005
405 * TODO rename to isea_coordtrans and take an struct pointer to the result
406 * struct TODO remove function entirely. we use the synder one. the only
407 * difference is a 180 degree longitude rotation
411 coordtrans(struct isea_geo * newNPold, struct isea_geo * ptold, double lon0)
413 double cosptnewlat, cosptnewlon;
414 struct isea_geo ptnew;
416 cosptnewlat = sin(newNPold->lat) * sin(ptold->lat) +
417 cos(newNPold->lat) * cos(ptold->lat) * cos(newNPold->lon - ptold->lon);
419 if (cosptnewlat > 1.)
421 if (cosptnewlat < -1.)
424 ptnew.lat = acos(cosptnewlat);
426 if (fabs(ptnew.lat - 0.0) < PRECISION * 100000)
428 else if (fabs(ptnew.lat - M_PI) < PRECISION * 100000)
431 cosptnewlon = (sin(ptold->lat) * cos(newNPold->lat) - cos(ptold->lat) *
432 sin(newNPold->lat) * cos(newNPold->lon - ptold->lon)) / sin(ptnew.lat);
433 if (cosptnewlon > 1.)
435 if (cosptnewlon < -1.)
438 ptnew.lon = acos(cosptnewlon);
440 if ((ptold->lon - newNPold->lon) >= 0 && (ptold->lon - newNPold->lon) < 180)
441 ptnew.lon = -ptnew.lon + lon0;
444 ptnew.lon = ptnew.lon + lon0;
446 if (ptnew.lon > M_PI)
447 ptnew.lon -= 2 * M_PI;
448 if (ptnew.lon < -M_PI)
449 ptnew.lon += 2 * M_PI;
451 ptnew.lat = M_PI / 2 - ptnew.lat;
455 /* formula from Snyder, Map Projections: A working manual, p31 */
457 * old north pole at np in new coordinates
458 * could be simplified a bit with fewer intermediates
460 * TODO take a result pointer
464 snyder_ctran(struct isea_geo * np, struct isea_geo * pt)
467 double alpha, phi, lambda, lambda0, beta, lambdap, phip;
469 double lp_b; /* lambda prime minus beta */
482 sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);
486 /* use the two argument form so we end up in the right quadrant */
487 lp_b = atan2(cos_p * sin(lambda - lambda0),
488 (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));
490 lambdap = lp_b + beta;
492 /* normalize longitude */
493 /* TODO can we just do a modulus ? */
494 lambdap = fmod(lambdap, 2 * M_PI);
495 while (lambdap > M_PI)
497 while (lambdap < -M_PI)
500 phip = asin(sin_phip);
510 isea_ctran(struct isea_geo * np, struct isea_geo * pt, double lon0)
515 npt = snyder_ctran(np, pt);
518 npt.lon -= (M_PI - lon0 + np->lon);
521 * snyder is down tri 3, isea is along side of tri1 from vertex 0 to
522 * vertex 1 these are 180 degrees apart
525 /* normalize longitude */
526 npt.lon = fmod(npt.lon, 2 * M_PI);
527 while (npt.lon > M_PI)
529 while (npt.lon < -M_PI)
536 #define ISEA_STD_LAT 1.01722196792335072101
537 #define ISEA_STD_LON .19634954084936207740
540 struct isea_dgg isea_standard_dgg = {
541 20, 58.28252559, 11.25, 0.0,
545 /* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */
549 isea_grid_init(struct isea_dgg * g)
555 g->o_lat = ISEA_STD_LAT;
556 g->o_lon = ISEA_STD_LON;
568 isea_orient_isea(struct isea_dgg * g)
572 g->o_lat = ISEA_STD_LAT;
573 g->o_lon = ISEA_STD_LON;
580 isea_orient_pole(struct isea_dgg * g)
584 g->o_lat = M_PI / 2.0;
592 isea_orient_dymax(struct isea_dgg * g)
596 g->o_lat = 2.300882 * M_PI / 180.0;
597 g->o_lon = -5.24539 * M_PI / 180.0;
598 g->o_az = 7.46658 * M_PI / 180.0;
604 isea_transform(struct isea_dgg * g, struct isea_geo * in,
605 struct isea_pt * out)
607 struct isea_geo i, pole;
613 i = isea_ctran(&pole, in, g->o_az);
615 tri = isea_snyder_forward(&i, out);
623 #define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1)
625 /* get which triangle a point is in, assumes radius = 1.0 */
627 int isea_xy_triangle(struct isea_pt *p, double radius) {
630 struct isea_pt tc, pt;
631 pt.x = p->x / radius;
632 pt.y = p->y / radius;
634 tc = isea_triangle_xy(1);
635 dist = (pt.x - tc.x) * (pt.x - tc.x) + (pt.y - tc.y) * (pt.y - tc.y);
638 /* TODO just calculate it directly, no need to actually do all this
641 for (tri = 2; tri <= 20; tri++) {
642 tc = isea_triangle_xy(tri);
643 test = (pt.x - tc.x) * (pt.x - tc.x) + (pt.y - tc.y) * (pt.y - tc.y);
653 int isea_projtri(struct isea_pt *xy, double radius) {
656 tri = isea_xy_triangle(xy, radius);
661 xy->y += 2.0 * .14433756729740644112;
668 isea_rotate(struct isea_pt * pt, double degrees)
674 rad = -degrees * M_PI / 180.0;
675 while (rad >= 2.0 * M_PI) rad -= 2.0 * M_PI;
676 while (rad <= -2.0 * M_PI) rad += 2.0 * M_PI;
678 x = pt->x * cos(rad) + pt->y * sin(rad);
679 y = -pt->x * sin(rad) + pt->y * cos(rad);
686 int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {
687 struct isea_pt tc; /* center of triangle */
690 isea_rotate(pt, 180.0);
692 tc = isea_triangle_xy(tri);
702 int isea_projection(struct isea_dgg *g, double xin, double yin, double *xout,
707 /* convert projected triangle coords to quad xy coords, return quad number */
710 isea_ptdd(int tri, struct isea_pt *pt) {
713 downtri = (((tri - 1) / 5) % 2 == 1);
714 quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
716 isea_rotate(pt, downtri ? 240.0 : 60.0);
719 /* pt->y += cos(30.0 * M_PI / 180.0); */
720 pt->y += .86602540378443864672;
727 isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di)
731 double sidelength; /* in hexes */
736 /* This is the number of hexes from apex to base of a triangle */
737 sidelength = (pow(2.0, g->resolution) + 1.0) / 2.0;
739 /* apex to base is cos(30deg) */
740 hexwidth = cos(M_PI / 6.0) / sidelength;
742 /* TODO I think sidelength is always x.5, so
743 * (int)sidelength * 2 + 1 might be just as good
745 maxcoord = (int) (sidelength * 2.0 + 0.5);
748 hexbin2(0, hexwidth, v.x, v.y, &h.x, &h.y);
756 * you want to test for max coords for the next quad in the same
757 * "row" first to get the case where both are max
760 if (d == 0 && i == maxcoord) {
765 } else if (i == maxcoord) {
766 /* upper right in next quad */
772 } else if (d == maxcoord) {
773 /* lower right in quad to lower right */
777 } else if (quad >= 6) {
778 if (i == 0 && d == maxcoord) {
783 } else if (d == maxcoord) {
784 /* lower right in next quad */
790 } else if (i == maxcoord) {
791 /* upper right in quad to upper right */
792 quad = (quad - 4) % 5;
806 isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di) {
809 int sidelength; /* in hexes */
812 if (g->aperture == 3 && g->resolution % 2 != 0) {
813 return isea_dddi_ap3odd(g, quad, pt, di);
815 /* todo might want to do this as an iterated loop */
816 if (g->aperture >0) {
817 sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
819 sidelength = g->resolution;
822 hexwidth = 1.0 / sidelength;
825 isea_rotate(&v, -30.0);
826 hexbin2(0, hexwidth, v.x, v.y, &h.x, &h.y);
830 /* we may actually be on another quad */
832 if (h.x == 0 && h.z == -sidelength) {
838 } else if (h.z == -sidelength) {
842 h.y = sidelength - h.x;
843 h.z = h.x - sidelength;
845 } else if (h.x == sidelength) {
850 } else if (quad >= 6) {
851 if (h.z == 0 && h.x == sidelength) {
857 } else if (h.x == sidelength) {
861 h.x = h.y + sidelength;
864 } else if (h.y == -sidelength) {
878 int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
879 struct isea_pt *di) {
884 quad = isea_ptdd(tri, &v);
885 quad = isea_dddi(g, quad, &v, di);
891 int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
900 /* hexes in a quad */
901 hexes = (int) (pow(g->aperture, g->resolution) + 0.5);
903 g->serial = 1 + 10 * hexes + 1;
906 if (g->aperture == 3 && g->resolution % 2 == 1) {
907 height = (int) (pow(g->aperture, (g->resolution - 1) / 2.0));
908 sn = ((int) di->x) * height;
909 sn += ((int) di->y) / height;
910 sn += (quad - 1) * hexes;
913 sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
914 sn = (quad - 1) * hexes + sidelength * di->x + di->y + 2;
921 /* TODO just encode the quad in the d or i coordinate
922 * quad is 0-11, which can be four bits.
923 * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf
925 /* convert a q2di to global hex coord */
927 int isea_hex(struct isea_dgg *g, int tri,
928 struct isea_pt *pt, struct isea_pt *hex) {
931 int d, i, x, y, quad;
934 quad = isea_ptdi(g, tri, pt, &v);
936 hex->x = ((int)v.x << 4) + quad;
944 /* Aperture 3 odd resolutions */
945 if (g->aperture == 3 && g->resolution % 2 != 0) {
946 int offset = (int)(pow(3.0, g->resolution - 1) + 0.5);
948 oddsl = (pow(2.0, g->resolution) + 1.0) / 2.0;
950 d += offset * ((g->quad-1) % 5);
951 i += offset * ((g->quad-1) % 5);
956 } else if (quad == 11) {
959 } else if (quad > 5) {
966 hex->x = x + offset / 3;
967 hex->y = y + 2 * offset / 3;
971 /* aperture 3 even resolutions and aperture 4 */
972 sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
976 } else if (g->quad == 11) {
977 hex->x = sidelength * 2;
980 hex->x = d + sidelength * ((g->quad-1) % 5);
981 if (g->quad > 5) hex->x += sidelength;
982 hex->y = i + sidelength * ((g->quad-1) % 5);
990 isea_forward(struct isea_dgg *g, struct isea_geo *in)
992 int tri, downtri, quad;
993 struct isea_pt out, coord;
995 tri = isea_transform(g, in, &out);
997 downtri = (((tri - 1) / 5) % 2 == 1);
998 quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
1000 if (g->output == ISEA_PLANE) {
1001 isea_tri_plane(tri, &out, g->radius);
1005 /* convert to isea standard triangle size */
1006 out.x = out.x / g->radius * ISEA_SCALE;
1007 out.y = out.y / g->radius * ISEA_SCALE;
1009 out.y += 2.0 * .14433756729740644112;
1011 switch (g->output) {
1013 /* nothing to do, already in projected triangle */
1015 case ISEA_VERTEX2DD:
1016 g->quad = isea_ptdd(tri, &out);
1019 /* Same as above, we just don't print as much */
1020 g->quad = isea_ptdd(tri, &out);
1023 g->quad = isea_ptdi(g, tri, &out, &coord);
1027 isea_ptdi(g, tri, &out, &coord);
1028 /* disn will set g->serial */
1029 isea_disn(g, g->quad, &coord);
1033 isea_hex(g, tri, &out, &coord);