Skip to content

Commit

Permalink
Initial work to increase accuracy with omerc projection
Browse files Browse the repository at this point in the history
  • Loading branch information
@dulldrums
dulldrums committed Apr 16, 2019
1 parent 6ebcf9c commit 1523223
Show file tree
Hide file tree
Showing 2 changed files with 192 additions and 131 deletions.
318 changes: 187 additions & 131 deletions lib/projections/omerc.js
View file
@@ -1,168 +1,224 @@
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;

/* Initialize the Oblique Mercator projection
------------------------------------------*/
export function init() {
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;
var no_off = 0;

this.no_off = this.no_off || false;
this.no_rot = this.no_rot || false;

if (isNaN(this.k0)) {
this.k0 = 1;

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);
}
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;

}
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);
}
else {
this.el = (dl - Math.sqrt(dl * dl - 1)) * Math.pow(t0, this.bl);

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;
}
}
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.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;
}

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

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.uc = -1 * this.al / this.bl * Math.atan2(Math.sqrt(dl * dl - 1), Math.cos(this.alpha));
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;
}

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 (no_off) {
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;
}
}

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;
}
else {
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) {
// TODO: this is sometimes resulting in a negative number to a positive power, which is NaN
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();
}
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.to_meter - this.x0) * (1.0 / this.a);
p.y = (p.y * this.to_meter - 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"];
// TODO: confirm all of these are the correct supported coordinates
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", "omerc"];
export default {
init: init,
forward: forward,
Expand Down
5 changes: 5 additions & 0 deletions test/testData.js
View file
Expand Up @@ -553,6 +553,11 @@ var testPoints = [
ll: [2, 1],
xy: [-959006.4926646841, 113457.31956265299]
},
{
code: 'PROJCS["CUSTOM_OBLIQUE_MERCATOR:37.50832038,-122.25064809,-361.2500,254.9150,37.50644797,-122.24946191,382.0320,65.0710,ft", GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], AXIS["Geodetic latitude", NORTH], AXIS["Geodetic longitude", EAST]], PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center", AUTHORITY["EPSG", "9815"]], PARAMETER["latitude_of_center", 37.50832038], PARAMETER["longitude_of_center", -122.25064809], PARAMETER["azimuth", 45.0], PARAMETER["rectified_grid_angle", -3.99], PARAMETER["scale_factor", 1.0], PARAMETER["false_easting", -361.25], PARAMETER["false_northing", 254.915], UNIT["foot", 0.3048], AXIS["Easting", EAST], AXIS["Northing", NORTH]]',
xy: [397.4667002959177, -3861.3150102049112],
ll: [-122.25963848377265, 37.49932998622735]
},
{
code: 'PROJCS["WGS 84 / Pseudo-Mercator", GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 0, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic latitude", NORTH], AXIS["Geodetic longitude", EAST], AUTHORITY["EPSG","4326"]], PROJECTION["Popular Visualisation Pseudo Mercator", AUTHORITY["EPSG","1024"]], PARAMETER["semi_minor", 6378137.0], PARAMETER["latitude_of_origin", 0.0], PARAMETER["central_meridian", 0.0], PARAMETER["scale_factor", 1.0], PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], UNIT["m", 1.0], AXIS["Easting", EAST], AXIS["Northing", NORTH], AUTHORITY["EPSG","3857"]]',
xy: [-12523490.49256873, 5166512.50707369],
Expand Down

0 comments on commit 1523223

Please sign in to comment.