aboutsummaryrefslogtreecommitdiff
path: root/src/gie.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/gie.c')
-rw-r--r--src/gie.c735
1 files changed, 729 insertions, 6 deletions
diff --git a/src/gie.c b/src/gie.c
index 09a33069..df4a0b07 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -435,8 +435,25 @@ static int operation (char *args) {
an operation is the general term describing something that can be
either a conversion or a transformation)
******************************************************************************/
+ int i, j, n;
T.op_id++;
+
+ /* compactify the args, so we can fit more info on a line in verbose mode */
+ n = (int) strlen (args);
+ for (i = j = 0; i < n; ) {
+ /* skip prefix whitespace */
+ while (isspace (args[i]))
+ i++;
+ /* move a whitespace delimited text string to the left, skipping over superfluous whitespace */
+ while ((0!=args[i]) && (!isspace (args[i])))
+ args[j++] = args[i++];
+ if (args[j+1]!=0)
+ args[j++] = ' ';
+ i++;
+ }
+ args[j++] = 0;
strcpy (&(T.operation[0]), args);
+
if (T.verbosity > 1) {
finish_previous_operation (args);
banner (args);
@@ -448,15 +465,67 @@ static int operation (char *args) {
direction ("forward");
tolerance ("0.5 mm");
+
if (T.P)
proj_destroy (T.P);
T.P = proj_create (0, args);
- if (0==T.P)
- errmsg(3, "Invalid operation definition!\n %s\n", args);
+
return 0;
}
+
+
+static int pj_unitconvert_selftest (void);
+static int pj_cart_selftest (void);
+static int pj_horner_selftest (void);
+/*****************************************************************************/
+static int builtins (char *args) {
+/*****************************************************************************
+ There are still a few tests that cannot be described using gie
+ primitives. Instead, they are implemented as builtins, and invoked
+ using the "builtins" command verb.
+******************************************************************************/
+ int i;
+ if (T.verbosity > 1) {
+ finish_previous_operation (args);
+ banner ("builtins: unitconvert, horner, cart");
+ }
+ T.op_ok = 0;
+ T.op_ko = 0;
+
+ i = pj_unitconvert_selftest ();
+ if (i!=0) {
+ printf ("pj_unitconvert_selftest fails with %d\n", i);
+ another_failure();
+ }
+ else
+ another_success ();
+
+
+ i = pj_cart_selftest ();
+ if (i!=0) {
+ printf ("pj_cart_selftest fails with %d\n", i);
+ another_failure();
+ }
+ else
+ another_success ();
+
+ i = pj_horner_selftest ();
+ if (i!=0) {
+ printf ("pj_horner_selftest fails with %d\n", i);
+ another_failure();
+ }
+ else
+ another_success ();
+
+ return 0;
+}
+
+
+
+
+
static PJ_COORD torad_coord (PJ_COORD a) {
PJ_COORD c = a;
c.lpz.lam = proj_torad (a.lpz.lam);
@@ -496,7 +565,7 @@ static int accept (char *args) {
******************************************************************************/
T.a = parse_coord (args);
if (T.verbosity > 3)
- printf ("# %s", args);
+ printf ("# %s\n", args);
return 0;
}
@@ -583,22 +652,58 @@ static int expect (char *args) {
******************************************************************************/
PJ_COORD ci, co, ce;
double d;
- if (0==T.P)
+ int expect_failure = 0;
+
+ if (0==strcmp (args, "failure"))
+ expect_failure = 1;
+
+ if (0==T.P && !expect_failure) {
+ errmsg(3, "Invalid operation definition!\n %s\n", args);
+ return another_failure ();
+ }
+
+
+ if (expect_failure) {
+ /* If we expect failure, and fail, then it's a success... */
+ if (0==T.P)
+ return another_success ();
+ /* We may still successfully fail even if the proj_create succeeded */
+ ci = proj_angular_input (T.P, T.dir)? torad_coord (T.a): T.a;
+ co = proj_trans (T.P, T.dir, ci);
+ if (co.xyz.x==HUGE_VAL)
+ return another_success ();
+ /* no - we didn't manage to purportedly fail */
return another_failure ();
+ }
+
+
+ if (T.verbosity > 3) {
+ puts (T.P->inverted? "INVERTED": "NOT INVERTED");
+ puts (T.dir== 1? "forward": "reverse");
+ puts (proj_angular_input (T.P, T.dir)? "angular in": "linear in");
+ puts (proj_angular_output (T.P, T.dir)? "angular out": "linear out");
+ }
T.e = parse_coord (args);
if (HUGE_VAL==T.e.v[0])
return expect_message_cannot_parse (args);
- /* expected angular values probably in degrees */
+ /* expected angular values, probably in degrees */
ce = proj_angular_output (T.P, T.dir)? torad_coord (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 */
+ /* input ("accepted") values, also probably in degrees */
ci = proj_angular_input (T.P, T.dir)? torad_coord (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 = proj_trans (T.P, T.dir, ci);
T.b = proj_angular_output (T.P, T.dir)? todeg_coord (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]);
/* but there are a few more possible input conventions... */
if (proj_angular_output (T.P, T.dir)) {
@@ -693,6 +798,8 @@ static int dispatch (char *cmnd, char *args) {
if (0==strcmp (cmnd, "direction")) return direction (args);
if (0==strcmp (cmnd, "TOLERANCE")) return tolerance (args);
if (0==strcmp (cmnd, "tolerance")) return tolerance (args);
+ if (0==strcmp (cmnd, "BUILTINS")) return builtins (args);
+ if (0==strcmp (cmnd, "builtins")) return builtins (args);
if (0==strcmp (cmnd, "ECHO")) return echo (args);
if (0==strcmp (cmnd, "echo")) return echo (args);
if (0==strcmp (cmnd, "END")) return finish_previous_operation (args), level++, 0;
@@ -806,3 +913,619 @@ static char *get_args (char *inp) {
return args;
return --args;
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+char tc32_utm32[] = {
+ " +proj=horner"
+ " +ellps=intl"
+ " +range=500000"
+ " +fwd_origin=877605.269066,6125810.306769"
+ " +inv_origin=877605.760036,6125811.281773"
+ " +deg=4"
+ " +fwd_v=6.1258112678e+06,9.9999971567e-01,1.5372750011e-10,5.9300860915e-15,2.2609497633e-19,4.3188227445e-05,2.8225130416e-10,7.8740007114e-16,-1.7453997279e-19,1.6877465415e-10,-1.1234649773e-14,-1.7042333358e-18,-7.9303467953e-15,-5.2906832535e-19,3.9984284847e-19"
+ " +fwd_u=8.7760574982e+05,9.9999752475e-01,2.8817299305e-10,5.5641310680e-15,-1.5544700949e-18,-4.1357045890e-05,4.2106213519e-11,2.8525551629e-14,-1.9107771273e-18,3.3615590093e-10,2.4380247154e-14,-2.0241230315e-18,1.2429019719e-15,5.3886155968e-19,-1.0167505000e-18"
+ " +inv_v=6.1258103208e+06,1.0000002826e+00,-1.5372762184e-10,-5.9304261011e-15,-2.2612705361e-19,-4.3188331419e-05,-2.8225549995e-10,-7.8529116371e-16,1.7476576773e-19,-1.6875687989e-10,1.1236475299e-14,1.7042518057e-18,7.9300735257e-15,5.2881862699e-19,-3.9990736798e-19"
+ " +inv_u=8.7760527928e+05,1.0000024735e+00,-2.8817540032e-10,-5.5627059451e-15,1.5543637570e-18,4.1357152105e-05,-4.2114813612e-11,-2.8523713454e-14,1.9109017837e-18,-3.3616407783e-10,-2.4382678126e-14,2.0245020199e-18,-1.2441377565e-15,-5.3885232238e-19,1.0167203661e-18"
+};
+
+
+char sb_utm32[] = {
+ " +proj=horner"
+ " +ellps=intl"
+ " +range=500000"
+ " +tolerance=0.0005"
+ " +fwd_origin=4.94690026817276e+05,6.13342113183056e+06"
+ " +inv_origin=6.19480258923588e+05,6.13258568148837e+06"
+ " +deg=3"
+ " +fwd_c=6.13258562111350e+06,6.19480105709997e+05,9.99378966275206e-01,-2.82153291753490e-02,-2.27089979140026e-10,-1.77019590701470e-09,1.08522286274070e-14,2.11430298751604e-15"
+ " +inv_c=6.13342118787027e+06,4.94690181709311e+05,9.99824464710368e-01,2.82279070814774e-02,7.66123542220864e-11,1.78425334628927e-09,-1.05584823306400e-14,-3.32554258683744e-15"
+};
+
+static int pj_horner_selftest (void) {
+ PJ *P;
+ PJ_COORD a, b, c;
+ double dist;
+
+ /* Real polynonia relating the technical coordinate system TC32 to "System 45 Bornholm" */
+ P = proj_create (PJ_DEFAULT_CTX, tc32_utm32);
+ if (0==P)
+ return 10;
+
+ a = b = proj_coord (0,0,0,0);
+ a.uv.v = 6125305.4245;
+ a.uv.u = 878354.8539;
+ c = a;
+
+ /* Check roundtrip precision for 1 iteration each way, starting in forward direction */
+ dist = proj_roundtrip (P, PJ_FWD, 1, &c);
+ if (dist > 0.01)
+ return 1;
+
+ /* The complex polynomial transformation between the "System Storebaelt" and utm32/ed50 */
+ P = proj_create (PJ_DEFAULT_CTX, sb_utm32);
+ if (0==P)
+ return 11;
+
+ /* Test value: utm32_ed50(620000, 6130000) = sb_ed50(495136.8544, 6130821.2945) */
+ a = b = c = proj_coord (0,0,0,0);
+ a.uv.v = 6130821.2945;
+ a.uv.u = 495136.8544;
+ c.uv.v = 6130000.0000;
+ c.uv.u = 620000.0000;
+
+ /* Forward projection */
+ b = proj_trans (P, PJ_FWD, a);
+ dist = proj_xy_dist (b.xy, c.xy);
+ if (dist > 0.001)
+ return 2;
+
+ /* Inverse projection */
+ b = proj_trans (P, PJ_INV, c);
+ dist = proj_xy_dist (b.xy, a.xy);
+ if (dist > 0.001)
+ return 3;
+
+ /* Check roundtrip precision for 1 iteration each way */
+ dist = proj_roundtrip (P, PJ_FWD, 1, &a);
+ if (dist > 0.01)
+ return 4;
+
+ proj_destroy(P);
+ return 0;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+/* Testing quite a bit of the pj_obs_api as a side effect (inspired by pj_obs_api_test.c) */
+static 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;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+static int test_time(char* args, double tol, double t_in, double t_exp) {
+ PJ_COORD in, out;
+ PJ *P = proj_create(PJ_DEFAULT_CTX, args);
+ int ret = 0;
+
+ if (P == 0)
+ return 5;
+
+ in.xyzt.t = t_in;
+
+ out = proj_trans(P, PJ_FWD, in);
+ if (fabs(out.xyzt.t - t_exp) > tol) {
+ proj_log_error(P, "out: %10.10g, expect: %10.10g", out.xyzt.t, t_exp);
+ ret = 1;
+ }
+ out = proj_trans(P, PJ_INV, out);
+ if (fabs(out.xyzt.t - t_in) > tol) {
+ proj_log_error(P, "out: %10.10g, expect: %10.10g", out.xyzt.t, t_in);
+ ret = 2;
+ }
+ pj_free(P);
+
+ proj_log_level(NULL, 0);
+ return ret;
+}
+
+static int test_xyz(char* args, double tol, PJ_TRIPLET in, PJ_TRIPLET exp) {
+ PJ_COORD out, obs_in;
+ PJ *P = proj_create(PJ_DEFAULT_CTX, args);
+ int ret = 0;
+
+ if (P == 0)
+ return 5;
+
+ obs_in.xyz = in.xyz;
+ out = proj_trans(P, PJ_FWD, obs_in);
+ if (proj_xyz_dist(out.xyz, exp.xyz) > tol) {
+ printf("exp: %10.10g, %10.10g, %10.10g\n", exp.xyz.x, exp.xyz.y, exp.xyz.z);
+ printf("out: %10.10g, %10.10g, %10.10g\n", out.xyz.x, out.xyz.y, out.xyz.z);
+ ret = 1;
+ }
+
+ out = proj_trans(P, PJ_INV, out);
+ if (proj_xyz_dist(out.xyz, in.xyz) > tol) {
+ printf("exp: %g, %g, %g\n", in.xyz.x, in.xyz.y, in.xyz.z);
+ printf("out: %g, %g, %g\n", out.xyz.x, out.xyz.y, out.xyz.z);
+ ret += 2;
+ }
+ proj_destroy(P);
+ proj_log_level(NULL, 0);
+ return ret;
+}
+
+
+static int pj_unitconvert_selftest (void) {
+ int ret = 0;
+ char args1[] = "+proj=unitconvert +t_in=decimalyear +t_out=decimalyear";
+ double in1 = 2004.25;
+
+ char args2[] = "+proj=unitconvert +t_in=gps_week +t_out=gps_week";
+ double in2 = 1782.0;
+
+ char args3[] = "+proj=unitconvert +t_in=mjd +t_out=mjd";
+ double in3 = 57390.0;
+
+ char args4[] = "+proj=unitconvert +t_in=gps_week +t_out=decimalyear";
+ double in4 = 1877.71428, exp4 = 2016.0;
+
+ char args5[] = "+proj=unitconvert +xy_in=m +xy_out=dm +z_in=cm +z_out=mm";
+ PJ_TRIPLET in5 = {{55.25, 23.23, 45.5}}, exp5 = {{552.5, 232.3, 455.0}};
+
+ char args6[] = "+proj=unitconvert +xy_in=m +xy_out=m +z_in=m +z_out=m";
+ PJ_TRIPLET in6 = {{12.3, 45.6, 7.89}};
+
+ ret = test_time(args1, 1e-6, in1, in1); if (ret) return ret + 10;
+ ret = test_time(args2, 1e-6, in2, in2); if (ret) return ret + 20;
+ ret = test_time(args3, 1e-6, in3, in3); if (ret) return ret + 30;
+ ret = test_time(args4, 1e-6, in4, exp4); if (ret) return ret + 40;
+ ret = test_xyz (args5, 1e-10, in5, exp5); if (ret) return ret + 50;
+ ret = test_xyz (args6, 1e-10, in6, in6); if (ret) return ret + 50;
+
+ return 0;
+
+}
+
+
+
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 11:39:21 +0000