aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am1
-rw-r--r--src/lib_proj.cmake1
-rw-r--r--src/pj_list.h1
-rw-r--r--src/projections/igh.cpp93
-rw-r--r--src/projections/igh_o.cpp268
5 files changed, 331 insertions, 33 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index 9f24a60d..3f5b2968 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -150,6 +150,7 @@ libproj_la_SOURCES = \
projections/gn_sinu.cpp \
projections/goode.cpp \
projections/igh.cpp \
+ projections/igh_o.cpp \
projections/hatano.cpp \
projections/loxim.cpp \
projections/mbt_fps.cpp \
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index c318a266..4eca0256 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -134,6 +134,7 @@ set(SRC_LIBPROJ_PROJECTIONS
projections/gn_sinu.cpp
projections/goode.cpp
projections/igh.cpp
+ projections/igh_o.cpp
projections/hatano.cpp
projections/loxim.cpp
projections/mbt_fps.cpp
diff --git a/src/pj_list.h b/src/pj_list.h
index cf81ae0e..c7e6d209 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -69,6 +69,7 @@ PROJ_HEAD(helmert, "3- and 7-parameter Helmert shift")
PROJ_HEAD(hgridshift, "Horizontal grid shift")
PROJ_HEAD(horner, "Horner polynomial evaluation")
PROJ_HEAD(igh, "Interrupted Goode Homolosine")
+PROJ_HEAD(igh_o, "Interrupted Goode Homolosine Oceanic View")
PROJ_HEAD(imw_p, "International Map of the World Polyconic")
PROJ_HEAD(isea, "Icosahedral Snyder Equal Area")
PROJ_HEAD(kav5, "Kavraisky V")
diff --git a/src/projections/igh.cpp b/src/projections/igh.cpp
index 8a41cea3..2be94889 100644
--- a/src/projections/igh.cpp
+++ b/src/projections/igh.cpp
@@ -8,15 +8,31 @@
PROJ_HEAD(igh, "Interrupted Goode Homolosine") "\n\tPCyl, Sph";
+/*
+This projection is a compilation of 12 separate sub-projections.
+Sinusoidal projections are found near the equator and Mollweide
+projections are found at higher latitudes. The transition between
+the two occurs at 40 degrees latitude and is represented by the
+constant `phi_boundary`.
+
+Each sub-projection is assigned an integer label
+numbered 1 through 12. Most of this code contains logic to assign
+the labels based on latitude (phi) and longitude (lam) regions.
+
+Original Reference:
+J. Paul Goode (1925) THE HOMOLOSINE PROJECTION: A NEW DEVICE FOR
+ PORTRAYING THE EARTH'S SURFACE ENTIRE, Annals of the Association of
+ American Geographers, 15:3, 119-125, DOI: 10.1080/00045602509356949
+*/
+
C_NAMESPACE PJ *pj_sinu(PJ *), *pj_moll(PJ *);
-/* 40d 44' 11.8" [degrees] */
-/*
-static const double d4044118 = (40 + 44/60. + 11.8/3600.) * DEG_TO_RAD;
-has been replaced by this define, to eliminate portability issue:
-Initializer element not computable at load time
+/*
+Transition from sinusoidal to Mollweide projection
+Latitude (phi): 40deg 44' 11.8"
*/
-#define d4044118 ((40 + 44/60. + 11.8/3600.) * DEG_TO_RAD)
+
+static const double phi_boundary = (40 + 44/60. + 11.8/3600.) * DEG_TO_RAD;
static const double d10 = 10 * DEG_TO_RAD;
static const double d20 = 20 * DEG_TO_RAD;
@@ -46,13 +62,13 @@ static PJ_XY igh_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward
struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
int z;
- if (lp.phi >= d4044118) { /* 1|2 */
+ if (lp.phi >= phi_boundary) { /* 1|2 */
z = (lp.lam <= -d40 ? 1: 2);
}
else if (lp.phi >= 0) { /* 3|4 */
z = (lp.lam <= -d40 ? 3: 4);
}
- else if (lp.phi >= -d4044118) { /* 5|6|7|8 */
+ else if (lp.phi >= -phi_boundary) { /* 5|6|7|8 */
if (lp.lam <= -d100) z = 5; /* 5 */
else if (lp.lam <= -d20) z = 6; /* 6 */
else if (lp.lam <= d80) z = 7; /* 7 */
@@ -82,11 +98,11 @@ static PJ_LP igh_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse
int z = 0;
if (xy.y > y90+EPSLN || xy.y < -y90+EPSLN) /* 0 */
z = 0;
- else if (xy.y >= d4044118) /* 1|2 */
+ else if (xy.y >= phi_boundary) /* 1|2 */
z = (xy.x <= -d40? 1: 2);
else if (xy.y >= 0) /* 3|4 */
z = (xy.x <= -d40? 3: 4);
- else if (xy.y >= -d4044118) { /* 5|6|7|8 */
+ else if (xy.y >= -phi_boundary) { /* 5|6|7|8 */
if (xy.x <= -d100) z = 5; /* 5 */
else if (xy.x <= -d20) z = 6; /* 6 */
else if (xy.x <= d80) z = 7; /* 7 */
@@ -100,7 +116,7 @@ static PJ_LP igh_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse
}
if (z) {
- int ok = 0;
+ bool ok = false;
xy.x -= Q->pj[z-1]->x0;
xy.y -= Q->pj[z-1]->y0;
@@ -145,9 +161,11 @@ static PJ *destructor (PJ *P, int errlev) {
if (nullptr==P->opaque)
return pj_default_destructor (P, errlev);
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+
for (i = 0; i < 12; ++i) {
- if (static_cast<struct pj_opaque*>(P->opaque)->pj[i])
- static_cast<struct pj_opaque*>(P->opaque)->pj[i]->destructor(static_cast<struct pj_opaque*>(P->opaque)->pj[i], errlev);
+ if (Q->pj[i])
+ Q->pj[i]->destructor(Q->pj[i], errlev);
}
return pj_default_destructor(P, errlev);
@@ -175,18 +193,21 @@ static PJ *destructor (PJ *P, int errlev) {
-180 -100 -20 80 180
*/
-#define SETUP(n, proj, x_0, y_0, lon_0) \
- if (!(Q->pj[n-1] = pj_##proj(nullptr))) return destructor(P, ENOMEM); \
- if (!(Q->pj[n-1] = pj_##proj(Q->pj[n-1]))) return destructor(P, ENOMEM); \
- Q->pj[n-1]->ctx = P->ctx; \
- Q->pj[n-1]->x0 = x_0; \
- Q->pj[n-1]->y0 = y_0; \
+static bool setup_zone(PJ *P, struct pj_opaque *Q, int n,
+ PJ*(*proj_ptr)(PJ*), double x_0,
+ double y_0, double lon_0) {
+ if (!(Q->pj[n-1] = proj_ptr(nullptr))) return false;
+ if (!(Q->pj[n-1] = proj_ptr(Q->pj[n-1]))) return false;
+ Q->pj[n-1]->ctx = P->ctx;
+ Q->pj[n-1]->x0 = x_0;
+ Q->pj[n-1]->y0 = y_0;
Q->pj[n-1]->lam0 = lon_0;
-
+ return true;
+}
PJ *PROJECTION(igh) {
PJ_XY xy1, xy3;
- PJ_LP lp = { 0, d4044118 };
+ PJ_LP lp = { 0, phi_boundary };
struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
if (nullptr==Q)
return pj_default_destructor (P, ENOMEM);
@@ -194,15 +215,18 @@ PJ *PROJECTION(igh) {
/* sinusoidal zones */
- SETUP(3, sinu, -d100, 0, -d100);
- SETUP(4, sinu, d30, 0, d30);
- SETUP(5, sinu, -d160, 0, -d160);
- SETUP(6, sinu, -d60, 0, -d60);
- SETUP(7, sinu, d20, 0, d20);
- SETUP(8, sinu, d140, 0, d140);
+ if (!setup_zone(P, Q, 3, pj_sinu, -d100, 0, -d100) ||
+ !setup_zone(P, Q, 4, pj_sinu, d30, 0, d30) ||
+ !setup_zone(P, Q, 5, pj_sinu, -d160, 0, -d160) ||
+ !setup_zone(P, Q, 6, pj_sinu, -d60, 0, -d60) ||
+ !setup_zone(P, Q, 7, pj_sinu, d20, 0, d20) ||
+ !setup_zone(P, Q, 8, pj_sinu, d140, 0, d140))
+ {
+ return destructor(P, ENOMEM);
+ }
/* mollweide zones */
- SETUP(1, moll, -d100, 0, -d100);
+ setup_zone(P, Q, 1, pj_moll, -d100, 0, -d100);
/* y0 ? */
xy1 = Q->pj[0]->fwd(lp, Q->pj[0]); /* zone 1 */
@@ -213,11 +237,14 @@ PJ *PROJECTION(igh) {
Q->pj[0]->y0 = Q->dy0;
/* mollweide zones (cont'd) */
- SETUP( 2, moll, d30, Q->dy0, d30);
- SETUP( 9, moll, -d160, -Q->dy0, -d160);
- SETUP(10, moll, -d60, -Q->dy0, -d60);
- SETUP(11, moll, d20, -Q->dy0, d20);
- SETUP(12, moll, d140, -Q->dy0, d140);
+ if (!setup_zone(P, Q, 2, pj_moll, d30, Q->dy0, d30) ||
+ !setup_zone(P, Q, 9, pj_moll, -d160, -Q->dy0, -d160) ||
+ !setup_zone(P, Q,10, pj_moll, -d60, -Q->dy0, -d60) ||
+ !setup_zone(P, Q,11, pj_moll, d20, -Q->dy0, d20) ||
+ !setup_zone(P, Q,12, pj_moll, d140, -Q->dy0, d140))
+ {
+ return destructor(P, ENOMEM);
+ }
P->inv = igh_s_inverse;
P->fwd = igh_s_forward;
diff --git a/src/projections/igh_o.cpp b/src/projections/igh_o.cpp
new file mode 100644
index 00000000..785a8784
--- /dev/null
+++ b/src/projections/igh_o.cpp
@@ -0,0 +1,268 @@
+#define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+
+#include "proj.h"
+#include "proj_internal.h"
+
+PROJ_HEAD(igh_o, "Interrupted Goode Homolosine Oceanic View") "\n\tPCyl, Sph";
+
+/*
+This projection is a variant of the Interrupted Goode Homolosine
+projection that emphasizes ocean areas. The projection is a
+compilation of 12 separate sub-projections. Sinusoidal projections
+are found near the equator and Mollweide projections are found at
+higher latitudes. The transition between the two occurs at 40 degrees
+latitude and is represented by `phi_boundary`.
+
+Each sub-projection is assigned an integer label
+numbered 1 through 12. Most of this code contains logic to assign
+the labels based on latitude (phi) and longitude (lam) regions.
+
+Original Reference:
+J. Paul Goode (1925) THE HOMOLOSINE PROJECTION: A NEW DEVICE FOR
+ PORTRAYING THE EARTH'S SURFACE ENTIRE, Annals of the Association of
+ American Geographers, 15:3, 119-125, DOI: 10.1080/00045602509356949
+*/
+
+C_NAMESPACE PJ *pj_sinu(PJ *), *pj_moll(PJ *);
+
+/*
+Transition from sinusoidal to Mollweide projection
+Latitude (phi): 40deg 44' 11.8"
+*/
+
+static const double phi_boundary = (40 + 44/60. + 11.8/3600.) * DEG_TO_RAD;
+
+static const double d10 = 10 * DEG_TO_RAD;
+static const double d20 = 20 * DEG_TO_RAD;
+static const double d50 = 50 * DEG_TO_RAD;
+static const double d60 = 60 * DEG_TO_RAD;
+static const double d90 = 90 * DEG_TO_RAD;
+static const double d110 = 110 * DEG_TO_RAD;
+static const double d140 = 140 * DEG_TO_RAD;
+static const double d150 = 150 * DEG_TO_RAD;
+static const double d160 = 160 * DEG_TO_RAD;
+static const double d130 = 130 * DEG_TO_RAD;
+static const double d180 = 180 * DEG_TO_RAD;
+
+static const double EPSLN = 1.e-10; /* allow a little 'slack' on zone edge positions */
+
+namespace { // anonymous namespace
+struct pj_opaque {
+ struct PJconsts* pj[12]; \
+ double dy0;
+};
+} // anonymous namespace
+
+
+/*
+Assign an integer index representing each of the 12
+sub-projection zones based on latitude (phi) and
+longitude (lam) ranges.
+*/
+
+static PJ_XY igh_o_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */
+ PJ_XY xy;
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+ int z;
+
+ if (lp.phi >= phi_boundary) {
+ if (lp.lam <= -d90) z = 1;
+ else if (lp.lam >= d60) z = 3;
+ else z = 2;
+ }
+ else if (lp.phi >= 0) {
+ if (lp.lam <= -d90) z = 4;
+ else if (lp.lam >= d60) z = 6;
+ else z = 5;
+ }
+ else if (lp.phi >= -phi_boundary) {
+ if (lp.lam <= -d60) z = 7;
+ else if (lp.lam >= d90) z = 9;
+ else z = 8;
+ }
+ else {
+ if (lp.lam <= -d60) z = 10;
+ else if (lp.lam >= d90) z = 12;
+ else z = 11;
+ }
+
+ lp.lam -= Q->pj[z-1]->lam0;
+ xy = Q->pj[z-1]->fwd(lp, Q->pj[z-1]);
+ xy.x += Q->pj[z-1]->x0;
+ xy.y += Q->pj[z-1]->y0;
+
+ return xy;
+}
+
+
+static PJ_LP igh_o_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
+ PJ_LP lp = {0.0,0.0};
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+ const double y90 = Q->dy0 + sqrt(2.0); /* lt=90 corresponds to y=y0+sqrt(2) */
+
+ int z = 0;
+ if (xy.y > y90+EPSLN || xy.y < -y90+EPSLN) /* 0 */
+ z = 0;
+ else if (xy.y >= phi_boundary)
+ if (xy.x <= -d90) z = 1;
+ else if (xy.x >= d60) z = 3;
+ else z = 2;
+ else if (xy.y >= 0)
+ if (xy.x <= -d90) z = 4;
+ else if (xy.x >= d60) z = 6;
+ else z = 5;
+ else if (xy.y >= -phi_boundary) {
+ if (xy.x <= -d60) z = 7;
+ else if (xy.x >= d90) z = 9;
+ else z = 8;
+ }
+ else {
+ if (xy.x <= -d60) z = 10;
+ else if (xy.x >= d90) z = 12;
+ else z = 11;
+ }
+
+ if (z) {
+ bool ok = false;
+
+ xy.x -= Q->pj[z-1]->x0;
+ xy.y -= Q->pj[z-1]->y0;
+ lp = Q->pj[z-1]->inv(xy, Q->pj[z-1]);
+ lp.lam += Q->pj[z-1]->lam0;
+
+ switch (z) {
+ /* Plot projectable ranges with exetension lobes in zones 1 & 3 */
+ case 1: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d90+EPSLN) ||
+ ((lp.lam >= d160-EPSLN && lp.lam <= d180+EPSLN) &&
+ (lp.phi >= d50-EPSLN && lp.phi <= d90+EPSLN)); break;
+ case 2: ok = (lp.lam >= -d90-EPSLN && lp.lam <= d60+EPSLN); break;
+ case 3: ok = (lp.lam >= d60-EPSLN && lp.lam <= d180+EPSLN) ||
+ ((lp.lam >= -d180-EPSLN && lp.lam <= -d160+EPSLN) &&
+ (lp.phi >= d50-EPSLN && lp.phi <= d90+EPSLN)); break;
+ case 4: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d90+EPSLN); break;
+ case 5: ok = (lp.lam >= -d90-EPSLN && lp.lam <= d60+EPSLN); break;
+ case 6: ok = (lp.lam >= d60-EPSLN && lp.lam <= d180+EPSLN); break;
+ case 7: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d60+EPSLN); break;
+ case 8: ok = (lp.lam >= -d60-EPSLN && lp.lam <= d90+EPSLN); break;
+ case 9: ok = (lp.lam >= d90-EPSLN && lp.lam <= d180+EPSLN); break;
+ case 10: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d60+EPSLN); break;
+ case 11: ok = (lp.lam >= -d60-EPSLN && lp.lam <= d90+EPSLN); break;
+ case 12: ok = (lp.lam >= d90-EPSLN && lp.lam <= d180+EPSLN); break;
+
+ }
+ z = (!ok? 0: z); /* projectable? */
+ }
+
+ if (!z) lp.lam = HUGE_VAL;
+ if (!z) lp.phi = HUGE_VAL;
+
+ return lp;
+}
+
+
+static PJ *destructor (PJ *P, int errlev) {
+ int i;
+ if (nullptr==P)
+ return nullptr;
+
+ if (nullptr==P->opaque)
+ return pj_default_destructor (P, errlev);
+
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+
+ for (i = 0; i < 12; ++i) {
+ if (Q->pj[i])
+ Q->pj[i]->destructor(Q->pj[i], errlev);
+ }
+
+ return pj_default_destructor(P, errlev);
+}
+
+
+
+/*
+ Zones:
+
+ -180 -90 60 180
+ +---------+----------------+-------------+ Zones 1,2,3,10,11 & 12:
+ |1 |2 |3 | Mollweide projection
+ | | | |
+ +---------+----------------+-------------+ Zones 4,5,6,7,8 & 9:
+ |4 |5 |6 | Sinusoidal projection
+ | | | |
+ 0 +---------+--+-------------+--+----------+
+ |7 |8 |9 |
+ | | | |
+ +------------+----------------+----------+
+ |10 |11 |12 |
+ | | | |
+ +------------+----------------+----------+
+ -180 -60 90 180
+*/
+
+static bool setup_zone(PJ *P, struct pj_opaque *Q, int n,
+ PJ*(*proj_ptr)(PJ*), double x_0,
+ double y_0, double lon_0) {
+ if (!(Q->pj[n-1] = proj_ptr(nullptr))) return false;
+ if (!(Q->pj[n-1] = proj_ptr(Q->pj[n-1]))) return false;
+ Q->pj[n-1]->ctx = P->ctx;
+ Q->pj[n-1]->x0 = x_0;
+ Q->pj[n-1]->y0 = y_0;
+ Q->pj[n-1]->lam0 = lon_0;
+ return true;
+}
+
+PJ *PROJECTION(igh_o) {
+ PJ_XY xy1, xy4;
+ PJ_LP lp = { 0, phi_boundary };
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
+ if (nullptr==Q)
+ return pj_default_destructor (P, ENOMEM);
+ P->opaque = Q;
+
+
+ /* sinusoidal zones */
+ if (!setup_zone(P, Q, 4, pj_sinu, -d140, 0, -d140) ||
+ !setup_zone(P, Q, 5, pj_sinu, d10, 0, d10) ||
+ !setup_zone(P, Q, 6, pj_sinu, d130, 0, d130) ||
+ !setup_zone(P, Q, 7, pj_sinu, -d110, 0, -d110) ||
+ !setup_zone(P, Q, 8, pj_sinu, d20, 0, d20) ||
+ !setup_zone(P, Q, 9, pj_sinu, d150, 0, d150))
+ {
+ return destructor(P, ENOMEM);
+ }
+
+
+ /* mollweide zones */
+ if (!setup_zone(P, Q, 1, pj_moll, -d140, 0, -d140)) {
+ return destructor(P, ENOMEM);
+ }
+
+ /* y0 ? */
+ xy1 = Q->pj[0]->fwd(lp, Q->pj[0]); /* zone 1 */
+ xy4 = Q->pj[3]->fwd(lp, Q->pj[3]); /* zone 4 */
+ /* y0 + xy1.y = xy4.y for lt = 40d44'11.8" */
+ Q->dy0 = xy4.y - xy1.y;
+
+ Q->pj[0]->y0 = Q->dy0;
+
+ /* mollweide zones (cont'd) */
+ if (!setup_zone(P, Q, 2, pj_moll, d10, Q->dy0, d10) ||
+ !setup_zone(P, Q, 3, pj_moll, d130, Q->dy0, d130) ||
+ !setup_zone(P, Q, 10, pj_moll, -d110, -Q->dy0, -d110) ||
+ !setup_zone(P, Q, 11, pj_moll, d20, -Q->dy0, d20) ||
+ !setup_zone(P, Q, 12, pj_moll, d150, -Q->dy0, d150))
+ {
+ return destructor(P, ENOMEM);
+ }
+
+ P->inv = igh_o_s_inverse;
+ P->fwd = igh_o_s_forward;
+ P->destructor = destructor;
+ P->es = 0.;
+
+ return P;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-16 08:30:09 +0000