Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Increase accuracy with omerc projection #317

Merged
merged 6 commits into from Jan 21, 2021
Merged

Increase accuracy with omerc projection #317

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
1523223
Initial work to increase accuracy with omerc projection
dulldrums Apr 15, 2019
9d488df
update omerc projection names, update failing omerc tests
dulldrums Apr 16, 2019
eab6854
fix issue where units other than meters would cause issues with scaling
dulldrums Apr 24, 2019
502aad0
Fix issues with Type A and Type B omerc projections, add more omerc …
dulldrums May 14, 2019
ded41da
add tests for EPSG:3468, update accuracy of some omerc Type A tests
dulldrums May 14, 2019
587ddd5
added test for omerc using unit: Foot_US
dulldrums May 14, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
327 changes: 195 additions & 132 deletions lib/projections/omerc.js
View file
Open in desktop
@@ -1,168 +1,231 @@
import tsfnz from '../common/tsfnz';
import adjust_lon from '../common/adjust_lon';
import phi2z from '../common/phi2z';
import {EPSLN, HALF_PI, FORTPI} from '../constants/values';
import { D2R, EPSLN, HALF_PI, TWO_PI, FORTPI } from '../constants/values';

var TOL = 1e-7;

function isTypeA(P) {
var typeAProjections = ['Hotine_Oblique_Mercator','Hotine_Oblique_Mercator_Azimuth_Natural_Origin'];
var projectionName = typeof P.PROJECTION === "object" ? Object.keys(P.PROJECTION)[0] : P.PROJECTION;

return 'no_uoff' in P || 'no_off' in P || typeAProjections.indexOf(projectionName) !== -1;
}


/* Initialize the Oblique Mercator projection
------------------------------------------*/
export function init() {
this.no_off = this.no_off || false;
this.no_rot = this.no_rot || false;

if (isNaN(this.k0)) {
this.k0 = 1;
export function init() {
var con, com, cosph0, D, F, H, L, sinph0, p, J, gamma = 0,
gamma0, lamc = 0, lam1 = 0, lam2 = 0, phi1 = 0, phi2 = 0, alpha_c = 0, AB;

// only Type A uses the no_off or no_uoff property
// https://github.com/OSGeo/proj.4/issues/104
this.no_off = isTypeA(this);
this.no_rot = 'no_rot' in this;

var alp = this.alpha !== 0;
var gam = (this.rectified_grid_angle && this.rectified_grid_angle !== 0) || false;

if (alp) {
alpha_c = this.alpha;
}
var sinlat = Math.sin(this.lat0);
var coslat = Math.cos(this.lat0);
var con = this.e * sinlat;

this.bl = Math.sqrt(1 + this.es / (1 - this.es) * Math.pow(coslat, 4));
this.al = this.a * this.bl * this.k0 * Math.sqrt(1 - this.es) / (1 - con * con);
var t0 = tsfnz(this.e, this.lat0, sinlat);
var dl = this.bl / coslat * Math.sqrt((1 - this.es) / (1 - con * con));
if (dl * dl < 1) {
dl = 1;

if (gam) {
gamma = (this.rectified_grid_angle * D2R);
}
var fl;
var gl;
if (!isNaN(this.longc)) {
//Central point and azimuth method

if (this.lat0 >= 0) {
fl = dl + Math.sqrt(dl * dl - 1);

if (alp || gam) {
lamc = this.longc;
} else {
lam1 = this.long1;
phi1 = this.lat1;
lam2 = this.long2;
phi2 = this.lat2;

if (Math.abs(phi1 - phi2) <= TOL || (con = Math.abs(phi1)) <= TOL ||
Math.abs(con - HALF_PI) <= TOL || Math.abs(Math.abs(this.lat0) - HALF_PI) <= TOL ||
Math.abs(Math.abs(phi2) - HALF_PI) <= TOL) {
throw new Error();
}
else {
fl = dl - Math.sqrt(dl * dl - 1);
}

var one_es = 1.0 - this.es;
com = Math.sqrt(one_es);

if (Math.abs(this.lat0) > EPSLN) {
sinph0 = Math.sin(this.lat0);
cosph0 = Math.cos(this.lat0);
con = 1 - this.es * sinph0 * sinph0;
this.B = cosph0 * cosph0;
this.B = Math.sqrt(1 + this.es * this.B * this.B / one_es);
this.A = this.B * this.k0 * com / con;
D = this.B * com / (cosph0 * Math.sqrt(con));
F = D * D -1;

if (F <= 0) {
F = 0;
} else {
F = Math.sqrt(F);
if (this.lat0 < 0) {
F = -F;
}
}
this.el = fl * Math.pow(t0, this.bl);
gl = 0.5 * (fl - 1 / fl);
this.gamma0 = Math.asin(Math.sin(this.alpha) / dl);
this.long0 = this.longc - Math.asin(gl * Math.tan(this.gamma0)) / this.bl;


this.E = F += D;
this.E *= Math.pow(tsfnz(this.e, this.lat0, sinph0), this.B);
} else {
this.B = 1 / com;
this.A = this.k0;
this.E = D = F = 1;
}
else {
//2 points method
var t1 = tsfnz(this.e, this.lat1, Math.sin(this.lat1));
var t2 = tsfnz(this.e, this.lat2, Math.sin(this.lat2));
if (this.lat0 >= 0) {
this.el = (dl + Math.sqrt(dl * dl - 1)) * Math.pow(t0, this.bl);

if (alp || gam) {
if (alp) {
gamma0 = Math.asin(Math.sin(alpha_c) / D);
if (!gam) {
gamma = alpha_c;
}
} else {
gamma0 = gamma;
alpha_c = Math.asin(D * Math.sin(gamma0));
}
else {
this.el = (dl - Math.sqrt(dl * dl - 1)) * Math.pow(t0, this.bl);
this.lam0 = lamc - Math.asin(0.5 * (F - 1 / F) * Math.tan(gamma0)) / this.B;
} else {
H = Math.pow(tsfnz(this.e, phi1, Math.sin(phi1)), this.B);
L = Math.pow(tsfnz(this.e, phi2, Math.sin(phi2)), this.B);
F = this.E / H;
p = (L - H) / (L + H);
J = this.E * this.E;
J = (J - L * H) / (J + L * H);
con = lam1 - lam2;

if (con < -Math.pi) {
lam2 -=TWO_PI;
} else if (con > Math.pi) {
lam2 += TWO_PI;
}
var hl = Math.pow(t1, this.bl);
var ll = Math.pow(t2, this.bl);
fl = this.el / hl;
gl = 0.5 * (fl - 1 / fl);
var jl = (this.el * this.el - ll * hl) / (this.el * this.el + ll * hl);
var pl = (ll - hl) / (ll + hl);
var dlon12 = adjust_lon(this.long1 - this.long2);
this.long0 = 0.5 * (this.long1 + this.long2) - Math.atan(jl * Math.tan(0.5 * this.bl * (dlon12)) / pl) / this.bl;
this.long0 = adjust_lon(this.long0);
var dlon10 = adjust_lon(this.long1 - this.long0);
this.gamma0 = Math.atan(Math.sin(this.bl * (dlon10)) / gl);
this.alpha = Math.asin(dl * Math.sin(this.gamma0));

this.lam0 = adjust_lon(0.5 * (lam1 + lam2) - Math.atan(J * Math.tan(0.5 * this.B * (lam1 - lam2)) / p) / this.B);
gamma0 = Math.atan(2 * Math.sin(this.B * adjust_lon(lam1 - this.lam0)) / (F - 1 / F));
gamma = alpha_c = Math.asin(D * Math.sin(gamma0));
}


this.singam = Math.sin(gamma0);
this.cosgam = Math.cos(gamma0);
this.sinrot = Math.sin(gamma);
this.cosrot = Math.cos(gamma);

this.rB = 1 / this.B;
this.ArB = this.A * this.rB;
this.BrA = 1 / this.ArB;
AB = this.A * this.B;

if (this.no_off) {
this.uc = 0;
this.u_0 = 0;
} else {
this.u_0 = Math.abs(this.ArB * Math.atan(Math.sqrt(D * D - 1) / Math.cos(alpha_c)));

if (this.lat0 < 0) {
this.u_0 = - this.u_0;
}
}
else {
if (this.lat0 >= 0) {
this.uc = this.al / this.bl * Math.atan2(Math.sqrt(dl * dl - 1), Math.cos(this.alpha));
}
else {
this.uc = -1 * this.al / this.bl * Math.atan2(Math.sqrt(dl * dl - 1), Math.cos(this.alpha));
}
}


F = 0.5 * gamma0;
this.v_pole_n = this.ArB * Math.log(Math.tan(FORTPI - F));
this.v_pole_s = this.ArB * Math.log(Math.tan(FORTPI + F));
}


/* Oblique Mercator forward equations--mapping lat,long to x,y
----------------------------------------------------------*/
export function forward(p) {
var lon = p.x;
var lat = p.y;
var dlon = adjust_lon(lon - this.long0);
var us, vs;
var con;
if (Math.abs(Math.abs(lat) - HALF_PI) <= EPSLN) {
if (lat > 0) {
con = -1;
var coords = {};
var S, T, U, V, W, temp, u, v;
p.x = p.x - this.lam0;

if (Math.abs(Math.abs(p.y) - HALF_PI) > EPSLN) {
W = this.E / Math.pow(tsfnz(this.e, p.y, Math.sin(p.y)), this.B);

temp = 1 / W;
S = 0.5 * (W - temp);
T = 0.5 * (W + temp);
V = Math.sin(this.B * p.x);
U = (S * this.singam - V * this.cosgam) / T;

if (Math.abs(Math.abs(U) - 1.0) < EPSLN) {
throw new Error();
}
else {
con = 1;
}
vs = this.al / this.bl * Math.log(Math.tan(FORTPI + con * this.gamma0 * 0.5));
us = -1 * con * HALF_PI * this.al / this.bl;

v = 0.5 * this.ArB * Math.log((1 - U)/(1 + U));
temp = Math.cos(this.B * p.x);

if (Math.abs(temp) < TOL) {
u = this.A * p.x;
} else {
u = this.ArB * Math.atan2((S * this.cosgam + V * this.singam), temp);
}
} else {
v = p.y > 0 ? this.v_pole_n : this.v_pole_s;
u = this.ArB * p.y;
}
else {
var t = tsfnz(this.e, lat, Math.sin(lat));
var ql = this.el / Math.pow(t, this.bl);
var sl = 0.5 * (ql - 1 / ql);
var tl = 0.5 * (ql + 1 / ql);
var vl = Math.sin(this.bl * (dlon));
var ul = (sl * Math.sin(this.gamma0) - vl * Math.cos(this.gamma0)) / tl;
if (Math.abs(Math.abs(ul) - 1) <= EPSLN) {
vs = Number.POSITIVE_INFINITY;
}
else {
vs = 0.5 * this.al * Math.log((1 - ul) / (1 + ul)) / this.bl;
}
if (Math.abs(Math.cos(this.bl * (dlon))) <= EPSLN) {
us = this.al * this.bl * (dlon);
}
else {
us = this.al * Math.atan2(sl * Math.cos(this.gamma0) + vl * Math.sin(this.gamma0), Math.cos(this.bl * dlon)) / this.bl;
}
}


if (this.no_rot) {
p.x = this.x0 + us;
p.y = this.y0 + vs;
coords.x = u;
coords.y = v;
} else {
u -= this.u_0;
coords.x = v * this.cosrot + u * this.sinrot;
coords.y = u * this.cosrot - v * this.sinrot;
}
else {

us -= this.uc;
p.x = this.x0 + vs * Math.cos(this.alpha) + us * Math.sin(this.alpha);
p.y = this.y0 + us * Math.cos(this.alpha) - vs * Math.sin(this.alpha);
}
return p;

coords.x = (this.a * coords.x + this.x0);
coords.y = (this.a * coords.y + this.y0);

return coords;
}

export function inverse(p) {
var us, vs;
var u, v, Qp, Sp, Tp, Vp, Up;
var coords = {};

p.x = (p.x - this.x0) * (1.0 / this.a);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are we going to need to add some tests for non meter omerc projections ?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would think so.

I'll get back to this PR sometime next week, got pulled into some other stuff.

p.y = (p.y - this.y0) * (1.0 / this.a);

if (this.no_rot) {
vs = p.y - this.y0;
us = p.x - this.x0;
}
else {
vs = (p.x - this.x0) * Math.cos(this.alpha) - (p.y - this.y0) * Math.sin(this.alpha);
us = (p.y - this.y0) * Math.cos(this.alpha) + (p.x - this.x0) * Math.sin(this.alpha);
us += this.uc;
v = p.y;
u = p.x;
} else {
v = p.x * this.cosrot - p.y * this.sinrot;
u = p.y * this.cosrot + p.x * this.sinrot + this.u_0;
}
var qp = Math.exp(-1 * this.bl * vs / this.al);
var sp = 0.5 * (qp - 1 / qp);
var tp = 0.5 * (qp + 1 / qp);
var vp = Math.sin(this.bl * us / this.al);
var up = (vp * Math.cos(this.gamma0) + sp * Math.sin(this.gamma0)) / tp;
var ts = Math.pow(this.el / Math.sqrt((1 + up) / (1 - up)), 1 / this.bl);
if (Math.abs(up - 1) < EPSLN) {
p.x = this.long0;
p.y = HALF_PI;
}
else if (Math.abs(up + 1) < EPSLN) {
p.x = this.long0;
p.y = -1 * HALF_PI;
}
else {
p.y = phi2z(this.e, ts);
p.x = adjust_lon(this.long0 - Math.atan2(sp * Math.cos(this.gamma0) - vp * Math.sin(this.gamma0), Math.cos(this.bl * us / this.al)) / this.bl);

Qp = Math.exp(-this.BrA * v);
Sp = 0.5 * (Qp - 1 / Qp);
Tp = 0.5 * (Qp + 1 / Qp);
Vp = Math.sin(this.BrA * u);
Up = (Vp * this.cosgam + Sp * this.singam) / Tp;

if (Math.abs(Math.abs(Up) - 1) < EPSLN) {
coords.x = 0;
coords.y = Up < 0 ? -HALF_PI : HALF_PI;
} else {
coords.y = this.E / Math.sqrt((1 + Up) / (1 - Up));
coords.y = phi2z(this.e, Math.pow(coords.y, 1 / this.B));

if (coords.y === Infinity) {
throw new Error();
}

coords.x = -this.rB * Math.atan2((Sp * this.cosgam - Vp * this.singam), Math.cos(this.BrA * u));
}
return p;

coords.x += this.lam0;

return coords;
}

export var names = ["Hotine_Oblique_Mercator", "Hotine Oblique Mercator", "Hotine_Oblique_Mercator_Azimuth_Natural_Origin", "Hotine_Oblique_Mercator_Azimuth_Center", "omerc"];
export var names = ["Hotine_Oblique_Mercator", "Hotine Oblique Mercator", "Hotine_Oblique_Mercator_Azimuth_Natural_Origin", "Hotine_Oblique_Mercator_Two_Point_Natural_Origin", "Hotine_Oblique_Mercator_Azimuth_Center", "Oblique_Mercator", "omerc"];
export default {
init: init,
forward: forward,
Expand Down