aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/PJ_helmert.c32
-rw-r--r--src/gie.c47
-rw-r--r--src/pj_fwd.c115
-rw-r--r--src/pj_geocent.c2
-rw-r--r--src/pj_init.c23
-rw-r--r--src/pj_inv.c123
-rw-r--r--src/pj_malloc.c34
-rw-r--r--src/proj.def2
-rw-r--r--src/proj_4D_api.c119
-rw-r--r--src/proj_internal.h4
-rw-r--r--src/projects.h10
-rw-r--r--test/gie/4D-API_cs2cs-style.gie161
-rw-r--r--test/gie/axisswap.gie11
-rw-r--r--test/gie/builtins.gie16
-rw-r--r--test/gie/more_builtins.gie10
15 files changed, 603 insertions, 106 deletions
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
index f726d946..229e30c2 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -486,6 +486,25 @@ PJ *TRANSFORMATION(helmert, 0) {
P->right = PJ_IO_UNITS_PROJECTED;
}
+ /* Support the classic PROJ towgs84 parameter, but allow later overrides.*/
+ /* Note that if towgs84 is specified, the datum_params array is set up */
+ /* for us automagically by the pj_datum_set call in pj_init_ctx */
+ if (pj_param_exists (P->params, "towgs84")) {
+ Q->xyz_0.x = P->datum_params[0];
+ Q->xyz_0.y = P->datum_params[1];
+ Q->xyz_0.z = P->datum_params[2];
+
+ Q->opk_0.o = P->datum_params[3];
+ Q->opk_0.p = P->datum_params[4];
+ Q->opk_0.k = P->datum_params[5];
+
+ /* We must undo conversion to absolute scale from pj_datum_set */
+ if (0==P->datum_params[6])
+ Q->scale_0 = 0;
+ else
+ Q->scale_0 = (P->datum_params[6] - 1) * 1e6;
+ }
+
/* Translations */
if (pj_param (P->ctx, P->params, "tx").i)
Q->xyz_0.x = pj_param (P->ctx, P->params, "dx").f;
@@ -570,14 +589,13 @@ PJ *TRANSFORMATION(helmert, 0) {
/* Let's help with debugging */
if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
proj_log_debug(P, "Helmert parameters:");
- proj_log_debug(P, "x= % 3.5f y= % 3.5f z= % 3.5f", Q->xyz.x, Q->xyz.y, Q->xyz.z);
- proj_log_debug(P, "rx= % 3.5f ry= % 3.5f rz= % 3.5f",
+ proj_log_debug(P, "x= %8.5f y= %8.5f z= %8.5f", Q->xyz.x, Q->xyz.y, Q->xyz.z);
+ proj_log_debug(P, "rx= %8.5f ry= %8.5f rz= %8.5f",
Q->opk.o / ARCSEC_TO_RAD, Q->opk.p / ARCSEC_TO_RAD, Q->opk.k / ARCSEC_TO_RAD);
- proj_log_debug(P, "s=% 3.5f exact=% d transpose=% d",
- Q->scale, Q->exact, Q->transpose);
- proj_log_debug(P, "dx= % 3.5f dy= % 3.5f dz= % 3.5f", Q->dxyz.x, Q->dxyz.y, Q->dxyz.z);
- proj_log_debug(P, "drx=% 3.5f dry=% 3.5f drz=% 3.5f", Q->dopk.o, Q->dopk.p, Q->dopk.k);
- proj_log_debug(P, "ds=% 3.5f t_epoch=% 5.5f t_obs=% 5.5f", Q->dscale, Q->t_epoch, Q->t_obs);
+ proj_log_debug(P, "s= %8.5f exact=%d transpose=%d", Q->scale, Q->exact, Q->transpose);
+ proj_log_debug(P, "dx= %8.5f dy= %8.5f dz= %8.5f", Q->dxyz.x, Q->dxyz.y, Q->dxyz.z);
+ proj_log_debug(P, "drx=%8.5f dry=%8.5f drz=%8.5f", Q->dopk.o, Q->dopk.p, Q->dopk.k);
+ proj_log_debug(P, "ds= %8.5f t_epoch=%8.5f t_obs=%8.5f", Q->dscale, Q->t_epoch, Q->t_obs);
}
if ((Q->opk.o==0) && (Q->opk.p==0) && (Q->opk.k==0) && (Q->scale==0) &&
diff --git a/src/gie.c b/src/gie.c
index e0045ebd..b2673d13 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -387,6 +387,7 @@ Interpret <args> as a numeric followed by a linear decadal prefix.
Return the properly scaled numeric
******************************************************************************/
double s;
+ const double GRS80_DEG = 111319.4908; /* deg-to-m at equator of GRS80 */
const char *endp = args;
s = proj_strtod (args, (char **) &endp);
if (args==endp)
@@ -408,6 +409,10 @@ Return the properly scaled numeric
s /= 1e6;
else if (0==strcmp(endp, "nm"))
s /= 1e9;
+ else if (0==strcmp(endp, "rad"))
+ s = GRS80_DEG * proj_todeg (s);
+ else if (0==strcmp(endp, "deg"))
+ s = GRS80_DEG * s;
else
s *= default_scale;
return s;
@@ -553,19 +558,29 @@ using the "builtins" command verb.
}
-static PJ_COORD torad_coord (PJ_COORD a) {
- PJ_COORD c = a;
- c.lpz.lam = proj_torad (a.lpz.lam);
- c.lpz.phi = proj_torad (a.lpz.phi);
- return c;
+static PJ_COORD torad_coord (PJ *P, PJ_DIRECTION dir, PJ_COORD a) {
+ size_t i, n;
+ char *axis = "enut";
+ paralist *l = pj_param_exists (P->params, "axis");
+ if (l && dir==PJ_INV)
+ axis = l->param + strlen ("axis=");
+ for (i = 0, n = strlen (axis); i < n; i++)
+ if (strchr ("news", axis[i]))
+ a.v[i] = proj_torad (a.v[i]);
+ return a;
}
-static PJ_COORD todeg_coord (PJ_COORD a) {
- PJ_COORD c = a;
- c.lpz.lam = proj_todeg (a.lpz.lam);
- c.lpz.phi = proj_todeg (a.lpz.phi);
- return c;
+static PJ_COORD todeg_coord (PJ *P, PJ_DIRECTION dir, PJ_COORD a) {
+ size_t i, n;
+ char *axis = "enut";
+ paralist *l = pj_param_exists (P->params, "axis");
+ if (l && dir==PJ_FWD)
+ axis = l->param + strlen ("axis=");
+ for (i = 0, n = strlen (axis); i < n; i++)
+ if (strchr ("news", axis[i]))
+ a.v[i] = proj_todeg (a.v[i]);
+ return a;
}
@@ -627,7 +642,7 @@ back/forward transformation pairs.
d = d==HUGE_VAL? T.tolerance: d;
/* input ("accepted") values - probably in degrees */
- coo = proj_angular_input (T.P, T.dir)? torad_coord (T.a): T.a;
+ coo = proj_angular_input (T.P, T.dir)? torad_coord (T.P, T.dir, T.a): T.a;
r = proj_roundtrip (T.P, T.dir, ntrips, &coo);
if (r <= d)
@@ -747,7 +762,7 @@ Tell GIE what to expect, when transforming the ACCEPTed input
proj_errno_reset (T.P);
/* Try to carry out the operation - and expect failure */
- ci = proj_angular_input (T.P, T.dir)? torad_coord (T.a): T.a;
+ ci = proj_angular_input (T.P, T.dir)? torad_coord (T.P, T.dir, T.a): T.a;
co = expect_trans_n_dim (ci);
/* Failed to fail? - that's a failure */
@@ -780,20 +795,20 @@ Tell GIE what to expect, when transforming the ACCEPTed input
return expect_message_cannot_parse (args);
/* expected angular values, probably in degrees */
- ce = proj_angular_output (T.P, T.dir)? torad_coord (T.e): T.e;
+ ce = proj_angular_output (T.P, T.dir)? torad_coord (T.P, T.dir, T.e): T.e;
if (T.verbosity > 3)
printf ("EXPECTS %.4f %.4f %.4f %.4f\n", ce.v[0],ce.v[1],ce.v[2],ce.v[3]);
/* input ("accepted") values, also probably in degrees */
- ci = proj_angular_input (T.P, T.dir)? torad_coord (T.a): T.a;
+ ci = proj_angular_input (T.P, T.dir)? torad_coord (T.P, T.dir, T.a): T.a;
if (T.verbosity > 3)
printf ("ACCEPTS %.4f %.4f %.4f %.4f\n", ci.v[0],ci.v[1],ci.v[2],ci.v[3]);
/* angular output from proj_trans comes in radians */
co = expect_trans_n_dim (ci);
- T.b = proj_angular_output (T.P, T.dir)? todeg_coord (co): co;
+ T.b = proj_angular_output (T.P, T.dir)? todeg_coord (T.P, T.dir, co): co;
if (T.verbosity > 3)
- printf ("GOT %.4f %.4f %.4f %.4f\n", ci.v[0],ci.v[1],ci.v[2],ci.v[3]);
+ printf ("GOT %.4f %.4f %.4f %.4f\n", co.v[0],co.v[1],co.v[2],co.v[3]);
/* but there are a few more possible input conventions... */
if (proj_angular_output (T.P, T.dir)) {
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index a1025aac..2293cf2e 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -32,58 +32,137 @@
#include "proj_internal.h"
#include "projects.h"
-PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
+#define INPUT_UNITS P->left
+#define OUTPUT_UNITS P->right
+
+
+static PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
if (HUGE_VAL==coo.v[0])
return proj_coord_error ();
+ /* The helmert datum shift will choke unless it gets a sensible 4D coordinate */
+ if (HUGE_VAL==coo.v[2] && P->helmert) coo.v[2] = 0.0;
+ if (HUGE_VAL==coo.v[3] && P->helmert) coo.v[3] = 0.0;
+
/* Check validity of angular input coordinates */
- if (P->left==PJ_IO_UNITS_RADIANS) {
+ if (INPUT_UNITS==PJ_IO_UNITS_RADIANS) {
double t;
- LP lp = coo.lp;
/* check for latitude or longitude over-range */
- t = (lp.phi < 0 ? -lp.phi : lp.phi) - M_HALFPI;
- if (t > PJ_EPS_LAT || lp.lam > 10 || lp.lam < -10) {
+ t = (coo.lp.phi < 0 ? -coo.lp.phi : coo.lp.phi) - M_HALFPI;
+ if (t > PJ_EPS_LAT || coo.lp.lam > 10 || coo.lp.lam < -10) {
proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
return proj_coord_error ();
}
/* Clamp latitude to -90..90 degree range */
- if (lp.phi > M_HALFPI)
- lp.phi = M_HALFPI;
- if (lp.phi < -M_HALFPI)
- lp.phi = -M_HALFPI;
+ if (coo.lp.phi > M_HALFPI)
+ coo.lp.phi = M_HALFPI;
+ if (coo.lp.phi < -M_HALFPI)
+ coo.lp.phi = -M_HALFPI;
/* If input latitude is geocentrical, convert to geographical */
if (P->geoc)
coo = proj_geoc_lat (P, PJ_INV, coo);
+ /* Ensure longitude is in the -pi:pi range */
+ if (0==P->over)
+ coo.lp.lam = adjlon(coo.lp.lam);
+
+ if (P->hgridshift)
+ coo = proj_trans (P->hgridshift, PJ_INV, coo);
+ else if (P->helmert) {
+ coo = proj_trans (P->cart_wgs84, PJ_FWD, coo); /* Go cartesian in WGS84 frame */
+ coo = proj_trans (P->helmert, PJ_FWD, coo); /* Step into local frame */
+ coo = proj_trans (P->cart, PJ_INV, coo); /* Go back to angular using local ellps */
+ }
+ if (P->vgridshift)
+ coo = proj_trans (P->vgridshift, PJ_FWD, coo);
+
/* Distance from central meridian, taking system zero meridian into account */
- lp.lam = (lp.lam - P->from_greenwich) - P->lam0;
+ coo.lp.lam = (coo.lp.lam - P->from_greenwich) - P->lam0;
/* Ensure longitude is in the -pi:pi range */
if (0==P->over)
- lp.lam = adjlon(lp.lam);
+ coo.lp.lam = adjlon(coo.lp.lam);
- coo.lp = lp;
+ return coo;
}
+
+ /* We do not support gridshifts on cartesian input */
+ if (INPUT_UNITS==PJ_IO_UNITS_CARTESIAN && P->helmert)
+ return proj_trans (P->helmert, PJ_FWD, coo);
return coo;
}
-PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo) {
+
+static PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo) {
+
+ switch (OUTPUT_UNITS) {
+
+ /* Handle false eastings/northings and non-metric linear units */
+ case PJ_IO_UNITS_CARTESIAN:
+
+ if (P->is_geocent) {
+ coo = proj_trans (P->cart, PJ_FWD, coo);
+ }
+
+ coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0);
+ coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0);
+ coo.xyz.z = P->fr_meter * (coo.xyz.z + P->z0);
+ break;
/* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
- if (P->right==PJ_IO_UNITS_CLASSIC) {
+ case PJ_IO_UNITS_CLASSIC:
coo.xy.x *= P->a;
coo.xy.y *= P->a;
+
+ /* Falls through */ /* (<-- GCC warning silencer) */
+ /* to continue processing in common with PJ_IO_UNITS_PROJECTED */
+ case PJ_IO_UNITS_PROJECTED:
+ coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0);
+ coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0);
+ coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
+ break;
+
+ case PJ_IO_UNITS_WHATEVER:
+ break;
+
+ case PJ_IO_UNITS_RADIANS:
+ if (INPUT_UNITS==PJ_IO_UNITS_RADIANS)
+ break;
+
+ /* adjust longitude to central meridian */
+ if (0==P->over)
+ coo.lpz.lam = adjlon(coo.lpz.lam);
+
+ if (P->vgridshift)
+ coo = proj_trans (P->vgridshift, PJ_INV, coo);
+ if (P->hgridshift)
+ coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+ else if (P->helmert) {
+ coo = proj_trans (P->cart_wgs84, PJ_FWD, coo); /* Go cartesian in WGS84 frame */
+ coo = proj_trans (P->helmert, PJ_FWD, coo); /* Step into local frame */
+ coo = proj_trans (P->cart, PJ_INV, coo); /* Go back to angular using local ellps */
+ }
+
+ /* If input latitude was geocentrical, convert back to geocentrical */
+ if (P->geoc)
+ coo = proj_geoc_lat (P, PJ_FWD, coo);
+
+
+ /* Distance from central meridian, taking system zero meridian into account */
+ coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
+
+ /* adjust longitude to central meridian */
+ if (0==P->over)
+ coo.lpz.lam = adjlon(coo.lpz.lam);
}
- /* Handle false eastings/northings and non-metric linear units */
- coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0);
- coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0);
- coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
+ if (P->axisswap)
+ coo = proj_trans (P->axisswap, PJ_FWD, coo);
return coo;
}
diff --git a/src/pj_geocent.c b/src/pj_geocent.c
index eca62080..e676b735 100644
--- a/src/pj_geocent.c
+++ b/src/pj_geocent.c
@@ -54,6 +54,8 @@ PJ *PROJECTION(geocent) {
P->y0 = 0.0;
P->inv = inverse;
P->fwd = forward;
+ P->left = PJ_IO_UNITS_RADIANS;
+ P->right = PJ_IO_UNITS_CARTESIAN;
return P;
}
diff --git a/src/pj_init.c b/src/pj_init.c
index f25c0344..51a742c7 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -769,26 +769,3 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
proj_errno_restore (PIN, err);
return PIN;
}
-
-
-
-/************************************************************************/
-/* pj_free() */
-/* */
-/* This is the application callable entry point for destroying */
-/* a projection definition. It does work generic to all */
-/* projection types, and then calls the projection specific */
-/* free function, P->destructor(), to do local work. */
-/* In most cases P->destructor()==pj_default_destructor. */
-/************************************************************************/
-
-void pj_free(PJ *P) {
- if (0==P)
- return;
- /* free projection parameters - all the hard work is done by */
- /* pj_default_destructor (in pj_malloc.c), which is supposed */
- /* to be called as the last step of the local destructor */
- /* pointed to by P->destructor. In most cases, */
- /* pj_default_destructor actually *is* what is pointed to */
- P->destructor (P, 0);
-}
diff --git a/src/pj_inv.c b/src/pj_inv.c
index e1fb05f6..02a2ca53 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -32,46 +32,131 @@
#include "proj_internal.h"
#include "projects.h"
+#define INPUT_UNITS P->right
+#define OUTPUT_UNITS P->left
-
-PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
+static PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
if (coo.xyz.x == HUGE_VAL) {
proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
return proj_coord_error ();
}
+ /* The helmert datum shift will choke unless it gets a sensible 4D coordinate */
+ if (HUGE_VAL==coo.v[2] && P->helmert) coo.v[2] = 0.0;
+ if (HUGE_VAL==coo.v[3] && P->helmert) coo.v[3] = 0.0;
+
+ if (P->axisswap)
+ coo = proj_trans (P->axisswap, PJ_INV, coo);
+
+ /* Check validity of angular input coordinates */
+ if (INPUT_UNITS==PJ_IO_UNITS_RADIANS) {
+ double t;
+
+ /* check for latitude or longitude over-range */
+ t = (coo.lp.phi < 0 ? -coo.lp.phi : coo.lp.phi) - M_HALFPI;
+ if (t > PJ_EPS_LAT || coo.lp.lam > 10 || coo.lp.lam < -10) {
+ proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
+ return proj_coord_error ();
+ }
+
+ /* Clamp latitude to -90..90 degree range */
+ if (coo.lp.phi > M_HALFPI)
+ coo.lp.phi = M_HALFPI;
+ if (coo.lp.phi < -M_HALFPI)
+ coo.lp.phi = -M_HALFPI;
+
+ /* If input latitude is geocentrical, convert to geographical */
+ if (P->geoc)
+ coo = proj_geoc_lat (P, PJ_INV, coo);
+
+ /* Distance from central meridian, taking system zero meridian into account */
+ coo.lp.lam = (coo.lp.lam + P->from_greenwich) - P->lam0;
+
+ /* Ensure longitude is in the -pi:pi range */
+ if (0==P->over)
+ coo.lp.lam = adjlon(coo.lp.lam);
+
+ if (P->hgridshift)
+ coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+ else if (P->helmert) {
+ coo = proj_trans (P->cart, PJ_FWD, coo); /* Go cartesian in local frame */
+ coo = proj_trans (P->helmert, PJ_INV, coo); /* Step into WGS84 */
+ coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
+ }
+ if (P->vgridshift)
+ coo = proj_trans (P->vgridshift, PJ_INV, coo);
+ return coo;
+ }
+
+ /* Handle remaining possible input types */
+ switch (INPUT_UNITS) {
+ case PJ_IO_UNITS_WHATEVER:
+ return coo;
+
/* de-scale and de-offset */
- coo.xyz.x = (coo.xyz.x * P->to_meter - P->x0);
- coo.xyz.y = (coo.xyz.y * P->to_meter - P->y0);
- coo.xyz.z = (coo.xyz.z * P->vto_meter - P->z0);
-
- /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis */
- /* Multiplying by ra, rather than dividing by a because the CALCOFI projection */
- /* stomps on a and hence (apparently) depends on this to roundtrip correctly */
- /* (CALCOFI avoids further scaling by stomping - a better solution must be possible) */
- if (P->right==PJ_IO_UNITS_CLASSIC) {
+ case PJ_IO_UNITS_CARTESIAN:
+ coo.xyz.x = P->to_meter * coo.xyz.x - P->x0;
+ coo.xyz.y = P->to_meter * coo.xyz.y - P->y0;
+ coo.xyz.z = P->to_meter * coo.xyz.z - P->z0;
+
+ if (P->is_geocent) {
+ coo = proj_trans (P->cart, PJ_INV, coo);
+ }
+
+ return coo;
+
+ case PJ_IO_UNITS_PROJECTED:
+ case PJ_IO_UNITS_CLASSIC:
+ coo.xyz.x = P->to_meter * coo.xyz.x - P->x0;
+ coo.xyz.y = P->to_meter * coo.xyz.y - P->y0;
+ coo.xyz.z = P->vto_meter * coo.xyz.z - P->z0;
+ if (INPUT_UNITS==PJ_IO_UNITS_PROJECTED)
+ return coo;
+
+ /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis */
+ /* Multiplying by ra, rather than dividing by a because the CALCOFI projection */
+ /* stomps on a and hence (apparently) depends on this to roundtrip correctly */
+ /* (CALCOFI avoids further scaling by stomping - but a better solution is possible) */
coo.xyz.x *= P->ra;
coo.xyz.y *= P->ra;
+ return coo;
+ /* Silence some compiler warnings about PJ_IO_UNITS_RADIANS not handled */
+ default:
+ break;
}
+ /* Should not happen, so we could return pj_coord_err here */
return coo;
}
-PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo) {
+static PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo) {
if (coo.xyz.x == HUGE_VAL) {
proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
return proj_coord_error ();
}
- if (P->left==PJ_IO_UNITS_RADIANS) {
- /* Distance from central meridian, taking system zero meridian into account */
- coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
-
- /* adjust longitude to central meridian */
- if (0==P->over)
- coo.lpz.lam = adjlon(coo.lpz.lam);
+ if (OUTPUT_UNITS==PJ_IO_UNITS_RADIANS) {
+
+ if (INPUT_UNITS!=PJ_IO_UNITS_RADIANS) {
+ /* Distance from central meridian, taking system zero meridian into account */
+ coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
+
+ /* adjust longitude to central meridian */
+ if (0==P->over)
+ coo.lpz.lam = adjlon(coo.lpz.lam);
+
+ if (P->vgridshift)
+ coo = proj_trans (P->vgridshift, PJ_INV, coo);
+ if (P->hgridshift)
+ coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+ else if (P->helmert) {
+ coo = proj_trans (P->cart, PJ_FWD, coo); /* Go cartesian in local frame */
+ coo = proj_trans (P->helmert, PJ_INV, coo); /* Step into WGS84 */
+ coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
+ }
+ }
/* If input latitude was geocentrical, convert back to geocentrical */
if (P->geoc)
diff --git a/src/pj_malloc.c b/src/pj_malloc.c
index 63278a56..cfb7f658 100644
--- a/src/pj_malloc.c
+++ b/src/pj_malloc.c
@@ -162,6 +162,32 @@ void *pj_dealloc_params (PJ_CONTEXT *ctx, paralist *start, int errlev) {
}
+
+
+/************************************************************************/
+/* pj_free() */
+/* */
+/* This is the application callable entry point for destroying */
+/* a projection definition. It does work generic to all */
+/* projection types, and then calls the projection specific */
+/* free function, P->destructor(), to do local work. */
+/* In most cases P->destructor()==pj_default_destructor. */
+/************************************************************************/
+
+void pj_free(PJ *P) {
+ if (0==P)
+ return;
+ /* free projection parameters - all the hard work is done by */
+ /* pj_default_destructor, which is supposed */
+ /* to be called as the last step of the local destructor */
+ /* pointed to by P->destructor. In most cases, */
+ /* pj_default_destructor actually *is* what is pointed to */
+ P->destructor (P, 0);
+}
+
+
+
+
/*****************************************************************************/
void *pj_default_destructor (PJ *P, int errlev) { /* Destructor */
/*****************************************************************************
@@ -196,6 +222,14 @@ void *pj_default_destructor (PJ *P, int errlev) { /* Destructor */
/* free parameter list elements */
pj_dealloc_params (pj_get_ctx(P), P->params, errlev);
+ /* free the cs2cs emulation elements */
+ pj_free (P->axisswap);
+ pj_free (P->helmert);
+ pj_free (P->cart);
+ pj_free (P->cart_wgs84);
+ pj_free (P->hgridshift);
+ pj_free (P->vgridshift);
+
pj_dealloc (P->opaque);
return pj_dealloc(P);
}
diff --git a/src/proj.def b/src/proj.def
index 6ab0c007..099a78a5 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -155,3 +155,5 @@ EXPORTS
pj_approx_2D_trans @138
pj_approx_3D_trans @139
pj_has_inverse @140
+ pj_param_exists @141
+
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 2c6bac2d..fa24c430 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -380,7 +380,114 @@ char* proj_rtodms(char *s, double r, int pos, int neg) {
return rtodms(s, r, pos, neg);
}
+/*************************************************************************************/
+static PJ* skip_prep_fin(PJ *P) {
+/**************************************************************************************
+Skip prepare and finalize function for the various "helper operations" added to P when
+in cs2cs compatibility mode.
+**************************************************************************************/
+ P->skip_fwd_prepare = 1;
+ P->skip_fwd_finalize = 1;
+ P->skip_inv_prepare = 1;
+ P->skip_inv_finalize = 1;
+ return P;
+}
+
+/*************************************************************************************/
+static int pj_cs2cs_emulation_setup (PJ *P) {
+/**************************************************************************************
+If any cs2cs style modifiers are given (axis=..., towgs84=..., ) create the 4D API
+equivalent operations, so the preparation and finalization steps in the pj_inv/pj_fwd
+invocators can emulate the behaviour of pj_transform and the cs2cs app.
+**************************************************************************************/
+ PJ *Q;
+ paralist *p;
+ char def[1000];
+ if (0==P)
+ return 0;
+
+ /* Don't recurse when calling proj_create (which calls us back) */
+ if (pj_param_exists (P->params, "break_cs2cs_recursion"))
+ return 1;
+
+ /* Swap axes? */
+ p = pj_param_exists (P->params, "axis");
+
+ /* Don't axisswap if data are already in "enu" order */
+ if (p && (0!=strcmp ("enu", p->param))) {
+ sprintf (def, "break_cs2cs_recursion proj=axisswap axis=%s", P->axis);
+ Q = proj_create (P->ctx, def);
+ if (0==Q)
+ return 0;
+ P->axisswap = skip_prep_fin(Q);
+ }
+
+ /* Geoid grid(s) given? */
+ p = pj_param_exists (P->params, "geoidgrids");
+ if (p && strlen (p->param) > strlen ("geoidgrids=")) {
+ char *gridnames = p->param + strlen ("geoidgrids=");
+ sprintf (def, "break_cs2cs_recursion proj=vgridshift grids=%s", gridnames);
+ Q = proj_create (P->ctx, def);
+ if (0==Q)
+ return 0;
+ P->vgridshift = skip_prep_fin(Q);
+ }
+
+ /* Datum shift grid(s) given? */
+ p = pj_param_exists (P->params, "nadgrids");
+ if (p && strlen (p->param) > strlen ("nadgrids=")) {
+ char *gridnames = p->param + strlen ("nadgrids=");
+ sprintf (def, "break_cs2cs_recursion proj=hgridshift grids=%s", gridnames);
+ Q = proj_create (P->ctx, def);
+ if (0==Q)
+ return 0;
+ P->hgridshift = skip_prep_fin(Q);
+ }
+
+ /* We ignore helmert if we have grid shift */
+ p = P->hgridshift ? 0 : pj_param_exists (P->params, "towgs84");
+ while (p) {
+ char *s = p->param;
+ double *d = P->datum_params;
+ size_t n = strlen (s);
+
+ /* We ignore null helmert shifts (common in auto-translated resource files, e.g. epsg) */
+ if (0==d[0] && 0==d[1] && 0==d[2] && 0==d[3] && 0==d[4] && 0==d[5] && 0==d[6])
+ break;
+
+ if (n > 900)
+ return 0;
+ if (n <= 8) /* 8==strlen ("towgs84=") */
+ return 0;
+ sprintf (def, "break_cs2cs_recursion proj=helmert %s", s);
+ Q = proj_create (P->ctx, def);
+ if (0==Q)
+ return 0;
+ P->helmert = skip_prep_fin(Q);
+
+ break;
+ }
+
+ /* We also need cartesian/geographical transformations if we are working in */
+ /* geocentric/cartesian space or we need to do a Helmert transform. */
+ if (P->is_geocent || P->helmert) {
+ char *wgs84 = "ellps=WGS84";
+ sprintf (def, "break_cs2cs_recursion proj=cart");
+ Q = proj_create (P->ctx, def);
+ if (0==Q)
+ return 0;
+ pj_inherit_ellipsoid_def(P, Q);
+ P->cart = skip_prep_fin(Q);
+
+ sprintf (def, "break_cs2cs_recursion proj=cart %s", wgs84);
+ Q = proj_create (P->ctx, def);
+ if (0==Q)
+ return 0;
+ P->cart_wgs84 = skip_prep_fin(Q);
+ }
+ return 1;
+}
/*************************************************************************************/
@@ -394,9 +501,10 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
It may even use free formatting "proj = utm; zone =32 ellps= GRS80".
Note that the semicolon separator is allowed, but not required.
**************************************************************************************/
- PJ *P;
- char *args, **argv;
+ PJ *P;
+ char *args, **argv;
size_t argc, n;
+ int ret;
if (0==ctx)
ctx = pj_get_default_ctx ();
@@ -421,6 +529,12 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
pj_dealloc (argv);
pj_dealloc (args);
+
+ /* Support cs2cs-style modifiers */
+ ret = pj_cs2cs_emulation_setup (P);
+ if (0==ret)
+ return proj_destroy (P);
+
return P;
}
@@ -449,6 +563,7 @@ indicator, as in {"+proj=utm", "+zone=32"}, or leave it out, as in {"proj=utm",
return 0;
P = proj_create (ctx, c);
+
pj_dealloc ((char *) c);
return P;
}
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 41f95c07..2c6cbee2 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -84,10 +84,6 @@ void proj_context_inherit (PJ *parent, PJ *child);
PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P);
PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P);
-PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo);
-PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo);
-PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo);
-PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo);
PJ_COORD pj_approx_2D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
PJ_COORD pj_approx_3D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
diff --git a/src/projects.h b/src/projects.h
index 86ca2524..aeaf09ae 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -349,6 +349,14 @@ struct PJconsts {
enum pj_io_units left; /* Flags for input/output coordinate types */
enum pj_io_units right;
+ /* These PJs are used for implementing cs2cs style coordinate handling in the 4D API */
+ PJ *axisswap;
+ PJ *cart;
+ PJ *cart_wgs84;
+ PJ *helmert;
+ PJ *hgridshift;
+ PJ *vgridshift;
+
/*************************************************************************************
@@ -394,7 +402,7 @@ struct PJconsts {
double from_greenwich; /* prime meridian offset (in radians) */
double long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/
int is_long_wrap_set;
- char axis[4]; /* TODO: Description needed */
+ char axis[4]; /* Axis order, pj_transform/pj_adjust_axis */
/* New Datum Shift Grid Catalogs */
char *catalog_name;
diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie
new file mode 100644
index 00000000..853df993
--- /dev/null
+++ b/test/gie/4D-API_cs2cs-style.gie
@@ -0,0 +1,161 @@
+===============================================================================
+
+Test the 4D API handling of cs2cs style transformation options.
+
+These tests are mostly based on the same material as those in
+more_builtins.gie, since we are testing the same kinds of things,
+but provided through a different interface.
+
+===============================================================================
+
+
+<gie>
+
+-------------------------------------------------------------------------------
+Test the handling of the +towgs84 parameter
+-------------------------------------------------------------------------------
+This example is from Lotti Jivall: "Simplified transformations from
+ITRF2008/IGS08 to ETRS89 for maritime applications". The original
+XYZ data (cf. more_builtins.gie) have been transformed to LPZ using
+this command: echo x_val y_val z_val | cct -It0 proj=cart ellps=GRS80 --
+
+NOTE: Here, the ellipsoid has been swapped to WGS84, to align with
+ the WGS84 ellipsoid used in the cs2cs emulation introduced by
+ pj_cs2cs_emulation_setup()
+-------------------------------------------------------------------------------
+operation proj=latlong ellps=WGS84
+ towgs84 = 0.676780, 0.654950, -0.528270,
+ -0.022742, 0.012667, 0.022704, -0.010700
+-------------------------------------------------------------------------------
+tolerance 0.05 mm
+direction forward
+
+accept 13.4999969828397 54.9999995493206 -0.6034
+expect 13.4999906103972 54.9999957469562 -0.6374
+
+direction inverse
+
+accept 13.4999906103972 54.9999957469562 -0.6374
+expect 13.4999969828397 54.9999995493206 -0.6034
+-------------------------------------------------------------------------------
+
+
+
+-------------------------------------------------------------------------------
+This example is a random point, transformed from ED50 to ETRS89 using KMStrans2.
+
+NOTE: Signs swapped wrt. KMStrans2, which apparently uses frame rotation.
+-------------------------------------------------------------------------------
+operation proj=latlong ellps=intl
+ towgs84 = 081.07030, 089.36030, 115.75260,
+ 000.48488, 000.02436, 000.41321, 0.540645
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+# echo SWAPPED and trimmed - expect succcess
+# accept 16.819999997 55.170000002 61.0
+accept 16.82 55.17 61.0
+expect 16.8210462130 55.1705688946 29.0317
+-------------------------------------------------------------------------------
+This commented-out version of the example above was used to detect the sign
+discrepancy between KMSTrans2 and PROJ.
+-------------------------------------------------------------------------------
+operation proj=latlong ellps=intl
+ towgs84 = -081.07030, -089.36030, -115.75260,
+ -000.48488, -000.02436, -000.41321, -0.540645
+-------------------------------------------------------------------------------
+# echo NOT SWAPPED - Expect failure
+# accept 16.8210462130 55.1705688946 29.0317
+# expect 16.819999997 55.170000002 61.0
+-------------------------------------------------------------------------------
+
+
+
+-------------------------------------------------------------------------------
+operation proj=latlong nadgrids=nzgd2kgrid0005.gsb ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 1 nm
+accept 173 -45 0 0
+expect 172.999892181021551 -45.001620431954613 0 0
+direction inverse
+accept 172.999892181021551 -45.001620431954613 0 0
+expect 173 -45 0 0
+-------------------------------------------------------------------------------
+
+
+
+-------------------------------------------------------------------------------
+operation proj=latlong geoidgrids=egm96_15.gtx ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 15 cm # lax tolerance due to widespread bad egm96 file
+accept 12.5 55.5 0
+expect 12.5 55.5 -36.0213
+direction inverse
+accept 12.5 55.5 -36.0213
+expect 12.5 55.5 0
+-------------------------------------------------------------------------------
+operation proj=merc geoidgrids=egm96_15.gtx ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept 12.5 55.5 0
+expect 1391493.63492 7424275.19462 -36.0213
+direction inverse
+accept 1391493.63492 7424275.19462 -36.0213
+expect 12.5 55.5 0
+-------------------------------------------------------------------------------
+
+
+
+-------------------------------------------------------------------------------
+Same as the two above, but also do axis swapping.
+-------------------------------------------------------------------------------
+operation proj=latlong geoidgrids=egm96_15.gtx axis=neu ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 15 cm # lax tolerance due to widely distributed, bad egm96 file
+accept 12.5 55.5 0
+expect 55.5 12.5 -36.0213
+direction inverse
+accept 55.5 12.5 -36.0213
+expect 12.5 55.5 0
+-------------------------------------------------------------------------------
+operation proj=latlong geoidgrids=egm96_15.gtx axis=dne ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 15 cm # lax tolerance due to widely distributed, bad egm96 file
+accept 12.5 55.5 0
+expect 36.0213 55.5 12.5
+direction inverse
+accept 36.0213 55.5 12.5
+expect 12.5 55.5 0
+-------------------------------------------------------------------------------
+operation proj=merc geoidgrids=egm96_15.gtx ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept 12.5 55.5 0
+expect 1391493.63492 7424275.19462 -36.0213
+direction inverse
+accept 1391493.63492 7424275.19462 -36.0213
+expect 12.5 55.5 0
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+Some more complex axis swapping.
+-------------------------------------------------------------------------------
+operation proj=latlong geoidgrids=egm96_15.gtx axis=nue ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 15 cm # lax tolerance due to widely distributed, bad egm96 file
+accept 12.5 55.5 0
+expect 55.5 -36.0213 12.5
+direction inverse
+accept 55.5 -36.0213 12.5
+expect 12.5 55.5 0
+-------------------------------------------------------------------------------
+operation proj=merc geoidgrids=egm96_15.gtx axis=sue ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 15 cm
+accept 12.5 55.5 0
+expect -7424275.1946 -36.0213 1391493.6349 0.0000
+direction inverse
+accept -7424275.1946 -36.0213 1391493.6349 0.0000
+expect 12.5 55.5 0
+-------------------------------------------------------------------------------
+</gie>
diff --git a/test/gie/axisswap.gie b/test/gie/axisswap.gie
index ac148a92..190ca3e7 100644
--- a/test/gie/axisswap.gie
+++ b/test/gie/axisswap.gie
@@ -53,9 +53,14 @@ expect 3 -2 1 4
roundtrip 100
operation proj=axisswap axis=neu
-tolerance 0.000001 m
-accept 1 2 3 4
-expect 2 1 3 4
+tolerance 0 m
+accept 1 2 3
+expect 2 1 3
+
+operation proj=axisswap axis=nue
+tolerance 0 m
+accept 1 2 3
+expect 2 3 1
operation proj=axisswap axis=swd
tolerance 0.000001 m
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index 3877b58a..e8f1bfe1 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -1233,14 +1233,14 @@ Geocentric
operation +proj=geocent +ellps=GRS80 +lat_1=0.5 +lat_2=2
-------------------------------------------------------------------------------
tolerance 0.00010 mm
-accept 2 1
-expect 222638.981586547 111319.490793274
-accept 2 -1
-expect 222638.981586547 -111319.490793274
-accept -2 1
-expect -222638.981586547 111319.490793274
-accept -2 -1
-expect -222638.981586547 -111319.490793274
+accept 2 1 0
+expect 6373287.27950247 222560.09599219 110568.77482092
+accept 2 -1 0
+expect 6373287.27950247 222560.09599219 -110568.77482092
+accept -2 1 0
+expect 6373287.27950247 -222560.09599219 110568.77482092
+accept -2 -1 0
+expect 6373287.27950247 -222560.09599219 -110568.77482092
direction inverse
accept 200 100
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 20bffbd9..e21f0321 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -131,11 +131,11 @@ expect 12 55 0 0
-------------------------------------------------------------------------------
Finally test a pipeline with more than one init step
-------------------------------------------------------------------------------
-operation proj=pipeline step
- init=epsg:25832 inv step
- init=epsg:25833 step
- init=epsg:25833 inv step
- init=epsg:25832
+operation proj=pipeline
+ step init=epsg:25832 inv
+ step init=epsg:25833
+ step init=epsg:25833 inv
+ step init=epsg:25832
-------------------------------------------------------------------------------
accept 691875.63214 6098907.82501 0 0
expect 691875.63214 6098907.82501 0 0
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-29 07:25:46 +0000