Skip to content

Commit

Permalink
Merge pull request #527 from kbevers/issue-526
Browse files Browse the repository at this point in the history 
Correct calculation of meridian convergence for non-conformal projections
  • Loading branch information
@kbevers
kbevers committed Jul 5, 2017
2 parents ee3ec9d + 1e8e527 commit ced55e8
Show file tree
Hide file tree
Showing 2 changed files with 138 additions and 107 deletions.
83 changes: 54 additions & 29 deletions src/pj_deriv.c
View file
Original file line number Diff line number Diff line change
@@ -1,33 +1,58 @@
/* dervative of (*P->fwd) projection */
#define PJ_LIB__
#include "projects.h"
int
pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {
XY t;

lp.lam += h;
lp.phi += h;
if (fabs(lp.phi) > M_HALFPI) return 1;
h += h;
t = (*P->fwd)(lp, P);
if (t.x == HUGE_VAL) return 1;
der->x_l = t.x; der->y_p = t.y; der->x_p = -t.x; der->y_l = -t.y;
lp.phi -= h;
if (fabs(lp.phi) > M_HALFPI) return 1;
t = (*P->fwd)(lp, P);
if (t.x == HUGE_VAL) return 1;
der->x_l += t.x; der->y_p -= t.y; der->x_p += t.x; der->y_l -= t.y;
lp.lam -= h;
t = (*P->fwd)(lp, P);
if (t.x == HUGE_VAL) return 1;
der->x_l -= t.x; der->y_p -= t.y; der->x_p += t.x; der->y_l += t.y;
lp.phi += h;
t = (*P->fwd)(lp, P);
if (t.x == HUGE_VAL) return 1;
der->x_l -= t.x; der->y_p += t.y; der->x_p -= t.x; der->y_l += t.y;
der->x_l /= (h += h);
der->y_p /= h;
der->x_p /= h;
der->y_l /= h;
return 0;

int pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {
XY t;

lp.lam += h;
lp.phi += h;
if (fabs(lp.phi) > M_HALFPI)
return 1;

h += h;
t = (*P->fwd)(lp, P);
if (t.x == HUGE_VAL)
return 1;

der->x_l = t.x;
der->y_p = t.y;
der->x_p = t.x;
der->y_l = t.y;
lp.phi -= h;
if (fabs(lp.phi) > M_HALFPI)
return 1;

t = (*P->fwd)(lp, P);
if (t.x == HUGE_VAL)
return 1;

der->x_l += t.x;
der->y_p -= t.y;
der->x_p -= t.x;
der->y_l += t.y;
lp.lam -= h;
t = (*P->fwd)(lp, P);
if (t.x == HUGE_VAL)
return 1;

der->x_l -= t.x;
der->y_p -= t.y;
der->x_p -= t.x;
der->y_l -= t.y;
lp.phi += h;
t = (*P->fwd)(lp, P);
if (t.x == HUGE_VAL)
return 1;

der->x_l -= t.x;
der->y_p += t.y;
der->x_p += t.x;
der->y_l -= t.y;
der->x_l /= (h += h);
der->y_p /= h;
der->x_p /= h;
der->y_l /= h;

return 0;
}
162 changes: 84 additions & 78 deletions src/pj_factors.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,90 +2,96 @@
#define PJ_LIB__
#include <projects.h>
#include <errno.h>

#ifndef DEFAULT_H
#define DEFAULT_H 1e-5 /* radian default for numeric h */
#endif

#define EPS 1.0e-12
int
pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
struct DERIVS der;
double cosphi, t, n, r;

der.x_l = 0.0;
der.x_p = 0.0;
der.y_l = 0.0;
der.y_p = 0.0;
int pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
struct DERIVS der;
double cosphi, t, n, r;

der.x_l = 0.0;
der.x_p = 0.0;
der.y_l = 0.0;
der.y_p = 0.0;

/* check for forward and latitude or longitude overange */
if ((t = fabs(lp.phi)-M_HALFPI) > EPS || fabs(lp.lam) > 10.) {
pj_ctx_set_errno( P->ctx, -14);
return 1;
} else { /* proceed */
errno = pj_errno = 0;
/* check for forward and latitude or longitude overange */
if ((t = fabs(lp.phi)-M_HALFPI) > EPS || fabs(lp.lam) > 10.) {
pj_ctx_set_errno( P->ctx, -14);
return 1;
} else { /* proceed */
errno = pj_errno = 0;
P->ctx->last_errno = 0;

if (h < EPS)
h = DEFAULT_H;
if (fabs(lp.phi) > (M_HALFPI - h))
/* adjust to value around pi/2 where derived still exists*/
lp.phi = lp.phi < 0. ? (-M_HALFPI+h) : (M_HALFPI-h);
else if (P->geoc)
lp.phi = atan(P->rone_es * tan(lp.phi));
lp.lam -= P->lam0; /* compute del lp.lam */
if (!P->over)
lp.lam = adjlon(lp.lam); /* adjust del longitude */
if (P->spc) /* get what projection analytic values */
P->spc(lp, P, fac);
if (((fac->code & (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) !=
(IS_ANAL_XL_YL+IS_ANAL_XP_YP)) &&
pj_deriv(lp, h, P, &der))
return 1;
if (!(fac->code & IS_ANAL_XL_YL)) {
fac->der.x_l = der.x_l;
fac->der.y_l = der.y_l;
}
if (!(fac->code & IS_ANAL_XP_YP)) {
fac->der.x_p = der.x_p;
fac->der.y_p = der.y_p;
}
cosphi = cos(lp.phi);
if (!(fac->code & IS_ANAL_HK)) {
fac->h = hypot(fac->der.x_p, fac->der.y_p);
fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi;
if (P->es != 0.0) {
t = sin(lp.phi);
t = 1. - P->es * t * t;
n = sqrt(t);
fac->h *= t * n / P->one_es;
fac->k *= n;
r = t * t / P->one_es;
} else
r = 1.;
} else if (P->es != 0.0) {
r = sin(lp.phi);
r = 1. - P->es * r * r;
r = r * r / P->one_es;
} else
r = 1.;
/* convergence */
if (!(fac->code & IS_ANAL_CONV)) {
fac->conv = - atan2(fac->der.y_l, fac->der.x_l);
if (fac->code & IS_ANAL_XL_YL)
fac->code |= IS_ANAL_CONV;
}
/* areal scale factor */
fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) *
r / cosphi;
/* meridian-parallel angle theta prime */
fac->thetap = aasin(P->ctx,fac->s / (fac->h * fac->k));
/* Tissot ellips axis */
t = fac->k * fac->k + fac->h * fac->h;
fac->a = sqrt(t + 2. * fac->s);
t = (t = t - 2. * fac->s) <= 0. ? 0. : sqrt(t);
fac->b = 0.5 * (fac->a - t);
fac->a = 0.5 * (fac->a + t);
/* omega */
fac->omega = 2. * aasin(P->ctx,(fac->a - fac->b)/(fac->a + fac->b));
}
return 0;
if (h < EPS)
h = DEFAULT_H;
if (fabs(lp.phi) > (M_HALFPI - h))
/* adjust to value around pi/2 where derived still exists*/
lp.phi = lp.phi < 0. ? (-M_HALFPI+h) : (M_HALFPI-h);
else if (P->geoc)
lp.phi = atan(P->rone_es * tan(lp.phi));
lp.lam -= P->lam0; /* compute del lp.lam */
if (!P->over)
lp.lam = adjlon(lp.lam); /* adjust del longitude */
if (P->spc) /* get what projection analytic values */
P->spc(lp, P, fac);
if (((fac->code & (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) !=
(IS_ANAL_XL_YL+IS_ANAL_XP_YP)) && pj_deriv(lp, h, P, &der))
return 1;

if (!(fac->code & IS_ANAL_XL_YL)) {
fac->der.x_l = der.x_l;
fac->der.y_l = der.y_l;
}
if (!(fac->code & IS_ANAL_XP_YP)) {
fac->der.x_p = der.x_p;
fac->der.y_p = der.y_p;
}
cosphi = cos(lp.phi);
if (!(fac->code & IS_ANAL_HK)) {
fac->h = hypot(fac->der.x_p, fac->der.y_p);
fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi;
if (P->es != 0.0) {
t = sin(lp.phi);
t = 1. - P->es * t * t;
n = sqrt(t);
fac->h *= t * n / P->one_es;
fac->k *= n;
r = t * t / P->one_es;
} else
r = 1.;
} else if (P->es != 0.0) {
r = sin(lp.phi);
r = 1. - P->es * r * r;
r = r * r / P->one_es;
} else
r = 1.;

/* convergence */
if (!(fac->code & IS_ANAL_CONV)) {
fac->conv = - atan2(fac->der.x_p, fac->der.y_p);
if (fac->code & IS_ANAL_XL_YL)
fac->code |= IS_ANAL_CONV;
}
/* areal scale factor */
fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) * r / cosphi;

/* meridian-parallel angle theta prime */
fac->thetap = aasin(P->ctx,fac->s / (fac->h * fac->k));

/* Tissot ellips axis */
t = fac->k * fac->k + fac->h * fac->h;
fac->a = sqrt(t + 2. * fac->s);
t = (t = t - 2. * fac->s) <= 0. ? 0. : sqrt(t);
fac->b = 0.5 * (fac->a - t);
fac->a = 0.5 * (fac->a + t);

/* omega */
fac->omega = 2. * aasin(P->ctx,(fac->a - fac->b)/(fac->a + fac->b));
}

return 0;
}

0 comments on commit ced55e8

Please sign in to comment.