From 1523223350086e3760661293c17117ceebd6ec5e Mon Sep 17 00:00:00 2001 From: Ryan Sims Date: Mon, 15 Apr 2019 18:15:35 -0400 Subject: [PATCH 1/6] Initial work to increase accuracy with omerc projection --- lib/projections/omerc.js | 318 +++++++++++++++++++++++---------------- test/testData.js | 5 + 2 files changed, 192 insertions(+), 131 deletions(-) diff --git a/lib/projections/omerc.js b/lib/projections/omerc.js index 3aa9c445..c5e6d485 100644 --- a/lib/projections/omerc.js +++ b/lib/projections/omerc.js @@ -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, diff --git a/test/testData.js b/test/testData.js index 78db2e28..8206844a 100644 --- a/test/testData.js +++ b/test/testData.js @@ -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], From 9d488df440f084dbb74c17220d745245dd758f60 Mon Sep 17 00:00:00 2001 From: Ryan Sims Date: Tue, 16 Apr 2019 15:02:18 -0400 Subject: [PATCH 2/6] update omerc projection names, update failing omerc tests --- lib/projections/omerc.js | 2 +- test/testData.js | 51 +++++++++++++++++++++++++++++----------- 2 files changed, 38 insertions(+), 15 deletions(-) diff --git a/lib/projections/omerc.js b/lib/projections/omerc.js index c5e6d485..74e83ff1 100644 --- a/lib/projections/omerc.js +++ b/lib/projections/omerc.js @@ -218,7 +218,7 @@ export function inverse(p) { } // 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 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, diff --git a/test/testData.js b/test/testData.js index 8206844a..bb719e25 100644 --- a/test/testData.js +++ b/test/testData.js @@ -8,17 +8,19 @@ var testPoints = [ ll: [37.33761240175515, 55.60447049026976] }, {code: 'PROJCS["CH1903 / LV03",GEOGCS["CH1903",DATUM["D_CH1903",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER["azimuth",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",600000],PARAMETER["false_northing",200000],UNIT["Meter",1]]', - xy: [660389.4751110513, 185731.68482649108], - ll: [8.23, 46.82], + xy: [660013.4882918689, 185172.17110117766], + ll: [8.225, 46.815], acc:{ - xy:0.1 + xy: 0.1, + ll: 5 } }, - {code: 'PROJCS["CH1903 / LV03",GEOGCS["CH1903",DATUM["CH1903",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[674.374,15.056,405.346,0,0,0,0],AUTHORITY["EPSG","6149"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4149"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Hotine_Oblique_Mercator"],PARAMETER["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER["azimuth",90],PARAMETER["rectified_grid_angle",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",600000],PARAMETER["false_northing",200000],AUTHORITY["EPSG","21781"],AXIS["Y",EAST],AXIS["X",NORTH]]', - xy: [660389.4751110513, 185731.68482649108], - ll: [8.23, 46.82], + {code: 'PROJCS["CH1903 / LV03",GEOGCS["CH1903",DATUM["CH1903",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[674.4,15.1,405.3,0,0,0,0],AUTHORITY["EPSG","6149"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4149"]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER["azimuth",90],PARAMETER["rectified_grid_angle",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",600000],PARAMETER["false_northing",200000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Y",EAST],AXIS["X",NORTH],AUTHORITY["EPSG","21781"]]', + xy: [660013.4882918689, 185172.17110117766], + ll: [8.225, 46.815], acc:{ - xy:0.1 + xy: 0.1, + ll: 5 } }, {code: 'PROJCS["NAD83 / Massachusetts Mainland",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",42.68333333333333],PARAMETER["standard_parallel_2",41.71666666666667],PARAMETER["latitude_of_origin",41],PARAMETER["central_meridian",-71.5],PARAMETER["false_easting",200000],PARAMETER["false_northing",750000],AUTHORITY["EPSG","26986"],AXIS["X",EAST],AXIS["Y",NORTH]]', @@ -192,12 +194,20 @@ var testPoints = [ xy:[2464770.343667, 6056137.861919] },{ code:'PROJCS["Rassadiran / Nakhl e Taqi",GEOGCS["Rassadiran",DATUM["Rassadiran",SPHEROID["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-133.63,-157.5,-158.62,0,0,0,0],AUTHORITY["EPSG","6153"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4153"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Hotine_Oblique_Mercator"],PARAMETER["latitude_of_center",27.51882880555555],PARAMETER["longitude_of_center",52.60353916666667],PARAMETER["azimuth",0.5716611944444444],PARAMETER["rectified_grid_angle",0.5716611944444444],PARAMETER["scale_factor",0.999895934],PARAMETER["false_easting",658377.437],PARAMETER["false_northing",3044969.194],AUTHORITY["EPSG","2057"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]', - ll:[52.6, 27.5], - xy:[658017.25458, 3043003.058818] + ll: [52.605, 27.5], + xy: [658511.2619462232, 3043003.0546802087], + acc:{ + ll: 8, + xy: 8 + } },{ code:'PROJCS["Rassadiran / Nakhl e Taqi",GEOGCS["Rassadiran",DATUM["D_Rassadiran",SPHEROID["International_1924",6378388,297]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Natural_Origin"],PARAMETER["latitude_of_center",27.51882880555555],PARAMETER["longitude_of_center",52.60353916666667],PARAMETER["azimuth",0.5716611944444444],PARAMETER["rectified_grid_angle",0.5716611944444444],PARAMETER["scale_factor",0.999895934],PARAMETER["false_easting",658377.437],PARAMETER["false_northing",3044969.194],UNIT["Meter",1]]', - ll:[52.6, 27.5], - xy:[658017.25458, 3043003.058818] + ll:[52.605, 27.5], + xy:[658511.2619462232, 3043003.0546802087], + acc:{ + ll: 8, + xy: 8 + } },{ code:'PROJCS["SAD69 / Brazil Polyconic",GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_SAD69",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Polyconic"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-54],PARAMETER["false_easting",5000000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]', ll:[-49.221772553812, -0.34551739237581], @@ -554,9 +564,22 @@ var testPoints = [ 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["CUSTOM_OBLIQUE_MERCATOR", 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: [-361.2499999983702, 254.91500000283122], + ll: [-122.25064809, 37.50832038], + acc:{ + ll: 3, + xy: 8 + } + }, + { + code: 'PROJCS["UNK / Oblique_Mercator",GEOGCS["UNK",DATUM["Unknown datum",SPHEROID["WGS 84", 6378137.0, 298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.017453292519943295]],PROJECTION["Oblique_Mercator"],PARAMETER["latitude_of_center",37.4769061],PARAMETER["longitude_of_center",141.0039618],PARAMETER["central_meridian",141.0039618],PARAMETER["azimuth",202.22],PARAMETER["scale_factor",1],PARAMETER["false_easting",138],PARAMETER["false_northing",77.65],UNIT["Meter",1]]', + xy: [168.2438, 64.1736], + ll: [141.003611, 37.476802], + acc:{ + ll: 9, + xy: 4 + } }, { 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"]]', From eab68540c770dc494ebd1cc4fe7cdb9474031edb Mon Sep 17 00:00:00 2001 From: Ryan Sims Date: Wed, 24 Apr 2019 19:42:09 -0400 Subject: [PATCH 3/6] fix issue where units other than meters would cause issues with scaling --- lib/projections/omerc.js | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/projections/omerc.js b/lib/projections/omerc.js index 74e83ff1..57719de6 100644 --- a/lib/projections/omerc.js +++ b/lib/projections/omerc.js @@ -181,8 +181,8 @@ export function inverse(p) { 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); + p.x = (p.x - this.x0) * (1.0 / this.a); + p.y = (p.y - this.y0) * (1.0 / this.a); if (this.no_rot) { v = p.y; From 502aad06962b915c386881d96e758e546280b443 Mon Sep 17 00:00:00 2001 From: Ryan Sims Date: Tue, 14 May 2019 12:57:07 -0400 Subject: [PATCH 4/6] Fix issues with Type A and Type B omerc projections, add more omerc tests --- lib/projections/omerc.js | 19 +++++++++----- test/testData.js | 57 +++++++++++++++++++++++++++------------- 2 files changed, 52 insertions(+), 24 deletions(-) diff --git a/lib/projections/omerc.js b/lib/projections/omerc.js index 57719de6..6b8bd3b1 100644 --- a/lib/projections/omerc.js +++ b/lib/projections/omerc.js @@ -5,15 +5,24 @@ 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() { 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; + // 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; @@ -112,7 +121,7 @@ export function init() { this.BrA = 1 / this.ArB; AB = this.A * this.B; - if (no_off) { + if (this.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))); @@ -136,7 +145,6 @@ export function forward(p) { 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; @@ -217,7 +225,6 @@ export function inverse(p) { return coords; } -// 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", "Oblique_Mercator", "omerc"]; export default { init: init, diff --git a/test/testData.js b/test/testData.js index bb719e25..91807d57 100644 --- a/test/testData.js +++ b/test/testData.js @@ -117,12 +117,12 @@ var testPoints = [ xy:-4 } }, - /*{ - code:'EPSG:3975', - ll:[-9.764450683, 25.751953], - xy:[-942135.525095996, 3178441.8667094777] - },*/ - { + // { + // code:'EPSG:3975', + // ll:[-9.764450683, 25.751953], + // xy:[-942135.525095996, 3178441.8667094777] + // }, + { code:'PROJCS["World Equidistant Cylindrical (Sphere)",GEOGCS["Unspecified datum based upon the GRS 1980 Authalic Sphere",DATUM["Not_specified_based_on_GRS_1980_Authalic_Sphere",SPHEROID["GRS 1980 Authalic Sphere",6371007,0,AUTHORITY["EPSG","7048"]],AUTHORITY["EPSG","6047"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4047"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Equirectangular"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3786"],AXIS["X",EAST],AXIS["Y",NORTH]]', ll:[-1.7539371169976, 12.632997701986], xy:[-195029.12334755991, 1395621.9368162225], @@ -193,20 +193,12 @@ var testPoints = [ ll:[172.465, -40.7], xy:[2464770.343667, 6056137.861919] },{ - code:'PROJCS["Rassadiran / Nakhl e Taqi",GEOGCS["Rassadiran",DATUM["Rassadiran",SPHEROID["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-133.63,-157.5,-158.62,0,0,0,0],AUTHORITY["EPSG","6153"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4153"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Hotine_Oblique_Mercator"],PARAMETER["latitude_of_center",27.51882880555555],PARAMETER["longitude_of_center",52.60353916666667],PARAMETER["azimuth",0.5716611944444444],PARAMETER["rectified_grid_angle",0.5716611944444444],PARAMETER["scale_factor",0.999895934],PARAMETER["false_easting",658377.437],PARAMETER["false_northing",3044969.194],AUTHORITY["EPSG","2057"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]', + code: 'PROJCS["Rassadiran / Nakhl e Taqi", GEOGCS["Rassadiran", DATUM["Rassadiran", SPHEROID["International 1924",6378388,297, AUTHORITY["EPSG","7022"]], TOWGS84[-133.63,-157.5,-158.62,0,0,0,0], AUTHORITY["EPSG","6153"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4153"]], PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"], PARAMETER["latitude_of_center",27.51882880555555], PARAMETER["longitude_of_center",52.60353916666667], PARAMETER["azimuth",0.5716611944444444], PARAMETER["rectified_grid_angle",0.5716611944444444], PARAMETER["scale_factor",0.999895934], PARAMETER["false_easting",658377.437], PARAMETER["false_northing",3044969.194], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AXIS["Easting",EAST], AXIS["Northing",NORTH], AUTHORITY["EPSG","2057"]]', ll: [52.605, 27.5], - xy: [658511.2619462232, 3043003.0546802087], - acc:{ - ll: 8, - xy: 8 - } - },{ - code:'PROJCS["Rassadiran / Nakhl e Taqi",GEOGCS["Rassadiran",DATUM["D_Rassadiran",SPHEROID["International_1924",6378388,297]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Natural_Origin"],PARAMETER["latitude_of_center",27.51882880555555],PARAMETER["longitude_of_center",52.60353916666667],PARAMETER["azimuth",0.5716611944444444],PARAMETER["rectified_grid_angle",0.5716611944444444],PARAMETER["scale_factor",0.999895934],PARAMETER["false_easting",658377.437],PARAMETER["false_northing",3044969.194],UNIT["Meter",1]]', - ll:[52.605, 27.5], - xy:[658511.2619462232, 3043003.0546802087], - acc:{ + xy: [658511.261946, 3043003.05468], + acc: { ll: 8, - xy: 8 + xy: 6 } },{ code:'PROJCS["SAD69 / Brazil Polyconic",GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_SAD69",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Polyconic"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-54],PARAMETER["false_easting",5000000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]', @@ -572,6 +564,35 @@ var testPoints = [ xy: 8 } }, + // Omerc Type A - #273 + { + code: '+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257964666666 +k=0.99984 +x_0=804671 +y_0=0 +no_uoff +gamma=323.1301023611111 +ellps=GRS80 +units=m +no_defs', + xy: [412597.532715, 338944.957259], + ll: [101.70979078430528, 3.06268465621428], + acc:{ + ll: 2, + xy: -4 + } + }, + { + code: 'PROJCS["GDM2000 / Peninsula RSO", GEOGCS["GDM2000", DATUM["Geodetic_Datum_of_Malaysia_2000", SPHEROID["GRS 1980",6378137,298.257222101, AUTHORITY["EPSG","7019"]], AUTHORITY["EPSG","6742"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4742"]], PROJECTION["Hotine_Oblique_Mercator"], PARAMETER["latitude_of_center",4], PARAMETER["longitude_of_center",102.25], PARAMETER["azimuth",323.0257964666666], PARAMETER["rectified_grid_angle",323.1301023611111], PARAMETER["scale_factor",0.99984], PARAMETER["false_easting",804671], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AXIS["Easting",EAST], AXIS["Northing",NORTH], AUTHORITY["EPSG","3375"]]', + xy: [412597.532715, 338944.957259], + ll: [101.70979078430528, 3.06268465621428], + acc:{ + ll: 2, + xy: -4 + } + }, + // Omerc Type B - #308 + { + code: '+proj=omerc +lat_0=37.4769061 +lonc=141.0039618 +alpha=202.22 +k=1 +x_0=138 +y_0=77.65 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs', + xy: [168.2438, 64.1736], + ll: [141.003611, 37.476802], + acc:{ + ll: 9, + xy: 4 + } + }, { code: 'PROJCS["UNK / Oblique_Mercator",GEOGCS["UNK",DATUM["Unknown datum",SPHEROID["WGS 84", 6378137.0, 298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.017453292519943295]],PROJECTION["Oblique_Mercator"],PARAMETER["latitude_of_center",37.4769061],PARAMETER["longitude_of_center",141.0039618],PARAMETER["central_meridian",141.0039618],PARAMETER["azimuth",202.22],PARAMETER["scale_factor",1],PARAMETER["false_easting",138],PARAMETER["false_northing",77.65],UNIT["Meter",1]]', xy: [168.2438, 64.1736], From ded41dafddbd57fede71dfdd104c78cd49560add Mon Sep 17 00:00:00 2001 From: Ryan Sims Date: Tue, 14 May 2019 13:28:48 -0400 Subject: [PATCH 5/6] add tests for EPSG:3468, update accuracy of some omerc Type A tests --- test/testData.js | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/test/testData.js b/test/testData.js index 91807d57..3dabfc25 100644 --- a/test/testData.js +++ b/test/testData.js @@ -571,7 +571,7 @@ var testPoints = [ ll: [101.70979078430528, 3.06268465621428], acc:{ ll: 2, - xy: -4 + xy: -3 } }, { @@ -579,8 +579,27 @@ var testPoints = [ xy: [412597.532715, 338944.957259], ll: [101.70979078430528, 3.06268465621428], acc:{ - ll: 2, - xy: -4 + ll: 7, + xy: 6 + } + }, + // EPSG:3468 + { + code: '+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000 +no_uoff +gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs', + xy: [1264314.74, -763162.04], + ll: [-128.115000029, 44.8150000066], + acc:{ + ll: 9, + xy: 4 + } + }, + { + code: 'PROJCS["NAD83(NSRS2007) / Alaska zone 1", GEOGCS["NAD83(NSRS2007)", DATUM["NAD83_National_Spatial_Reference_System_2007", SPHEROID["GRS 1980",6378137,298.257222101, AUTHORITY["EPSG","7019"]], TOWGS84[0,0,0,0,0,0,0], AUTHORITY["EPSG","6759"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4759"]], PROJECTION["Hotine_Oblique_Mercator"], PARAMETER["latitude_of_center",57], PARAMETER["longitude_of_center",-133.6666666666667], PARAMETER["azimuth",323.1301023611111], PARAMETER["rectified_grid_angle",323.1301023611111], PARAMETER["scale_factor",0.9999], PARAMETER["false_easting",5000000], PARAMETER["false_northing",-5000000], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AXIS["X",EAST], AXIS["Y",NORTH], AUTHORITY["EPSG","3468"]]', + xy: [1264314.74, -763162.04], + ll: [-128.115000029, 44.8150000066], + acc:{ + ll: 9, + xy: 4 } }, // Omerc Type B - #308 From 587ddd50fd78bbee1f63f79f33606ac307da0c89 Mon Sep 17 00:00:00 2001 From: Ryan Sims Date: Tue, 14 May 2019 13:52:15 -0400 Subject: [PATCH 6/6] added test for omerc using unit: Foot_US --- test/testData.js | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/testData.js b/test/testData.js index 3dabfc25..34d29f07 100644 --- a/test/testData.js +++ b/test/testData.js @@ -621,6 +621,12 @@ var testPoints = [ xy: 4 } }, + // Test with Feet + { + code: 'PROJCS["UNK / Oblique_Mercator",GEOGCS["UNK",DATUM["Unknown datum",SPHEROID["WGS 84", 6378137.0, 298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.017453292519943295]],PROJECTION["Oblique_Mercator"],PARAMETER["latitude_of_center",37.4769061],PARAMETER["longitude_of_center",141.0039618],PARAMETER["central_meridian",141.0039618],PARAMETER["azimuth",202.22],PARAMETER["scale_factor",1],PARAMETER["false_easting",138],PARAMETER["false_northing",77.65],UNIT["Foot_US",0.3048006096012192]]', + xy: [237.22488871325027, 33.43626458451221], + ll: [141.003611, 37.476802], + }, { 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],