diff options
| -rw-r--r-- | src/PJ_cart.c | 10 | ||||
| -rw-r--r-- | src/PJ_hgridshift.c | 3 | ||||
| -rw-r--r-- | src/PJ_horner.c | 5 | ||||
| -rw-r--r-- | src/PJ_molodensky.c | 8 | ||||
| -rw-r--r-- | src/PJ_vgridshift.c | 7 | ||||
| -rw-r--r-- | src/gie.c | 53 | ||||
| -rw-r--r-- | src/proj.def | 47 | ||||
| -rw-r--r-- | src/proj.h | 5 | ||||
| -rw-r--r-- | src/proj_4D_api.c | 46 | ||||
| -rw-r--r-- | test/gie/axisswap.gie | 8 | ||||
| -rw-r--r-- | test/gie/deformation.gie | 30 |
11 files changed, 130 insertions, 92 deletions
diff --git a/src/PJ_cart.c b/src/PJ_cart.c index 52237f24..e2768c55 100644 --- a/src/PJ_cart.c +++ b/src/PJ_cart.c @@ -318,8 +318,8 @@ int pj_cart_selftest (void) { 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); + dist = proj_roundtrip (P, PJ_FWD, 10000, &a); + dist = proj_roundtrip (P, PJ_INV, 10000, &b); if (dist > 2e-9) return 7; @@ -331,7 +331,7 @@ int pj_cart_selftest (void) { a.lpz.z = 100; /* Forward projection: Ellipsoidal-to-3D-Cartesian */ - dist = proj_roundtrip (P, PJ_FWD, 1, a); + dist = proj_roundtrip (P, PJ_FWD, 1, &a); if (dist > 1e-12) return 8; @@ -340,12 +340,14 @@ int pj_cart_selftest (void) { 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); + 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); diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c index 26a2f471..d57af697 100644 --- a/src/PJ_hgridshift.c +++ b/src/PJ_hgridshift.c @@ -105,8 +105,9 @@ int pj_hgridshift_selftest (void) { a = proj_coord (0,0,0,0); a.lpz.lam = PJ_TORAD(173); a.lpz.phi = PJ_TORAD(-45); + b = a; - dist = proj_roundtrip (P, PJ_FWD, 1, a); + dist = proj_roundtrip (P, PJ_FWD, 1, &b); if (dist > 0.00000001) { printf("dist: %f\n",dist); return 1; diff --git a/src/PJ_horner.c b/src/PJ_horner.c index c28e3907..51c271ae 100644 --- a/src/PJ_horner.c +++ b/src/PJ_horner.c @@ -529,9 +529,10 @@ int pj_horner_selftest (void) { 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, a); + dist = proj_roundtrip (P, PJ_FWD, 1, &c); if (dist > 0.01) return 1; @@ -560,7 +561,7 @@ int pj_horner_selftest (void) { return 3; /* Check roundtrip precision for 1 iteration each way */ - dist = proj_roundtrip (P, PJ_FWD, 1, a); + dist = proj_roundtrip (P, PJ_FWD, 1, &a); if (dist > 0.01) return 4; diff --git a/src/PJ_molodensky.c b/src/PJ_molodensky.c index a35cabe4..73d0e5c2 100644 --- a/src/PJ_molodensky.c +++ b/src/PJ_molodensky.c @@ -318,7 +318,7 @@ int pj_molodensky_selftest (void) {return 0;} #else int pj_molodensky_selftest (void) { - PJ_COORD in, res, exp; + PJ_COORD in, ni, res, exp; PJ *P; /* Test the abridged Molodensky first. Example from appendix 3 of Deakin (2004). */ @@ -346,7 +346,8 @@ int pj_molodensky_selftest (void) { } /* let's try a roundtrip */ - if (proj_roundtrip(P, PJ_FWD, 100, in) > 1) { + ni = in; + if (proj_roundtrip(P, PJ_FWD, 100, &ni) > 1) { proj_destroy(P); return 12; } @@ -375,7 +376,8 @@ int pj_molodensky_selftest (void) { } /* let's try a roundtrip */ - if (proj_roundtrip(P, PJ_FWD, 100, in) > 1) { + ni = in; + if (proj_roundtrip(P, PJ_FWD, 100, &ni) > 1) { proj_destroy(P); return 22; } diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c index e3f3cbd3..5ab4a162 100644 --- a/src/PJ_vgridshift.c +++ b/src/PJ_vgridshift.c @@ -108,8 +108,8 @@ int pj_vgridshift_selftest (void) { a = proj_coord(0,0,0,0); a.lpz.lam = PJ_TORAD(12.5); a.lpz.phi = PJ_TORAD(55.5); - - dist = proj_roundtrip (P, PJ_FWD, 1, a); + b = a; + dist = proj_roundtrip (P, PJ_FWD, 1, &b); if (dist > 0.00000001) return 1; @@ -125,6 +125,9 @@ int pj_vgridshift_selftest (void) { if (proj_xyz_dist(expect.xyz, b.xyz) > 1e-4) failures++; expect.lpz.z = -35.880001068115234000; if (proj_xyz_dist(expect.xyz, b.xyz) > 1e-4) failures++; + /* manual roundtrip a->b, b<-b, b==a */ + b = proj_trans(P, PJ_INV, b); + if (proj_xyz_dist(a.xyz, b.xyz) > 1e-9) failures++; if (failures > 1) return 2; @@ -377,8 +377,6 @@ static double strtod_scaled (char *args, double default_scale) { static int banner (char *args) { char dots[] = {"..."}, nodots[] = {""}, *thedots = nodots; - if (T.total_ko > 0 && T.op_ko==0) - printf ("\n\n"); if (strlen(args) > 70) thedots = dots; fprintf (T.fout, "%s%-70.70s%s\n", delim, args, thedots); @@ -513,29 +511,36 @@ static int roundtrip (char *args) { int ntrips; double d, r, ans; char *endp; + PJ_COORD coo; + if (0==T.P) return another_failure (); + ans = proj_strtod (args, &endp); ntrips = (int) (endp==args? 100: fabs(ans)); d = strtod_scaled (endp, 1); d = d==HUGE_VAL? T.tolerance: d; - r = proj_roundtrip (T.P, PJ_FWD, ntrips, T.a); - if (r > d) { - if (T.verbosity > -1) { - if (0==T.op_ko && T.verbosity < 2) - banner (T.operation); - fprintf (T.fout, "%s", T.op_ko? " -----\n": delim); - fprintf (T.fout, " FAILURE in %s(%d):\n", opt_strip_path (T.curr_file), (int) lineno); - fprintf (T.fout, " roundtrip deviation: %.3f mm, expected: %.3f mm\n", 1000*r, 1000*d); - } - another_failure (); - } - else - another_success (); + coo = T.a; - return 0; + /* input ("accepted") values - probably in degrees */ + coo = proj_angular_input (T.P, T.dir)? torad_coord (T.a): T.a; + + r = proj_roundtrip (T.P, T.dir, ntrips, &coo); + if (r <= d) + return another_success (); + + if (T.verbosity > -1) { + if (0==T.op_ko && T.verbosity < 2) + banner (T.operation); + fprintf (T.fout, "%s", T.op_ko? " -----\n": delim); + fprintf (T.fout, " FAILURE in %s(%d):\n", opt_strip_path (T.curr_file), (int) lineno); + fprintf (T.fout, " roundtrip deviation: %.3f mm, expected: %.3f mm\n", 1000*r, 1000*d); + } + return another_failure (); } + + static int expect_message (double d, char *args) { another_failure (); @@ -589,7 +594,7 @@ static int expect (char *args) { ce = proj_angular_output (T.P, T.dir)? torad_coord (T.e): T.e; /* 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.a): T.a; /* angular output from proj_trans comes in radians */ co = proj_trans (T.P, T.dir, ci); @@ -598,21 +603,22 @@ static int expect (char *args) { /* but there are a few more possible input conventions... */ if (proj_angular_output (T.P, T.dir)) { double e = HUGE_VAL; - d = hypot (proj_lp_dist (T.P, ce.lp, co.lp), ce.lpz.z - co.lpz.z); + d = proj_lpz_dist (T.P, ce.lpz, co.lpz); /* check whether input was already in radians */ if (d > T.tolerance) - e = hypot (proj_lp_dist (T.P, T.e.lp, co.lp), T.e.lpz.z - co.lpz.z); + e = proj_lpz_dist (T.P, T.e.lpz, co.lpz); if (e < d) d = e; + /* or the tolerance may be based on euclidean distance */ if (d > T.tolerance) e = proj_xyz_dist (T.b.xyz, T.e.xyz); if (e < d) d = e; + } else d = proj_xyz_dist (T.b.xyz, T.e.xyz); - if (d > T.tolerance) return expect_message (d, args); @@ -663,11 +669,14 @@ fprintf (T.fout, "%s\n", args); static int dispatch (char *cmnd, char *args) { + int last_errno = proj_errno_reset (T.P); + if (0==level%2) { if (0==strcmp (cmnd, "BEGIN")) level++; return 0; } + if (0==strcmp (cmnd, "OPERATION")) return operation (args); if (0==strcmp (cmnd, "operation")) return operation (args); if (0==strcmp (cmnd, "ACCEPT")) return accept (args); @@ -688,6 +697,10 @@ static int dispatch (char *cmnd, char *args) { if (0==strcmp (cmnd, "echo")) return echo (args); if (0==strcmp (cmnd, "END")) return finish_previous_operation (args), level++, 0; if ('#'==cmnd[0]) return comment (args); + + if (proj_errno(T.P)) + printf ("#####***** ERRNO=%d\n", proj_errno(T.P)); + proj_errno_restore (T.P, last_errno); return 0; } diff --git a/src/proj.def b/src/proj.def index 8b1882ee..cc1793e1 100644 --- a/src/proj.def +++ b/src/proj.def @@ -117,32 +117,33 @@ EXPORTS proj_context_destroy @109 proj_lp_dist @110 - proj_xy_dist @111 - proj_xyz_dist @112 + proj_lpz_dist @111 + proj_xy_dist @112 + proj_xyz_dist @113 - proj_log_level @113 - proj_log_func @114 - proj_log_error @115 - proj_log_debug @116 - proj_log_trace @117 + proj_log_level @114 + proj_log_func @115 + proj_log_error @116 + proj_log_debug @117 + proj_log_trace @118 - proj_info @118 - proj_pj_info @119 - proj_grid_info @120 - proj_init_info @121 + proj_info @119 + proj_pj_info @120 + proj_grid_info @121 + proj_init_info @122 - proj_torad @122 - proj_todeg @123 - proj_rtodms @124 - proj_dmstor @125 + proj_torad @123 + proj_todeg @124 + proj_rtodms @125 + proj_dmstor @126 - proj_derivatives @126 - proj_factors @127 + proj_derivatives @127 + proj_factors @128 - proj_list_operations @128 - proj_list_ellps @129 - proj_list_units @130 - proj_list_prime_meridians @131 + proj_list_operations @129 + proj_list_ellps @130 + proj_list_units @131 + proj_list_prime_meridians @132 - proj_angular_input @132 - proj_angular_output @133 + proj_angular_input @133 + proj_angular_output @134 @@ -393,11 +393,14 @@ int proj_trans_array (PJ *P, PJ_DIRECTION direction, size_t n, PJ_COORD *coord); PJ_COORD proj_coord (double x, double y, double z, double t); /* Measure internal consistency - in forward or inverse direction */ -double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo); +double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo); /* Geodesic distance between two points with angular 2D coordinates */ double proj_lp_dist (const PJ *P, LP a, LP b); +/* The geodesic distance AND the vertical offset */ +double proj_lpz_dist (const PJ *P, LPZ a, LPZ b); + /* Euclidean distance between two points with linear 2D coordinates */ double proj_xy_dist (XY a, XY b); diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c index 66a272f2..ee4cb20f 100644 --- a/src/proj_4D_api.c +++ b/src/proj_4D_api.c @@ -76,6 +76,14 @@ double proj_lp_dist (const PJ *P, LP a, LP b) { return s12; } +/* The geodesic distance AND the vertical offset */ +double proj_lpz_dist (const PJ *P, LPZ a, LPZ b) { + PJ_COORD aa, bb; + aa.lpz = a; + bb.lpz = b; + return hypot (proj_lp_dist (P, aa.lp, bb.lp), a.z - b.z); +} + /* Euclidean distance between two points with linear 2D coordinates */ double proj_xy_dist (XY a, XY b) { return hypot (a.x - b.x, a.y - b.y); @@ -89,9 +97,9 @@ double proj_xyz_dist (XYZ a, XYZ b) { /* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */ -double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo) { +double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) { int i; - PJ_COORD o, u; + PJ_COORD o, u, org; if (0==P) return HUGE_VAL; @@ -101,34 +109,28 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo) { return HUGE_VAL; } - o = coo; + /* in the first half-step, we generate the output value */ + u = org = *coo; + o = *coo = proj_trans (P, direction, u); - switch (direction) { - case PJ_FWD: - for (i = 0; i < n; i++) { - u = pj_fwd4d (o, P); - o = pj_inv4d (u, P); - } - break; - case PJ_INV: - for (i = 0; i < n; i++) { - u = pj_inv4d (o, P); - o = pj_fwd4d (u, P); - } - break; - default: - proj_errno_set (P, EINVAL); - return HUGE_VAL; + /* now we take n-1 full steps */ + for (i = 0; i < n - 1; i++) { + u = proj_trans (P, -direction, o); + o = proj_trans (P, direction, u); } - /* checking for angular input since we do a roundtrip, and end where we begin */ + /* finally, we take the last half-step */ + u = proj_trans (P, -direction, o); + + /* checking for angular *input* since we do a roundtrip, and end where we begin */ if (proj_angular_input (P, direction)) - return hypot (proj_lp_dist (P, coo.lp, o.lp), coo.lpz.z - o.lpz.z); + return proj_lpz_dist (P, org.lpz, u.lpz); - return proj_xyz_dist (coo.xyz, coo.xyz); + return proj_xyz_dist (org.xyz, u.xyz); } + /* Apply the transformation P to the coordinate coo */ PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo) { if (0==P) diff --git a/test/gie/axisswap.gie b/test/gie/axisswap.gie index 2264a056..b9281ef7 100644 --- a/test/gie/axisswap.gie +++ b/test/gie/axisswap.gie @@ -63,4 +63,12 @@ TOLERANCE 0.00001 m ACCEPT 12 55 0 0 EXPECT -55 -12 0 0 +------------------------------------------------------------------------------- +operation +proj=aea +ellps=GRS80 +lat_1=0 +lat_2=2 +------------------------------------------------------------------------------- +tolerance 0.00010 mm +accept 2 1 +expect 222571.608757106 110653.326743030 +ROUNDTRIP 100 + END diff --git a/test/gie/deformation.gie b/test/gie/deformation.gie index 4173fa31..ce9aca75 100644 --- a/test/gie/deformation.gie +++ b/test/gie/deformation.gie @@ -1,4 +1,3 @@ -BEGIN =============================================================================== Test for the deformation operation - Kinematic Gridshifting @@ -10,26 +9,27 @@ The input coordinate is located at lon=60, lam=-160 - somewhere in Alaska. =============================================================================== +BEGIN ------------------------------------------------------------------------------- Test using only horizontal grid and +tobs parameter ------------------------------------------------------------------------------- OPERATION +proj=deformation +xy_grids=alaska +t_epoch=2016.0 +t_obs=2000.0 ------------------------------------------------------------------------------- -TOLERANCE 0.000001 m +TOLERANCE 1 um ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457 -EXPECT -3004295.5888766116 -1093474.1688513425 5500477.1338251457 -ROUNDTRIP 1000 +EXPECT -3004295.5888766116 -1093474.1688513425 5500477.1338251457 +# ROUNDTRIP 1000 ------------------------------------------------------------------------------- Test using only vertical grid and +tobs parameter ------------------------------------------------------------------------------- OPERATION +proj=deformation +z_grids=egm96_15.gtx +t_epoch=2016.0 +t_obs=2000.0 ------------------------------------------------------------------------------- -TOLERANCE 0.000001 m +TOLERANCE 1 um ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457 EXPECT -3004295.5882503074 -1093474.1690603832 5500234.008855661 -ROUNDTRIP 1000 +# ROUNDTRIP 1000 ------------------------------------------------------------------------------- Test using both horizontal and vertical grids as well as the +tobs parameter @@ -39,35 +39,37 @@ OPERATION +proj=deformation +xy_grids=alaska +z_grids=egm96_15.gtx +t_epoch=20 TOLERANCE 0.000001 m ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457 EXPECT -3004295.5888766116 -1093474.1688513425 5500234.008855661 -ROUNDTRIP 1000 +# ROUNDTRIP 1000 ------------------------------------------------------------------------------- Test using only horizontal grid ------------------------------------------------------------------------------- OPERATION +proj=deformation +xy_grids=alaska +t_epoch=2016.0 ------------------------------------------------------------------------------- -TOLERANCE 0.000001 m +TOLERANCE 1 um ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457 2000.0 EXPECT -3004295.5888766116 -1093474.1688513425 5500477.1338251457 2000.0 -ROUNDTRIP 1000 +# ROUNDTRIP 1000 ------------------------------------------------------------------------------- Test using only vertical grid ------------------------------------------------------------------------------- OPERATION +proj=deformation +z_grids=egm96_15.gtx +t_epoch=2016.0 ------------------------------------------------------------------------------- -TOLERANCE 0.000001 m +TOLERANCE 1 um ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457 2000.0 EXPECT -3004295.5882503074 -1093474.1690603832 5500234.008855661 2000.0 -ROUNDTRIP 1000 +# ROUNDTRIP 1000 ------------------------------------------------------------------------------- Test using both horizontal and vertical grids ------------------------------------------------------------------------------- -OPERATION proj=deformation xy_grids=alaska z_grids=egm96_15.gtx t_epoch=2016.0 +OPERATION +proj=deformation +xy_grids=alaska +z_grids=egm96_15.gtx +t_epoch=2016.0 +ellps=GRS80 ------------------------------------------------------------------------------- -TOLERANCE 0.000001 m +TOLERANCE 1 um ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457 2000.0 EXPECT -3004295.5888766116 -1093474.1688513425 5500234.008855661 2000.0 -ROUNDTRIP 1000 +# ROUNDTRIP 1000 + + END |