aboutsummaryrefslogtreecommitdiff
path: root/src/PJ_cart.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/PJ_cart.c')
-rw-r--r--src/PJ_cart.c412
1 files changed, 0 insertions, 412 deletions
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index e2768c55..2bcf3b4b 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -219,415 +219,3 @@ PJ *CONVERSION(cart,1) {
return P;
}
-#ifndef PJ_SELFTEST
-/* selftest stub */
-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_CONTEXT *ctx;
- PJ *P;
- PJ_COORD a, b, obs[2];
- PJ_COORD coord[2];
-
- PJ_INFO info;
- PJ_PROJ_INFO pj_info;
- PJ_GRID_INFO grid_info;
- PJ_INIT_INFO init_info;
-
- PJ_DERIVS derivs;
- PJ_FACTORS factors;
-
- const PJ_OPERATIONS *oper_list;
- const PJ_ELLPS *ellps_list;
- const PJ_UNITS *unit_list;
- const PJ_PRIME_MERIDIANS *pm_list;
-
- int err;
- size_t n, sz;
- double dist, h, t;
- char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"};
- char *arg = {"+proj=utm +zone=32 +ellps=GRS80"};
- char buf[40];
-
- /* An utm projection on the GRS80 ellipsoid */
- P = proj_create (PJ_DEFAULT_CTX, arg);
- if (0==P)
- return 1;
-
-
- /* Clean up */
- proj_destroy (P);
-
- /* Same projection, now using argc/argv style initialization */
- P = proj_create_argv (PJ_DEFAULT_CTX, 3, args);
- if (0==P)
- return 2;
-
- /* zero initialize everything, then set (longitude, latitude) to (12, 55) */
- a = proj_coord (0,0,0,0);
- /* a.lp: The coordinate part of a, interpreted as a classic LP pair */
- a.lp.lam = PJ_TORAD(12);
- a.lp.phi = PJ_TORAD(55);
-
- /* Forward projection */
- b = proj_trans (P, PJ_FWD, a);
-
- /* Inverse projection */
- a = proj_trans (P, PJ_INV, b);
-
- /* Null projection */
- a = proj_trans (P, PJ_IDENT, a);
-
- /* Forward again, to get two linear items for comparison */
- a = proj_trans (P, PJ_FWD, a);
-
- dist = proj_xy_dist (a.xy, b.xy);
- if (dist > 2e-9)
- return 3;
-
- /* Clear any previous error */
- proj_errno_set (P, 0);
-
- /* Invalid projection */
- a = proj_trans (P, 42, a);
- if (a.lpz.lam!=HUGE_VAL)
- return 4;
- err = proj_errno (P);
- if (0==err)
- return 5;
-
- /* Clear error again */
- proj_errno_set (P, 0);
-
- /* Clean up */
- proj_destroy (P);
-
- /* Now do some 3D transformations */
- P = proj_create (PJ_DEFAULT_CTX, "+proj=cart +ellps=GRS80");
- if (0==P)
- return 6;
-
- /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
- a = b = proj_coord (0,0,0,0);
- a.lpz.lam = PJ_TORAD(12);
- a.lpz.phi = PJ_TORAD(55);
- a.lpz.z = 100;
-
- /* Forward projection: 3D-Cartesian-to-Ellipsoidal */
- b = proj_trans (P, PJ_FWD, a);
-
- /* Check roundtrip precision for 10000 iterations each way */
- dist = proj_roundtrip (P, PJ_FWD, 10000, &a);
- dist = proj_roundtrip (P, PJ_INV, 10000, &b);
- if (dist > 2e-9)
- return 7;
-
-
- /* Test at the North Pole */
- a = b = proj_coord (0,0,0,0);
- a.lpz.lam = PJ_TORAD(0);
- a.lpz.phi = PJ_TORAD(90);
- a.lpz.z = 100;
-
- /* Forward projection: Ellipsoidal-to-3D-Cartesian */
- dist = proj_roundtrip (P, PJ_FWD, 1, &a);
- if (dist > 1e-12)
- return 8;
-
- /* Test at the South Pole */
- a = b = proj_coord (0,0,0,0);
- a.lpz.lam = PJ_TORAD(0);
- a.lpz.phi = PJ_TORAD(-90);
- a.lpz.z = 100;
- b = a;
-
- /* Forward projection: Ellipsoidal-to-3D-Cartesian */
- dist = proj_roundtrip (P, PJ_FWD, 1, &a);
- if (dist > 1e-12)
- return 9;
-
-
- /* Inverse projection: 3D-Cartesian-to-Ellipsoidal */
- b = proj_trans (P, PJ_INV, b);
-
- /* Move p to another context */
- ctx = proj_context_create ();
- if (ctx==pj_get_default_ctx())
- return 10;
- proj_context_set (P, ctx);
- if (ctx != P->ctx)
- return 11;
- b = proj_trans (P, PJ_FWD, b);
-
- /* Move it back to the default context */
- proj_context_set (P, 0);
- if (pj_get_default_ctx() != P->ctx)
- return 12;
- proj_context_destroy (ctx);
-
- /* We go on with the work - now back on the default context */
- b = proj_trans (P, PJ_INV, b);
- proj_destroy (P);
-
-
- /* Testing proj_trans_generic () */
-
- /* An utm projection on the GRS80 ellipsoid */
- P = proj_create (PJ_DEFAULT_CTX, "+proj=utm +zone=32 +ellps=GRS80");
- if (0==P)
- return 13;
-
- obs[0] = proj_coord (PJ_TORAD(12), PJ_TORAD(55), 45, 0);
- obs[1] = proj_coord (PJ_TORAD(12), PJ_TORAD(56), 50, 0);
- sz = sizeof (PJ_COORD);
-
- /* Forward projection */
- a = proj_trans (P, PJ_FWD, obs[0]);
- b = proj_trans (P, PJ_FWD, obs[1]);
-
- n = proj_trans_generic (
- P, PJ_FWD,
- &(obs[0].lpz.lam), sz, 2,
- &(obs[0].lpz.phi), sz, 2,
- &(obs[0].lpz.z), sz, 2,
- 0, sz, 0
- );
- if (2!=n)
- return 14;
- if (a.lpz.lam != obs[0].lpz.lam) return 15;
- if (a.lpz.phi != obs[0].lpz.phi) return 16;
- if (a.lpz.z != obs[0].lpz.z) return 17;
- if (b.lpz.lam != obs[1].lpz.lam) return 18;
- if (b.lpz.phi != obs[1].lpz.phi) return 19;
- if (b.lpz.z != obs[1].lpz.z) return 20;
-
- /* now test the case of constant z */
- obs[0] = proj_coord (PJ_TORAD(12), PJ_TORAD(55), 45, 0);
- obs[1] = proj_coord (PJ_TORAD(12), PJ_TORAD(56), 50, 0);
- h = 27;
- t = 33;
- n = proj_trans_generic (
- P, PJ_FWD,
- &(obs[0].lpz.lam), sz, 2,
- &(obs[0].lpz.phi), sz, 2,
- &h, 0, 1,
- &t, 0, 1
- );
- if (2!=n)
- return 21;
- if (a.lpz.lam != obs[0].lpz.lam) return 22;
- if (a.lpz.phi != obs[0].lpz.phi) return 23;
- if (45 != obs[0].lpz.z) return 24;
- if (b.lpz.lam != obs[1].lpz.lam) return 25;
- if (b.lpz.phi != obs[1].lpz.phi) return 26;
- if (50 != obs[1].lpz.z) return 27; /* NOTE: unchanged */
- if (50==h) return 28;
-
- /* test proj_trans_array () */
-
- coord[0] = proj_coord (PJ_TORAD(12), PJ_TORAD(55), 45, 0);
- coord[1] = proj_coord (PJ_TORAD(12), PJ_TORAD(56), 50, 0);
- if (proj_trans_array (P, PJ_FWD, 2, coord))
- return 40;
-
- if (a.lpz.lam != coord[0].lpz.lam) return 41;
- if (a.lpz.phi != coord[0].lpz.phi) return 42;
- if (a.lpz.z != coord[0].lpz.z) return 43;
- if (b.lpz.lam != coord[1].lpz.lam) return 44;
- if (b.lpz.phi != coord[1].lpz.phi) return 45;
- if (b.lpz.z != coord[1].lpz.z) return 46;
-
- /* Clean up after proj_trans_* tests */
- proj_destroy (P);
-
- /* test proj_create_crs_to_crs() */
- P = proj_create_crs_to_crs(PJ_DEFAULT_CTX, "epsg:25832", "epsg:25833", NULL);
- if (P==0)
- return 50;
-
- a.xy.x = 700000.0;
- a.xy.y = 6000000.0;
- b.xy.x = 307788.8761171057;
- b.xy.y = 5999669.3036037628;
-
- a = proj_trans(P, PJ_FWD, a);
- if (dist > 1e-7)
- return 51;
- proj_destroy(P);
-
- /* let's make sure that only entries in init-files results in a usable PJ */
- P = proj_create_crs_to_crs(PJ_DEFAULT_CTX, "proj=utm +zone=32 +datum=WGS84", "proj=utm +zone=33 +datum=WGS84", NULL);
- if (P != 0) {
- proj_destroy(P);
- return 52;
- }
- proj_destroy(P);
-
- /* ********************************************************************** */
- /* Test info functions */
- /* ********************************************************************** */
-
- /* proj_info() */
- /* this one is difficult to test, since the output changes with the setup */
- info = proj_info();
- if (info.version[0] != '\0' ) {
- char tmpstr[64];
- sprintf(tmpstr, "%d.%d.%d", info.major, info.minor, info.patch);
- if (strcmp(info.version, tmpstr)) return 55;
- }
- if (info.release[0] == '\0') return 56;
- if (info.searchpath[0] == '\0') return 57;
-
- /* proj_pj_info() */
- P = proj_create(PJ_DEFAULT_CTX, "+proj=august"); /* august has no inverse */
- if (proj_pj_info(P).has_inverse) { proj_destroy(P); return 60; }
- proj_destroy(P);
-
- P = proj_create(PJ_DEFAULT_CTX, arg);
- pj_info = proj_pj_info(P);
- if ( !pj_info.has_inverse ) { proj_destroy(P); return 61; }
- if ( strcmp(pj_info.definition, arg) ) { proj_destroy(P); return 62; }
- if ( strcmp(pj_info.id, "utm") ) { proj_destroy(P); return 63; }
- proj_destroy(P);
-
- /* proj_grid_info() */
- grid_info = proj_grid_info("egm96_15.gtx");
- if ( strlen(grid_info.filename) == 0 ) return 64;
- if ( strcmp(grid_info.gridname, "egm96_15.gtx") ) return 65;
- grid_info = proj_grid_info("nonexistinggrid");
- if ( strlen(grid_info.filename) > 0 ) return 66;
-
- /* proj_init_info() */
- init_info = proj_init_info("unknowninit");
- if ( strlen(init_info.filename) != 0 ) return 67;
-
- init_info = proj_init_info("epsg");
- /* Need to allow for "Unknown" until all commonly distributed EPSG-files comes with a metadata section */
- if ( strcmp(init_info.origin, "EPSG") && strcmp(init_info.origin, "Unknown") ) return 69;
- if ( strcmp(init_info.name, "epsg") ) return 68;
-
-
-
- /* test proj_rtodms() and proj_dmstor() */
- if (strcmp("180dN", proj_rtodms(buf, M_PI, 'N', 'S')))
- return 70;
-
- if (proj_dmstor(&buf[0], NULL) != M_PI)
- return 71;
-
- if (strcmp("114d35'29.612\"S", proj_rtodms(buf, -2.0, 'N', 'S')))
- return 72;
-
- /* we can't expect perfect numerical accuracy so testing with a tolerance */
- if (fabs(-2.0 - proj_dmstor(&buf[0], NULL)) > 1e-7)
- return 73;
-
-
- /* test proj_derivatives_retrieve() and proj_factors_retrieve() */
- P = proj_create(PJ_DEFAULT_CTX, "+proj=merc");
- a = proj_coord (0,0,0,0);
- a.lp.lam = PJ_TORAD(12);
- a.lp.phi = PJ_TORAD(55);
-
- derivs = proj_derivatives(P, a.lp);
- if (proj_errno(P))
- return 80; /* derivs not created correctly */
-
- if ( fabs(derivs.x_l - 1.0) > 1e-5 ) return 81;
- if ( fabs(derivs.x_p - 0.0) > 1e-5 ) return 82;
- if ( fabs(derivs.y_l - 0.0) > 1e-5 ) return 83;
- if ( fabs(derivs.y_p - 1.73959) > 1e-5 ) return 84;
-
-
- factors = proj_factors(P, a.lp);
- if (proj_errno(P))
- return 85; /* factors not created correctly */
-
- /* check a few key characteristics of the Mercator projection */
- if (factors.omega != 0.0) return 86; /* angular distortion should be 0 */
- if (factors.thetap != M_PI_2) return 87; /* Meridian/parallel angle should be 90 deg */
- if (factors.conv != 0.0) return 88; /* meridian convergence should be 0 */
-
-
- proj_destroy(P);
-
- /* Check that proj_list_* functions work by looping through them */
- n = 0;
- for (oper_list = proj_list_operations(); oper_list->id; ++oper_list) n++;
- if (n == 0) return 90;
-
- n = 0;
- for (ellps_list = proj_list_ellps(); ellps_list->id; ++ellps_list) n++;
- if (n == 0) return 91;
-
- n = 0;
- for (unit_list = proj_list_units(); unit_list->id; ++unit_list) n++;
- if (n == 0) return 92;
-
- n = 0;
- for (pm_list = proj_list_prime_meridians(); pm_list->id; ++pm_list) n++;
- if (n == 0) return 93;
-
-
- /* check io-predicates */
-
- /* angular in on fwd, linear out */
- P = proj_create (PJ_DEFAULT_CTX, "+proj=cart +ellps=GRS80");
- if (0==P) return 0;
- if (!proj_angular_input (P, PJ_FWD)) return 100;
- if ( proj_angular_input (P, PJ_INV)) return 101;
- if ( proj_angular_output (P, PJ_FWD)) return 102;
- if (!proj_angular_output (P, PJ_INV)) return 103;
- P->inverted = 1;
- if ( proj_angular_input (P, PJ_FWD)) return 104;
- if (!proj_angular_input (P, PJ_INV)) return 105;
- if (!proj_angular_output (P, PJ_FWD)) return 106;
- if ( proj_angular_output (P, PJ_INV)) return 107;
- proj_destroy(P);
-
- /* angular in and out */
- P = proj_create(PJ_DEFAULT_CTX,
- "+proj=molodensky +a=6378160 +rf=298.25 "
- "+da=-23 +df=-8.120449e-8 +dx=-134 +dy=-48 +dz=149 "
- "+abridged "
- );
- if (0==P) return 0;
- if (!proj_angular_input (P, PJ_FWD)) return 108;
- if (!proj_angular_input (P, PJ_INV)) return 109;
- if (!proj_angular_output (P, PJ_FWD)) return 110;
- if (!proj_angular_output (P, PJ_INV)) return 111;
- P->inverted = 1;
- if (!proj_angular_input (P, PJ_FWD)) return 112;
- if (!proj_angular_input (P, PJ_INV)) return 113;
- if (!proj_angular_output (P, PJ_FWD)) return 114;
- if (!proj_angular_output (P, PJ_INV)) return 115;
- proj_destroy(P);
-
- /* linear in and out */
- P = proj_create(PJ_DEFAULT_CTX,
- " +proj=helmert +ellps=GRS80"
- " +x=0.0127 +y=0.0065 +z=-0.0209 +s=0.00195"
- " +rx=-0.00039 +ry=0.00080 +rz=-0.00114"
- " +dx=-0.0029 +dy=-0.0002 +dz=-0.0006 +ds=0.00001"
- " +drx=-0.00011 +dry=-0.00019 +drz=0.00007"
- " +epoch=1988.0 +transpose"
- );
- if (0==P) return 0;
- if (proj_angular_input (P, PJ_FWD)) return 116;
- if (proj_angular_input (P, PJ_INV)) return 117;
- if (proj_angular_output (P, PJ_FWD)) return 118;
- if (proj_angular_output (P, PJ_INV)) return 119;
- P->inverted = 1;
- if (proj_angular_input (P, PJ_FWD)) return 120;
- if (proj_angular_input (P, PJ_INV)) return 121;
- if (proj_angular_output (P, PJ_FWD)) return 122;
- if (proj_angular_output (P, PJ_INV)) return 123;
- proj_destroy(P);
-
-
- return 0;
-}
-
-
-#endif