diff options
| author | Kristian Evers <kristianevers@gmail.com> | 2019-02-24 20:01:28 +0100 |
|---|---|---|
| committer | GitHub <noreply@github.com> | 2019-02-24 20:01:28 +0100 |
| commit | 450ea14050701c2487d7e600e668088808248794 (patch) | |
| tree | 74a2c16b0ab6539e8da9d3141cbb289e2ad89e7b | |
| parent | 8c607c8963b3879ac1e7473dc015b4d8194278f0 (diff) | |
| parent | 911852371a6a109445cd79e5176a89c539b2e768 (diff) | |
| download | PROJ-450ea14050701c2487d7e600e668088808248794.tar.gz PROJ-450ea14050701c2487d7e600e668088808248794.zip | |
Merge pull request #1292 from rouault/fix_robinson
Robinson: fix wrong values for forward path for latitudes >= 87.5 (fixes #1172), and fix inaccurace inverse method
| -rw-r--r-- | src/projections/robin.cpp | 10 | ||||
| -rw-r--r-- | test/gie/builtins.gie | 26 |
2 files changed, 29 insertions, 7 deletions
diff --git a/src/projections/robin.cpp b/src/projections/robin.cpp index 8f142aad..9f7908f6 100644 --- a/src/projections/robin.cpp +++ b/src/projections/robin.cpp @@ -7,7 +7,7 @@ PROJ_HEAD(robin, "Robinson") "\n\tPCyl, Sph"; #define V(C,z) (C.c0 + z * (C.c1 + z * (C.c2 + z * C.c3))) -#define DV(C,z) (C.c1 + z * (C.c2 + C.c2 + z * 3. * C.c3)) +#define DV(C,z) (C.c1 + 2 * z * C.c2 + z * z * 3. * C.c3) /* note: following terms based upon 5 deg. intervals in degrees. @@ -74,7 +74,7 @@ static const struct COEFS Y[] = { #define RC1 0.08726646259971647884 #define NODES 18 #define ONEEPS 1.000001 -#define EPS 1e-8 +#define EPS 1e-10 /* Not sure at all of the appropriate number for MAX_ITER... */ #define MAX_ITER 100 @@ -90,7 +90,7 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */ proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); return xy; } - if (i >= NODES) i = NODES - 1; + if (i >= NODES) i = NODES; dphi = RAD_TO_DEG * (dphi - RC1 * i); xy.x = V(X[i], dphi) * FXC * lp.lam; xy.y = V(Y[i], dphi) * FYC; @@ -133,10 +133,8 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ T = Y[i]; /* first guess, linear interp */ t = 5. * (lp.phi - T.c0)/(Y[i+1].c0 - T.c0); - /* make into root */ - T.c0 = (float)(T.c0 - lp.phi); for (iters = MAX_ITER; iters ; --iters) { /* Newton-Raphson */ - t -= t1 = V(T,t) / DV(T,t); + t -= t1 = (V(T,t) - lp.phi) / DV(T,t); if (fabs(t1) < EPS) break; } diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 37bdfc69..89970f07 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -4438,7 +4438,7 @@ Robinson =============================================================================== ------------------------------------------------------------------------------- -operation +proj=robin +a=6400000 +lat_1=0.5 +lat_2=2 +operation +proj=robin +a=6400000 ------------------------------------------------------------------------------- tolerance 0.1 mm accept 2 1 @@ -4450,6 +4450,18 @@ expect -189588.423282508 107318.530350703 accept -2 -1 expect -189588.423282508 -107318.530350703 +accept 0 89.5 +expect 0.000000000000 8639799.718722090125 + +accept 0 90 +expect 0.000000000000 8654720.000000000000 + +accept 0 -89.5 +expect 0.000000000000 -8639799.718722090125 + +accept 0 -90 +expect 0.000000000000 -8654720.000000000000 + direction inverse accept 200 100 expect 0.002109689 0.000931806 @@ -4460,6 +4472,18 @@ expect -0.002109689 0.000931806 accept -200 -100 expect -0.002109689 -0.000931806 +accept 0.000000000000 8639799.718722090125 +expect 0 89.5 + +accept 0.000000000000 8654720.000000000000 +expect 0 90 + +accept 0.000000000000 -8639799.718722090125 +expect 0 -89.5 + +accept 0.000000000000 -8654720.000000000000 +expect 0 -90 + =============================================================================== Roussilhe Stereographic |