Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Plumbing for pipelines #453

Merged
merged 8 commits into from Nov 20, 2016
Merged

Plumbing for pipelines #453

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
d56a940
re-enter pipeline
Nov 14, 2016
75c78ce
pipeline code cleaned up
Nov 16, 2016
dcaf85b
Geodesics + minor changes
Nov 17, 2016
35df792
Extended Ellipsoidal Parameters
Nov 17, 2016
4711daf
Rename pj_debug_set -> pj_log_level
Nov 18, 2016
6249b47
Clean up missing pj_apply->pj_trans
Nov 18, 2016
8cdc347
Clean up missing pj_obs_dist_2d rename
Nov 18, 2016
f518cc6
pj_strerrno bug fixed. Some doc/comments added
Nov 19, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
160 changes: 136 additions & 24 deletions examples/pj_obs_api_test.c
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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