aboutsummaryrefslogtreecommitdiff
path: root/src/projections/ob_tran.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/projections/ob_tran.cpp')
-rw-r--r--src/projections/ob_tran.cpp245
1 files changed, 245 insertions, 0 deletions
diff --git a/src/projections/ob_tran.cpp b/src/projections/ob_tran.cpp
new file mode 100644
index 00000000..d34059a9
--- /dev/null
+++ b/src/projections/ob_tran.cpp
@@ -0,0 +1,245 @@
+#define PJ_LIB__
+#include <errno.h>
+#include <math.h>
+#include <stddef.h>
+#include <string.h>
+
+#include "proj.h"
+#include "projects.h"
+
+namespace { // anonymous namespace
+struct pj_opaque {
+ struct PJconsts *link;
+ double lamp;
+ double cphip, sphip;
+};
+} // anonymous namespace
+
+PROJ_HEAD(ob_tran, "General Oblique Transformation") "\n\tMisc Sph"
+"\n\to_proj= plus parameters for projection"
+"\n\to_lat_p= o_lon_p= (new pole) or"
+"\n\to_alpha= o_lon_c= o_lat_c= or"
+"\n\to_lon_1= o_lat_1= o_lon_2= o_lat_2=";
+
+#define TOL 1e-10
+
+
+static XY o_forward(LP lp, PJ *P) { /* spheroid */
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+ double coslam, sinphi, cosphi;
+
+ coslam = cos(lp.lam);
+ sinphi = sin(lp.phi);
+ cosphi = cos(lp.phi);
+ lp.lam = adjlon(aatan2(cosphi * sin(lp.lam), Q->sphip * cosphi * coslam +
+ Q->cphip * sinphi) + Q->lamp);
+ lp.phi = aasin(P->ctx,Q->sphip * sinphi - Q->cphip * cosphi * coslam);
+
+ return Q->link->fwd(lp, Q->link);
+}
+
+
+static XY t_forward(LP lp, PJ *P) { /* spheroid */
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+ double cosphi, coslam;
+
+ cosphi = cos(lp.phi);
+ coslam = cos(lp.lam);
+ lp.lam = adjlon(aatan2(cosphi * sin(lp.lam), sin(lp.phi)) + Q->lamp);
+ lp.phi = aasin(P->ctx, - cosphi * coslam);
+
+ return Q->link->fwd(lp, Q->link);
+}
+
+
+static LP o_inverse(XY xy, PJ *P) { /* spheroid */
+
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+ double coslam, sinphi, cosphi;
+
+ LP lp = Q->link->inv(xy, Q->link);
+ if (lp.lam != HUGE_VAL) {
+ coslam = cos(lp.lam -= Q->lamp);
+ sinphi = sin(lp.phi);
+ cosphi = cos(lp.phi);
+ lp.phi = aasin(P->ctx,Q->sphip * sinphi + Q->cphip * cosphi * coslam);
+ lp.lam = aatan2(cosphi * sin(lp.lam), Q->sphip * cosphi * coslam -
+ Q->cphip * sinphi);
+ }
+ return lp;
+}
+
+
+static LP t_inverse(XY xy, PJ *P) { /* spheroid */
+
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+ double cosphi, t;
+
+ LP lp = Q->link->inv(xy, Q->link);
+ if (lp.lam != HUGE_VAL) {
+ cosphi = cos(lp.phi);
+ t = lp.lam - Q->lamp;
+ lp.lam = aatan2(cosphi * sin(t), - sin(lp.phi));
+ lp.phi = aasin(P->ctx,cosphi * cos(t));
+ }
+ return lp;
+}
+
+
+static PJ *destructor(PJ *P, int errlev) {
+ if (nullptr==P)
+ return nullptr;
+ if (nullptr==P->opaque)
+ return pj_default_destructor (P, errlev);
+
+ if (static_cast<struct pj_opaque*>(P->opaque)->link)
+ static_cast<struct pj_opaque*>(P->opaque)->link->destructor (static_cast<struct pj_opaque*>(P->opaque)->link, errlev);
+
+ return pj_default_destructor(P, errlev);
+}
+
+
+
+
+/***********************************************************************
+
+These functions are modified versions of the functions "argc_params"
+and "argv_params" from PJ_pipeline.c
+
+Basically, they do the somewhat backwards stunt of turning the paralist
+representation of the +args back into the original +argv, +argc
+representation accepted by pj_init_ctx().
+
+This, however, also begs the question of whether we really need the
+paralist linked list representation, or if we could do with a simpler
+null-terminated argv style array? This would simplfy some code, and
+keep memory allocations more localized.
+
+***********************************************************************/
+
+typedef struct {int argc; char **argv;} ARGS;
+
+/* count the number of args in the linked list <params> */
+static size_t paralist_params_argc (paralist *params) {
+ size_t argc = 0;
+ for (; params != nullptr; params = params->next)
+ argc++;
+ return argc;
+}
+
+
+/* turn paralist into argc/argv style argument list */
+static ARGS ob_tran_target_params (paralist *params) {
+ int i = 0;
+ ARGS args = {0, nullptr};
+ size_t argc = paralist_params_argc (params);
+ if (argc < 2)
+ return args;
+
+ /* all args except the proj_ob_tran */
+ args.argv = static_cast<char**>(pj_calloc (argc - 1, sizeof (char *)));
+ if (nullptr==args.argv)
+ return args;
+
+ /* Copy all args *except* the proj=ob_tran arg to the argv array */
+ for (i = 0; params != nullptr; params = params->next) {
+ if (0==strcmp (params->param, "proj=ob_tran"))
+ continue;
+ args.argv[i++] = params->param;
+ }
+ args.argc = i;
+
+ /* Then convert the o_proj=xxx element to proj=xxx */
+ for (i = 0; i < args.argc; i++) {
+ if (0!=strncmp (args.argv[i], "o_proj=", 7))
+ continue;
+ args.argv[i] += 2;
+ break;
+ }
+
+ return args;
+}
+
+
+
+PJ *PROJECTION(ob_tran) {
+ double phip;
+ char *name;
+ ARGS args;
+ PJ *R; /* projection to rotate */
+
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
+ if (nullptr==Q)
+ return destructor(P, ENOMEM);
+
+ P->opaque = Q;
+ P->destructor = destructor;
+
+ /* get name of projection to be translated */
+ if (!(name = pj_param(P->ctx, P->params, "so_proj").s))
+ return destructor(P, PJD_ERR_NO_ROTATION_PROJ);
+
+ /* avoid endless recursion */
+ if( strcmp(name, "ob_tran") == 0 )
+ return destructor(P, PJD_ERR_FAILED_TO_FIND_PROJ);
+
+ /* Create the target projection object to rotate */
+ args = ob_tran_target_params (P->params);
+ R = pj_init_ctx (pj_get_ctx(P), args.argc, args.argv);
+ pj_dealloc (args.argv);
+
+ if (nullptr==R)
+ return destructor (P, PJD_ERR_UNKNOWN_PROJECTION_ID);
+ Q->link = R;
+
+ if (pj_param(P->ctx, P->params, "to_alpha").i) {
+ double lamc, phic, alpha;
+
+ lamc = pj_param(P->ctx, P->params, "ro_lon_c").f;
+ phic = pj_param(P->ctx, P->params, "ro_lat_c").f;
+ alpha = pj_param(P->ctx, P->params, "ro_alpha").f;
+
+ if (fabs(fabs(phic) - M_HALFPI) <= TOL)
+ return destructor(P, PJD_ERR_LAT_0_OR_ALPHA_EQ_90);
+
+ Q->lamp = lamc + aatan2(-cos(alpha), -sin(alpha) * sin(phic));
+ phip = aasin(P->ctx,cos(phic) * sin(alpha));
+ } else if (pj_param(P->ctx, P->params, "to_lat_p").i) { /* specified new pole */
+ Q->lamp = pj_param(P->ctx, P->params, "ro_lon_p").f;
+ phip = pj_param(P->ctx, P->params, "ro_lat_p").f;
+ } else { /* specified new "equator" points */
+ double lam1, lam2, phi1, phi2, con;
+
+ lam1 = pj_param(P->ctx, P->params, "ro_lon_1").f;
+ phi1 = pj_param(P->ctx, P->params, "ro_lat_1").f;
+ lam2 = pj_param(P->ctx, P->params, "ro_lon_2").f;
+ phi2 = pj_param(P->ctx, P->params, "ro_lat_2").f;
+ if (fabs(phi1 - phi2) <= TOL || (con = fabs(phi1)) <= TOL ||
+ fabs(con - M_HALFPI) <= TOL || fabs(fabs(phi2) - M_HALFPI) <= TOL)
+ return destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
+
+ Q->lamp = atan2(cos(phi1) * sin(phi2) * cos(lam1) -
+ sin(phi1) * cos(phi2) * cos(lam2),
+ sin(phi1) * cos(phi2) * sin(lam2) -
+ cos(phi1) * sin(phi2) * sin(lam1));
+ phip = atan(-cos(Q->lamp - lam1) / tan(phi1));
+ }
+
+ if (fabs(phip) > TOL) { /* oblique */
+ Q->cphip = cos(phip);
+ Q->sphip = sin(phip);
+ P->fwd = Q->link->fwd ? o_forward : nullptr;
+ P->inv = Q->link->inv ? o_inverse : nullptr;
+ } else { /* transverse */
+ P->fwd = Q->link->fwd ? t_forward : nullptr;
+ P->inv = Q->link->inv ? t_inverse : nullptr;
+ }
+
+ /* Support some rather speculative test cases, where the rotated projection */
+ /* is actually latlong. We do not want scaling in that case... */
+ if (Q->link->right==PJ_IO_UNITS_ANGULAR)
+ P->right = PJ_IO_UNITS_PROJECTED;
+
+
+ return P;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 14:46:32 +0000