diff options
| author | Even Rouault <even.rouault@spatialys.com> | 2018-12-19 13:00:37 +0100 |
|---|---|---|
| committer | Even Rouault <even.rouault@spatialys.com> | 2018-12-26 10:08:55 +0100 |
| commit | df574ae332d57f556fd56314883b3354cab1d0ff (patch) | |
| tree | 63a68d40d7ed7932d6329d9c7baa340bb6294f7f /src/projections/ob_tran.cpp | |
| parent | e6de172371ea203f6393d745641d66c82b5b13e2 (diff) | |
| download | PROJ-df574ae332d57f556fd56314883b3354cab1d0ff.tar.gz PROJ-df574ae332d57f556fd56314883b3354cab1d0ff.zip | |
cpp conversion: remove useless pj_, PJ_ and proj_ filename prefixes
Diffstat (limited to 'src/projections/ob_tran.cpp')
| -rw-r--r-- | src/projections/ob_tran.cpp | 245 |
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; +} |
