Skip to content

Commit

Permalink
Merge 5c8b26d into 53b30f1
Browse files Browse the repository at this point in the history
  • Loading branch information
@busstoptaktik
busstoptaktik authored Jan 14, 2018
2 parents 53b30f1 + 5c8b26d commit 71f0b84
Show file tree
Hide file tree
Showing 13 changed files with 525 additions and 82 deletions.
1 change: 1 addition & 0 deletions src/PJ_axisswap.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,6 +239,7 @@ PJ *CONVERSION(axisswap,0) {
}
}
n = 3;
P->inverted = 1;
}

/* check for duplicate axes */
Expand Down
29 changes: 22 additions & 7 deletions src/PJ_helmert.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -486,6 +486,22 @@ 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 */
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;
Expand Down Expand Up @@ -570,14 +586,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) &&
Expand Down
76 changes: 52 additions & 24 deletions src/gie.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,8 +133,6 @@ typedef struct ffio {
size_t level;
} ffio;

FILE *test = 0;

static int get_inp (ffio *F);
static int skip_to_next_tag (ffio *F);
static int step_into_gie_block (ffio *F);
Expand Down Expand Up @@ -179,6 +177,7 @@ typedef struct {
PJ_COORD a, b, c, e;
PJ_DIRECTION dir;
int verbosity;
int gigs_infinity_norm;
int op_id;
int op_ok, op_ko;
int total_ok, total_ko;
Expand Down Expand Up @@ -231,7 +230,7 @@ static const char usage[] = {

int main (int argc, char **argv) {
int i;
const char *longflags[] = {"v=verbose", "q=quiet", "h=help", "l=list", 0};
const char *longflags[] = {"v=verbose", "q=quiet", "h=help", "l=list", "g=gigs", 0};
const char *longkeys[] = {"o=output", 0};
OPTARGS *o;

Expand All @@ -240,8 +239,7 @@ int main (int argc, char **argv) {
T.verbosity = 1;
T.tolerance = 5e-4;

test = fopen ("test.c", "wt");
o = opt_parse (argc, argv, "hlvq", "o", longflags, longkeys);
o = opt_parse (argc, argv, "ghlvq", "o", longflags, longkeys);
if (0==o)
return 0;

Expand All @@ -253,6 +251,7 @@ int main (int argc, char **argv) {
if (opt_given (o, "l"))
return list_err_codes ();

T.gigs_infinity_norm = opt_given (o, "g");

T.verbosity = opt_given (o, "q");
if (T.verbosity)
Expand Down Expand Up @@ -553,19 +552,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_FWD)
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_INV)
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;
}


Expand Down Expand Up @@ -626,8 +635,13 @@ back/forward transformation pairs.
d = strtod_scaled (endp, 1);
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;

/* We approximate the infinity norm by L2, and turn tolerance d from degree to meter */
if (T.gigs_infinity_norm)
d *= 111000;

r = proj_roundtrip (T.P, T.dir, ntrips, &coo);
if (r <= d)
Expand Down Expand Up @@ -747,7 +761,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 */
Expand Down Expand Up @@ -780,21 +794,37 @@ 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]);

/* not sure whether we should do the check in degrees or radians... */
if (T.gigs_infinity_norm) {
int ic;
double dd;
d = 0;
for (ic = 0; ic < 4; ic++) {
dd = fabs (co.v[ic] - ce.v[ic]);
if (dd > d)
d = dd;
}
if (d > T.tolerance)
return expect_message (d, args);
return another_success ();
}


/* but there are a few more possible input conventions... */
if (proj_angular_output (T.P, T.dir)) {
double e = HUGE_VAL;
Expand All @@ -813,12 +843,10 @@ Tell GIE what to expect, when transforming the ACCEPTed input
}
else
d = proj_xyz_dist (T.b.xyz, T.e.xyz);

if (d > T.tolerance)
return expect_message (d, args);

another_success ();

return 0;
return another_success ();
}


Expand Down
74 changes: 65 additions & 9 deletions src/pj_fwd.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,12 +32,19 @@
#include "proj_internal.h"
#include "projects.h"

#define INPUT_UNITS P->left
#define OUTPUT_UNITS P->right


PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
if (HUGE_VAL==coo.v[0])
return proj_coord_error ();

if (P->axisswap)
coo = proj_trans (P->axisswap, 1, coo);

/* 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;

Expand All @@ -54,36 +61,85 @@ PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
if (lp.phi < -M_HALFPI)
lp.phi = -M_HALFPI;


coo.lp = lp;

/* 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 */
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 = lp;
coo.lp.lam = adjlon(coo.lp.lam);

if (P->hgridshift)
coo = proj_trans (P->hgridshift, -1, coo);
else if (P->helmert) {
coo = proj_trans (P->cart, 1, coo);
coo = proj_trans (P->helmert, 1, coo);
coo = proj_trans (P->cart, -1, coo);
}
if (P->vgridshift)
coo = proj_trans (P->vgridshift, 1, coo);
return coo;
}


/* We do not support gridshifts on cartesian input */
if (INPUT_UNITS==PJ_IO_UNITS_CARTESIAN && P->helmert)
return proj_trans (P->helmert, 1, coo);
return coo;
}



PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo) {

/* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
if (P->right==PJ_IO_UNITS_CLASSIC) {
if (OUTPUT_UNITS==PJ_IO_UNITS_CLASSIC) {
coo.xy.x *= P->a;
coo.xy.y *= P->a;
}

/* 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 (OUTPUT_UNITS==PJ_IO_UNITS_CARTESIAN) {
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);
}
if (OUTPUT_UNITS==PJ_IO_UNITS_PROJECTED || OUTPUT_UNITS==PJ_IO_UNITS_CLASSIC) {
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 (OUTPUT_UNITS==PJ_IO_UNITS_RADIANS && 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, -1, coo);
if (P->hgridshift)
coo = proj_trans (P->hgridshift, 1, coo);
else if (P->helmert) {
coo = proj_trans (P->cart, 1, coo);
coo = proj_trans (P->helmert, -1, coo);
coo = proj_trans (P->cart, -1, coo);
}

/* If input latitude was geocentrical, convert back to geocentrical */
if (P->geoc)
coo = proj_geoc_lat (P, PJ_FWD, coo);
}


return coo;
}
Expand Down
23 changes: 0 additions & 23 deletions src/pj_init.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -768,26 +768,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);
}
Loading

0 comments on commit 71f0b84

Please sign in to comment.