Skip to content

Commit

Permalink
Plumbing for pipelines (#453)
Browse files Browse the repository at this point in the history 
* re-enter pipeline

The pipeline interface is now internally based on the pj_obs_api, which
simplifies the implementation significantly.

This is the first mock up - it compiles fine, but is currently untested

* pipeline code cleaned up

The pipeline code is now based on the PJ_OBS api (although you can still
invoke a pipeline through pj_fwd/pj_inv and their 3D brethren).

This has made it possible to eliminate scores of funky casts and
convoluted workarounds. The code is now way more straightforward and
mostly conforming with common C idioms..

Also, the proj.h / obs_api interface to the logging system has been
streamlined through the introduction of the pj_log_error, pj_log_debug,
and pj_log_trace functions.

* Geodesics + minor changes

First proj.h style interface to Charles Karney's geodesics code:
pj_lp_dist.
Also, renamed pj_apply -> pj_trans

* Extended Ellipsoidal Parameters

Second eccentricity, second and third flattening etc.

* Rename pj_debug_set -> pj_log_level

... and add self test code for PJ_pipeline

* Clean up missing pj_apply->pj_trans

* Clean up missing pj_obs_dist_2d rename

* pj_strerrno bug fixed. Some doc/comments added

(In response to a review by @kbevers)
  • Loading branch information
@busstoptaktik
busstoptaktik committed Nov 20, 2016
1 parent 3b1b439 commit fab2015
Show file tree
Hide file tree
Showing 16 changed files with 1,150 additions and 274 deletions.
160 changes: 136 additions & 24 deletions examples/pj_obs_api_test.c
View file
Expand Up @@ -17,7 +17,7 @@
pj_create (char *desc):
Create a new PJ object from a text description.
pj_apply (PJ *P, int direction, PJ_OBS obs):
pj_trans (PJ *P, int direction, PJ_OBS obs):
Forward or inverse transformation of obs.
pj_error (PJ *P):
Check error status of P.
Expand All @@ -42,19 +42,32 @@
*******************************************************************************/
#include <proj.h>

int pj_pipeline_selftest (void);



int main (void) {
PJ *p;
PJ *P;
PJ_OBS a, b;
char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"};
int err;
double dist;
XY cph_utm32;

/* Log everything libproj offers to log for you */
pj_debug_set (0, PJ_LOG_DEBUG_MINOR);
pj_log_level (0, PJ_LOG_TRACE);

/* An utm projection on the GRS80 ellipsoid */
p = pj_create ("+proj=utm +zone=32 +ellps=GRS80");
if (0==p)
P = pj_create ("+proj=utm +zone=32 +ellps=GRS80");
if (0==P)
return puts ("Oops"), 0;

/* Clean up */
pj_free (P);

/* Same projection, now using argc/argv style initialization */
P = pj_create_argv (3, args);
if (0==P)
return puts ("Oops"), 0;

/* zero initialize everything, then set (longitude, latitude) to (12, 55) */
Expand All @@ -64,37 +77,43 @@ int main (void) {
a.coo.lp.phi = TORAD(55);

/* Forward projection */
b = pj_apply (p, PJ_FWD, a);
b = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f\n", b.coo.enh.e, b.coo.enh.n);
cph_utm32 = b.coo.xy;

/* Inverse projection */
a = pj_apply (p, PJ_INV, b);
a = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi));

/* Null projection */
a = pj_apply (p, PJ_IDENT, a);
a = pj_trans (P, PJ_IDENT, a);
printf ("IDENT: %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi));

/* Forward again, to get two linear items for comparison */
a = pj_apply (p, PJ_FWD, a);
a = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f\n", b.coo.enh.e, b.coo.enh.n);

dist = pj_obs_dist_2d (a, b);
dist = pj_xy_dist (a.coo.xy, b.coo.xy);
printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9);

/* should be identical - checking whether null-init is done */
dist = pj_xyz_dist (a.coo.xyz, b.coo.xyz);
printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9);


/* Invalid projection */
a = pj_apply (p, 42, a);
a = pj_trans (P, 42, a);
if (a.coo.lpz.lam!=DBL_MAX)
printf ("%15.9f %15.9f\n", a.coo.lpz.lam, a.coo.lpz.phi);
err = pj_error (p);
err = pj_error (P);
printf ("pj_error: %d\n", err);

/* Clean up */
pj_free (p);
pj_free (P);

/* Now do some 3D transformations */
p = pj_create ("+proj=cart +ellps=GRS80");
if (0==p)
P = pj_create ("+proj=cart +ellps=GRS80");
if (0==P)
return puts ("Oops"), 0;

/* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
Expand All @@ -104,29 +123,122 @@ int main (void) {
a.coo.lpz.z = 100;

/* Forward projection: 3D-Cartesian-to-Ellipsoidal */
b = pj_apply (p, PJ_FWD, a);
b = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f %15.9f\n", b.coo.xyz.x, b.coo.xyz.y, b.coo.xyz.z);

/* Check roundtrip precision for 10000 iterations each way */
dist = pj_roundtrip (p, PJ_FWD, 10000, a);
dist = pj_roundtrip (P, PJ_FWD, 10000, a);
printf ("Roundtrip deviation, fwd (nm): %15.9f\n", dist*1e9*1e5);
dist = pj_roundtrip (p, PJ_INV, 10000, b);
dist = pj_roundtrip (P, PJ_INV, 10000, b);
printf ("Roundtrip deviation, inv (nm): %15.9f\n", dist*1e9);

/* Inverse projection: Ellipsoidal-to-3D-Cartesian */
b = pj_apply (p, PJ_INV, b);
b = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi), b.coo.lpz.z);

/* Move p to another context */
pj_context_renew (p);
b = pj_apply (p, PJ_FWD, b);
pj_context_renew (P);
b = pj_trans (P, PJ_FWD, b);
printf ("CTX1: %15.9f %15.9f %15.9f\n", b.coo.xyz.x, b.coo.xyz.y, b.coo.xyz.z);

/* Move it back to the default context */
pj_context_free (p);
b = pj_apply (p, PJ_INV, b);
pj_context_free (P);
b = pj_trans (P, PJ_INV, b);
printf ("CTX0: %15.9f %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi), b.coo.lpz.z);

pj_free (p);
pj_free (P);

/***************************************************************************
P I P E L I N E T E S T S
****************************************************************************
***************************************************************************/

/* forward-reverse geo->utm->geo */
P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm +step +proj=utm +inv");
if (0==P)
return puts ("Oops"), 0;
/* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
a = b = pj_obs_null;
a.coo.lpz.lam = TORAD(12);
a.coo.lpz.phi = TORAD(55);
a.coo.lpz.z = 100;
printf ("PRE: %15.9f %15.9f\n", a.coo.lpz.lam, a.coo.lpz.phi);

/* Forward projection */
b = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi));

/* Inverse projection */
a = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi));

pj_free (P);


/* And now the back-to-back situation utm->geo->utm */
P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm +inv +step +proj=utm");
if (0==P)
return puts ("Oops"), 0;

/* zero initialize everything, then set (easting, northing) to utm(12, 55) */
a = b = pj_obs_null;
a.coo.xy = cph_utm32;
printf ("PRE: %15.9f %15.9f\n", a.coo.xy.x, a.coo.xy.y);

/* Forward projection */
b = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f\n", b.coo.xy.x, b.coo.xy.y);

/* Inverse projection */
a = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f\n", a.coo.xy.x, a.coo.xy.y);

pj_free (P);



/* Finally testing a corner case: A rather pointless one-step pipeline geo->utm */
P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm");
if (0==P)
return puts ("Oops"), 0;

/* zero initialize everything, then set (easting, northing) to utm(12, 55) */
a = b = pj_obs_null;
a.coo.lpz.lam = TORAD(12);
a.coo.lpz.phi = TORAD(55);
printf ("PRE: %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi));
printf ("EXP: %15.9f %15.9f\n", cph_utm32.x, cph_utm32.y);

/* Forward projection */
b = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f\n", b.coo.xy.x, b.coo.xy.y);

/* Inverse projection */
a = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi));


/* Geodesic distance between two points with angular 2D coordinates */
a.coo.lp.lam = TORAD(12);
a.coo.lp.phi = TORAD(60);
b.coo.lp.lam = TORAD(12);
b.coo.lp.phi = TORAD(61);
dist = pj_lp_dist (P, a.coo.lp, b.coo.lp);
printf ("1 deg at 60N: %15.9f\n", dist);

a.coo.lp.lam = TORAD(12);
a.coo.lp.phi = TORAD(0.);
b.coo.lp.lam = TORAD(12);
b.coo.lp.phi = TORAD(1.);
dist = pj_lp_dist (P, a.coo.lp, b.coo.lp);
printf ("1 deg at Equator: %15.9f\n", dist);


pj_free (P);
pj_pipeline_selftest ();
return 0;
}
2 changes: 1 addition & 1 deletion src/Makefile.am
View file
Expand Up @@ -74,7 +74,7 @@ libproj_la_SOURCES = \
jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c \
pj_strtod.c \
\
pj_obs_api.c PJ_cart.c
pj_obs_api.c PJ_cart.c PJ_pipeline.c

install-exec-local:
rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
Expand Down
65 changes: 38 additions & 27 deletions src/PJ_cart.c
View file
Expand Up @@ -197,8 +197,9 @@ static LPZ geodetic (XYZ cartesian, PJ *P) {
c = cos(lpz.phi);
if (fabs(c) < 1e-6) {
/* poleward of 89.99994 deg, we avoid division by zero */
/* by computing the height as the cartesian z value */
/* minus the semiminor axis length */
/* by computing the height as the cartesian z value */
/* minus the semiminor axis length */
/* (potential improvement: approx. by 2nd order poly) */
lpz.z = cartesian.z - (cartesian.z > 0? b: -b);
}
else
Expand Down Expand Up @@ -247,21 +248,31 @@ int pj_cart_selftest (void) {return 0;}
#else
/* Testing quite a bit of the pj_obs_api as a side effect (inspired by pj_obs_api_test.c) */
int pj_cart_selftest (void) {
PJ *p;
PJ *P;
PJ_OBS a, b;
int err;
double dist;
char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"};


/* Log everything libproj offers to log for you */
pj_debug_set (0, PJ_LOG_DEBUG_MINOR);
pj_log_level (0, PJ_LOG_TRACE);

/* An utm projection on the GRS80 ellipsoid */
p = pj_create ("+proj=utm +zone=32 +ellps=GRS80");
if (0==p)
P = pj_create ("+proj=utm +zone=32 +ellps=GRS80");
if (0==P)
return 1;

/* Clean up */
pj_free (P);

/* Same projection, now using argc/argv style initialization */
P = pj_create_argv (3, args);
if (0==P)
return puts ("Oops"), 0;

/* Turn off logging */
pj_debug_set (0, PJ_LOG_NONE);
pj_log_level (0, PJ_LOG_NONE);

/* zero initialize everything, then set (longitude, latitude) to (12, 55) */
a = pj_obs_null;
Expand All @@ -270,38 +281,38 @@ int pj_cart_selftest (void) {
a.coo.lp.phi = TORAD(55);

/* Forward projection */
b = pj_apply (p, PJ_FWD, a);
b = pj_trans (P, PJ_FWD, a);

/* Inverse projection */
a = pj_apply (p, PJ_INV, b);
a = pj_trans (P, PJ_INV, b);

/* Null projection */
a = pj_apply (p, PJ_IDENT, a);
a = pj_trans (P, PJ_IDENT, a);

/* Forward again, to get two linear items for comparison */
a = pj_apply (p, PJ_FWD, a);
a = pj_trans (P, PJ_FWD, a);

dist = pj_obs_dist_2d (a, b);
dist = pj_xy_dist (a.coo.xy, b.coo.xy);
if (dist > 2e-9)
return 2;

/* Invalid projection */
a = pj_apply (p, 42, a);
a = pj_trans (P, 42, a);
if (a.coo.lpz.lam!=DBL_MAX)
return 3;
err = pj_error (p);
err = pj_error (P);
if (0==err)
return 4;

/* Clear error */
pj_error_set (p, 0);
pj_error_set (P, 0);

/* Clean up */
pj_free (p);
pj_free (P);

/* Now do some 3D transformations */
p = pj_create ("+proj=cart +ellps=GRS80");
if (0==p)
P = pj_create ("+proj=cart +ellps=GRS80");
if (0==P)
return 5;

/* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
Expand All @@ -311,26 +322,26 @@ int pj_cart_selftest (void) {
a.coo.lpz.z = 100;

/* Forward projection: 3D-Cartesian-to-Ellipsoidal */
b = pj_apply (p, PJ_FWD, a);
b = pj_trans (P, PJ_FWD, a);

/* Check roundtrip precision for 10000 iterations each way */
dist = pj_roundtrip (p, PJ_FWD, 10000, a);
dist = pj_roundtrip (p, PJ_INV, 10000, b);
dist = pj_roundtrip (P, PJ_FWD, 10000, a);
dist = pj_roundtrip (P, PJ_INV, 10000, b);
if (dist > 2e-9)
return 6;

/* Inverse projection: Ellipsoidal-to-3D-Cartesian */
b = pj_apply (p, PJ_INV, b);
b = pj_trans (P, PJ_INV, b);

/* Move p to another context */
pj_context_renew (p);
b = pj_apply (p, PJ_FWD, b);
pj_context_renew (P);
b = pj_trans (P, PJ_FWD, b);

/* Move it back to the default context */
pj_context_free (p);
b = pj_apply (p, PJ_INV, b);
pj_context_free (P);
b = pj_trans (P, PJ_INV, b);

pj_free (p);
pj_free (P);
return 0;
}
#endif

0 comments on commit fab2015

Please sign in to comment.