diff --git a/Gruntfile.js b/Gruntfile.js index d64b071..6553cc5 100644 --- a/Gruntfile.js +++ b/Gruntfile.js @@ -29,7 +29,8 @@ var projs = [ 'tpers', 'geos', 'eqearth', - 'bonne' + 'bonne', + 'ob_tran' ]; module.exports = function (grunt) { grunt.initConfig({ diff --git a/lib/defs.js b/lib/defs.js index 53332c0..e9809cc 100644 --- a/lib/defs.js +++ b/lib/defs.js @@ -38,6 +38,7 @@ import wkt from 'wkt-parser'; * @property {number} [rectified_grid_angle] * @property {boolean} [approx] * @property {boolean} [over] + * @property {string} [projStr] * @property {(coordinates: T, enforceAxis?: boolean) => T} inverse * @property {(coordinates: T, enforceAxis?: boolean) => T} forward */ diff --git a/lib/includedProjections.js b/lib/includedProjections.js index fa4da90..c6a740e 100644 --- a/lib/includedProjections.js +++ b/lib/includedProjections.js @@ -28,6 +28,7 @@ import tpers from './projections/tpers'; import geos from './projections/geos'; import eqearth from './projections/eqearth'; import bonne from './projections/bonne'; +import ob_tran from './projections/ob_tran'; var projs = [ tmerc, @@ -59,7 +60,8 @@ var projs = [ tpers, geos, eqearth, - bonne + bonne, + ob_tran ]; export default function (proj4) { diff --git a/lib/projString.js b/lib/projString.js index c9a6854..04eb138 100644 --- a/lib/projString.js +++ b/lib/projString.js @@ -146,5 +146,6 @@ export default function (defData) { if (typeof self.datumCode === 'string' && self.datumCode !== 'WGS84') { self.datumCode = self.datumCode.toLowerCase(); } + self['projStr'] = defData; return self; } diff --git a/lib/projections/eqearth.js b/lib/projections/eqearth.js index de0a80e..1a330b0 100644 --- a/lib/projections/eqearth.js +++ b/lib/projections/eqearth.js @@ -38,6 +38,8 @@ var A1 = 1.340264, export function init() { this.es = 0; this.long0 = this.long0 !== undefined ? this.long0 : 0; + this.x0 = this.x0 !== undefined ? this.x0 : 0; + this.y0 = this.y0 !== undefined ? this.y0 : 0; } export function forward(p) { diff --git a/lib/projections/moll.js b/lib/projections/moll.js index ee860c9..912d028 100644 --- a/lib/projections/moll.js +++ b/lib/projections/moll.js @@ -1,6 +1,13 @@ import adjust_lon from '../common/adjust_lon'; -export function init() {} import { EPSLN } from '../constants/values'; + +/** @this {import('../defs.js').ProjectionDefinition} */ +export function init() { + this.x0 = this.x0 !== undefined ? this.x0 : 0; + this.y0 = this.y0 !== undefined ? this.y0 : 0; + this.long0 = this.long0 !== undefined ? this.long0 : 0; +} + /* Mollweide forward equations--mapping lat,long to x,y ---------------------------------------------------- */ export function forward(p) { diff --git a/lib/projections/ob_tran.js b/lib/projections/ob_tran.js new file mode 100644 index 0000000..42369bf --- /dev/null +++ b/lib/projections/ob_tran.js @@ -0,0 +1,347 @@ +import adjust_lon from '../common/adjust_lon'; +import { D2R, HALF_PI, EPSLN } from '../constants/values'; +import Proj from '../Proj'; + +/** + Original projection implementation: + https://github.com/OSGeo/PROJ/blob/46c47e9adf6376ae06afabe5d24a0016a05ced82/src/projections/ob_tran.cpp + + Documentation: + https://proj.org/operations/projections/ob_tran.html + + References/Formulas: + https://pubs.usgs.gov/pp/1395/report.pdf + + Examples: + +proj=ob_tran +o_proj=moll +o_lat_p=45 +o_lon_p=-90 + +proj=ob_tran +o_proj=moll +o_lat_p=45 +o_lon_p=-90 +lon_0=60 + +proj=ob_tran +o_proj=moll +o_lat_p=45 +o_lon_p=-90 +lon_0=-90 +*/ + +const projectionType = { + OBLIQUE: { + forward: forwardOblique, + inverse: inverseOblique + }, + TRANSVERSE: { + forward: forwardTransverse, + inverse: inverseTransverse + } +}; + +/** + * @typedef {Object} LocalThis + * @property {number} lamp + * @property {number} cphip + * @property {number} sphip + * @property {Object} projectionType + * @property {string | undefined} o_proj + * @property {string | undefined} o_lon_p + * @property {string | undefined} o_lat_p + * @property {string | undefined} o_alpha + * @property {string | undefined} o_lon_c + * @property {string | undefined} o_lat_c + * @property {string | undefined} o_lon_1 + * @property {string | undefined} o_lat_1 + * @property {string | undefined} o_lon_2 + * @property {string | undefined} o_lat_2 + * @property {number | undefined} oLongP + * @property {number | undefined} oLatP + * @property {number | undefined} oAlpha + * @property {number | undefined} oLongC + * @property {number | undefined} oLatC + * @property {number | undefined} oLong1 + * @property {number | undefined} oLat1 + * @property {number | undefined} oLong2 + * @property {number | undefined} oLat2 + * @property {import('..').Converter} obliqueProjection + * + */ + +/** + * Parameters can be from the following sets: + * New pole --> o_lat_p, o_lon_p + * Rotate about point --> o_alpha, o_lon_c, o_lat_c + * New equator points --> lon_1, lat_1, lon_2, lat_2 + * + * Per the original source code, the parameter sets are + * checked in the order of the object below. + */ +const paramSets = { + ROTATE: { + o_alpha: 'oAlpha', + o_lon_c: 'oLongC', + o_lat_c: 'oLatC' + }, + NEW_POLE: { + o_lat_p: 'oLatP', + o_lon_p: 'oLongP' + }, + NEW_EQUATOR: { + o_lon_1: 'oLong1', + o_lat_1: 'oLat1', + o_lon_2: 'oLong2', + o_lat_2: 'oLat2' + } +}; + +/** @this {import('../defs.js').ProjectionDefinition & LocalThis} */ +export function init() { + this.x0 = this.x0 || 0; + this.y0 = this.y0 || 0; + this.long0 = this.long0 || 0; + this.title = this.title || 'General Oblique Transformation'; + + /** Verify required parameters exist */ + if (!this.o_proj) { + throw new Error('Missing parameter: o_proj'); + } + + if (this.o_proj === `ob_tran`) { + throw new Error('Invalid value for o_proj: ' + this.o_proj); + } + + const newProjStr = this.projStr.replace('+proj=ob_tran', '').replace('+o_proj=', '+proj=').trim(); + + /** @type {import('../defs.js').ProjectionDefinition} */ + const oProj = Proj(newProjStr); + if (!oProj) { + throw new Error('Invalid parameter: o_proj. Unknown projection ' + this.o_proj); + } + oProj.long0 = 0; // we handle long0 before/after forward/inverse + this.obliqueProjection = oProj; + + let matchedSet; + const paramSetsKeys = Object.keys(paramSets); + + /** + * parse strings, convert to radians, throw on NaN + * @param {string} name + * @returns {number | undefined} + */ + const parseParam = (name) => { + if (typeof this[name] === `undefined`) { + return undefined; + } + const val = parseFloat(this[name]) * D2R; + if (isNaN(val)) { + throw new Error('Invalid value for ' + name + ': ' + this[name]); + } + return val; + }; + + for (let i = 0; i < paramSetsKeys.length; i++) { + const setKey = paramSetsKeys[i]; + const set = paramSets[setKey]; + const params = Object.entries(set); + const setHasParams = params.some( + ([p]) => typeof this[p] !== 'undefined' + ); + if (!setHasParams) { + continue; + } + matchedSet = set; + for (let ii = 0; ii < params.length; ii++) { + const [inputParam, param] = params[ii]; + const val = parseParam(inputParam); + if (typeof val === 'undefined') { + throw new Error('Missing parameter: ' + inputParam + '.'); + } + this[param] = val; + } + break; + } + + if (!matchedSet) { + throw new Error('No valid parameters provided for ob_tran projection.'); + } + + const { lamp, phip } = createRotation(this, matchedSet); + this.lamp = lamp; + + if (Math.abs(phip) > EPSLN) { + this.cphip = Math.cos(phip); + this.sphip = Math.sin(phip); + this.projectionType = projectionType.OBLIQUE; + } else { + this.projectionType = projectionType.TRANSVERSE; + } +} + +// ob_tran forward equations--mapping (lat,long) to (x,y) +// transverse (90 degrees from normal orientation) - forwardTransverse +// or oblique (arbitrary angle) used based on parameters - forwardOblique +// ----------------------------------------------------------------- +/** @this {import('../defs.js').ProjectionDefinition & LocalThis} */ +export function forward(p) { + return this.projectionType.forward(this, p); +} + +// inverse equations--mapping (x,y) to (lat,long) +// transverse: inverseTransverse +// oblique: inverseOblique +// ----------------------------------------------------------------- +/** @this {import('../defs.js').ProjectionDefinition & LocalThis} */ +export function inverse(p) { + return this.projectionType.inverse(this, p); +} + +/** + * @param {import('../defs.js').ProjectionDefinition & LocalThis} params - Initialized projection definition + * @param {Object} how - Transformation method + * @returns {{phip: number, lamp: number}} + */ +function createRotation(params, how) { + let phip, lamp; + if (how === paramSets.ROTATE) { + let lamc = params.oLongC; + let phic = params.oLatC; + let alpha = params.oAlpha; + if (Math.abs(Math.abs(phic) - HALF_PI) <= EPSLN) { + throw new Error('Invalid value for o_lat_c: ' + params.o_lat_c + ' should be < 90°'); + } + lamp = lamc + Math.atan2(-1 * Math.cos(alpha), -1 * Math.sin(alpha) * Math.sin(phic)); + phip = Math.asin(Math.cos(phic) * Math.sin(alpha)); + } else if (how === paramSets.NEW_POLE) { + lamp = params.oLongP; + phip = params.oLatP; + } else { + let lam1 = params.oLong1; + let phi1 = params.oLat1; + let lam2 = params.oLong2; + let phi2 = params.oLat2; + let con = Math.abs(phi1); + + if (Math.abs(phi1) > HALF_PI - EPSLN) { + throw new Error('Invalid value for o_lat_1: ' + params.o_lat_1 + ' should be < 90°'); + } + + if (Math.abs(phi2) > HALF_PI - EPSLN) { + throw new Error('Invalid value for o_lat_2: ' + params.o_lat_2 + ' should be < 90°'); + } + + if (Math.abs(phi1 - phi2) < EPSLN) { + throw new Error('Invalid value for o_lat_1 and o_lat_2: o_lat_1 should be different from o_lat_2'); + } + if (con < EPSLN) { + throw new Error('Invalid value for o_lat_1: o_lat_1 should be different from zero'); + } + + lamp = Math.atan2( + (Math.cos(phi1) * Math.sin(phi2) * Math.cos(lam1)) + - (Math.sin(phi1) * Math.cos(phi2) * Math.cos(lam2)), + (Math.sin(phi1) * Math.cos(phi2) * Math.sin(lam2)) + - (Math.cos(phi1) * Math.sin(phi2) * Math.sin(lam1)) + ); + + phip = Math.atan(-1 * Math.cos(lamp - lam1) / Math.tan(phi1)); + } + + return { lamp, phip }; +} + +/** + * Forward (lng, lat) to (x, y) for oblique case + * @param {import('../defs.js').ProjectionDefinition & LocalThis} self + * @param {{x: number, y: number}} lp - lambda, phi + */ +function forwardOblique(self, lp) { + let { x: lam, y: phi } = lp; + lam += self.long0; + const coslam = Math.cos(lam); + const sinphi = Math.sin(phi); + const cosphi = Math.cos(phi); + + lp.x = adjust_lon( + Math.atan2( + cosphi * Math.sin(lam), + (self.sphip * cosphi * coslam) + (self.cphip * sinphi) + ) + self.lamp + ); + lp.y = Math.asin( + (self.sphip * sinphi) - (self.cphip * cosphi * coslam) + ); + + return self.obliqueProjection.forward(lp); +} + +/** + * Forward (lng, lat) to (x, y) for transverse case + * @param {import('../defs.js').ProjectionDefinition & LocalThis} self + * @param {{x: number, y: number}} lp - lambda, phi + */ +function forwardTransverse(self, lp) { + let { x: lam, y: phi } = lp; + lam += self.long0; + const cosphi = Math.cos(phi); + const coslam = Math.cos(lam); + lp.x = adjust_lon( + Math.atan2( + cosphi * Math.sin(lam), + Math.sin(phi) + ) + self.lamp + ); + lp.y = Math.asin(-1 * cosphi * coslam); + + return self.obliqueProjection.forward(lp); +} + +/** + * Inverse (x, y) to (lng, lat) for oblique case + * @param {import('../defs.js').ProjectionDefinition & LocalThis} self + * @param {{x: number, y: number}} lp - lambda, phi + */ +function inverseOblique(self, lp) { + const innerLp = self.obliqueProjection.inverse(lp); + let { x: lam, y: phi } = innerLp; + + if (lam < Number.MAX_VALUE) { + lam -= self.lamp; + const coslam = Math.cos(lam); + const sinphi = Math.sin(phi); + const cosphi = Math.cos(phi); + lp.x = Math.atan2( + cosphi * Math.sin(lam), + (self.sphip * cosphi * coslam) - (self.cphip * sinphi) + ); + lp.y = Math.asin( + (self.sphip * sinphi) + (self.cphip * cosphi * coslam) + ); + } + + lp.x = adjust_lon(lp.x + self.long0); + return lp; +} + +/** + * Inverse (x, y) to (lng, lat) for transverse case + * @param {import('../defs.js').ProjectionDefinition & LocalThis} self + * @param {{x: number, y: number}} lp - lambda, phi + */ +function inverseTransverse(self, lp) { + const innerLp = self.obliqueProjection.inverse(lp); + let { x: lam, y: phi } = innerLp; + + if (lam < Number.MAX_VALUE) { + const cosphi = Math.cos(phi); + lam -= self.lamp; + lp.x = Math.atan2( + cosphi * Math.sin(lam), + -1 * Math.sin(phi) + ); + lp.y = Math.asin( + cosphi * Math.cos(lam) + ); + } + + lp.x = adjust_lon(lp.x + self.long0); + return lp; +} + +export var names = ['General Oblique Transformation', 'General_Oblique_Transformation', 'ob_tran']; +export default { + init: init, + forward: forward, + inverse: inverse, + names: names +}; diff --git a/projs.js b/projs.js index c685cdc..50b0e9c 100644 --- a/projs.js +++ b/projs.js @@ -29,6 +29,7 @@ import tpers from './lib/projections/tpers'; import geos from './lib/projections/geos'; import eqearth from './lib/projections/eqearth'; import bonne from './lib/projections/bonne'; +import ob_tran from './lib/projections/ob_tran'; export default function (proj4) { proj4.Proj.projections.add(tmerc); proj4.Proj.projections.add(etmerc); @@ -61,4 +62,5 @@ export default function (proj4) { proj4.Proj.projections.add(geos); proj4.Proj.projections.add(eqearth); proj4.Proj.projections.add(bonne); + proj4.Proj.projections.add(ob_tran); } diff --git a/test/test.js b/test/test.js index 8181d21..14236f0 100644 --- a/test/test.js +++ b/test/test.js @@ -773,6 +773,120 @@ function startTests(chai, proj4, fromArrayBuffer, testPoints) { assert.closeTo(resultWithOver[0], resultWithoutOver[0], 1e-6, 'x is close'); }); }); + + describe('ob_tran projection', function () { + it('should error no +o_proj supplied', function () { + assert.throws(function () { + proj4('+proj=ob_tran'); + }, 'Missing parameter: o_proj', 'should work'); + }); + + it('should error if +o_proj is also ob_tran', function () { + assert.throws(function () { + proj4('+proj=ob_tran +o_proj=ob_tran'); + }, 'Invalid value for o_proj: ob_tran', 'should work'); + }); + + it('should error if +o_proj is not a valid projection name', function () { + assert.throws(function () { + proj4('+proj=ob_tran +o_proj=something_that_will_never_exist'); + }, 'Could not get projection name from: +proj=something_that_will_never_exist', 'should work'); + }); + + it('should error if +o_lat_p is provided, but not +o_lon_p', function () { + assert.throws(function () { + proj4('+proj=ob_tran +o_proj=moll +o_lat_p=30'); + }, 'Missing parameter: o_lon_p', 'should work'); + }); + + it('should error if +o_lon_p is provided, but not +o_lat_p', function () { + assert.throws(function () { + proj4('+proj=ob_tran +o_proj=moll +o_lon_p=30'); + }, 'Missing parameter: o_lat_p.', 'should work'); + }); + + it('should error if parameters are NaN', function () { + [ + ['+o_lat_p=30 +o_lon_p=foo', 'o_lon_p: foo'], + ['+o_lat_p=foo +o_lon_p=30', 'o_lat_p: foo'], + ['+o_lon_p=30 +o_lat_p=foo', 'o_lat_p: foo'], + ['+o_lon_p=foo +o_lat_p=30', 'o_lon_p: foo'], + ['+o_alpha=1.0.0 +o_lon_c=0 +o_lat_c=foo', 'o_lat_c: foo'], + ['+o_lon_1=-0- +o_lon_2=1.1e2 +o_lat_1=-10 +o_lat_2=FF2', 'o_lat_2: FF2'] + ].forEach(function ([params, errorMsg]) { + assert.throws(function () { + proj4('+proj=ob_tran +o_proj=moll +' + params); + }, 'Invalid value for ' + errorMsg, 'should work'); + }); + }); + + it('shouldn\'t matter the order in which the user inputs parameters', function () { + const tests = [ + [ + 'o_lon_p', 'x_0', 'o_lat_p' + ], + [ + 'lon_0', 'y_0', 'o_lon_1', 'o_lat_2', 'o_lat_1', 'o_lon_2' + ], + [ + 'o_lon_1', 'o_lon_c', 'o_alpha', 'o_lat_c' + ] + ]; + tests.forEach(function (test) { + const str = '+proj=ob_tran +o_proj=moll ' + + test.map(function (p, i) { return '+' + p + '=' + (30 + i); }) + .join(' '); + + assert.doesNotThrow(function () { + proj4(str); + }, undefined, undefined, 'should work'); + }); + }); + + it('should throw an error if o_lat_c is 90', function () { + const tests = [90, -90]; + tests.forEach(function (test) { + const str = '+proj=ob_tran +o_proj=moll +o_lat_c=' + test + ' +o_lon_c=0 +o_alpha=1'; + assert.throws(function () { + proj4(str); + }, 'Invalid value for o_lat_c: ' + test + ' should be < 90°', 'should work'); + }); + }); + + it('should throw an error if o_lat_1 is 90', function () { + const tests = [90, -90]; + tests.forEach(function (test) { + const str = '+proj=ob_tran +o_proj=moll +o_lat_1=' + test + ' +o_lon_1=0 +o_lat_2=45 +o_lon_2=90'; + assert.throws(function () { + proj4(str); + }, 'Invalid value for o_lat_1: ' + test + ' should be < 90°', 'should work'); + }); + }); + + it('should throw an error if o_lat_2 is 90', function () { + const tests = [90, -90]; + tests.forEach(function (test) { + const str = '+proj=ob_tran +o_proj=moll +o_lat_2=' + test + ' +o_lon_1=0 +o_lat_1=45 +o_lon_2=90'; + assert.throws(function () { + proj4(str); + }, 'Invalid value for o_lat_2: ' + test + ' should be < 90°', 'should work'); + }); + }); + + it('should throw an error if o_lat_1 = o_lat_2', function () { + const str = '+proj=ob_tran +o_proj=moll +o_lat_2=45 +o_lon_1=0 +o_lat_1=45 +o_lon_2=90'; + assert.throws(function () { + proj4(str); + }, 'Invalid value for o_lat_1 and o_lat_2: o_lat_1 should be different from o_lat_2', 'should work'); + }); + + it('should throw an error if o_lat_1 is 0', function () { + const str = '+proj=ob_tran +o_proj=moll +o_lat_2=45 +o_lon_1=0 +o_lat_1=0 +o_lon_2=90'; + assert.throws(function () { + proj4(str); + }, 'Invalid value for o_lat_1: o_lat_1 should be different from zero', 'should work'); + }); + }); }); } if (typeof module !== 'undefined') { diff --git a/test/testData.js b/test/testData.js index 6379eb7..09b7722 100644 --- a/test/testData.js +++ b/test/testData.js @@ -699,6 +699,87 @@ var testPoints = [ ll: [50.977303830208, 30.915260093747], xy: [5112279.911077, -4143196.76625] }, + { + code: '+proj=moll +lon_0=10 +R=6400000', + ll: [-90, 85], + xy: [-2080158.10954, 8855229.1452], + acc: { + ll: 3, + xy: 3 + } + }, + { + code: '+proj=ob_tran +o_proj=moll +o_lat_p=45 +o_lon_p=-90', + ll: [-2, -1], + xy: [-7421459.0847, -5444548.62239], + acc: { + ll: 5, + xy: 0 + } + }, + { + code: '+proj=ob_tran +o_proj=eqearth +o_lat_p=85 +o_lon_p=10', + ll: [20, 11], + xy: [2841069.7339, 808313.2811], + acc: { + ll: 1, + xy: -4 + } + }, + { + code: '+proj=ob_tran +o_proj=moll +o_alpha=5 +o_lon_c=40 +o_lat_c=-10', + ll: [10, 5], + xy: [-154995.9625, -8241537.7451], + acc: { + ll: 3, + xy: 3 + } + }, + { + code: '+proj=ob_tran +o_proj=moll +o_lon_1=-180 +o_lon_2=180 +o_lat_1=-3 +o_lat_2=3', + ll: [10, 5], + xy: [-938419.6738, -8448989.1020], + acc: { + ll: 3, + xy: 3 + } + }, + { + code: '+proj=ob_tran +o_proj=moll +o_lon_1=-11 +o_lon_2=6 +o_lat_1=-3 +o_lat_2=3 +x_0=10000 +y_0=50000', + ll: [-90, 85], + xy: [3725830.5914, -7713738.5789], + acc: { + ll: 3, + xy: 3 + } + }, + { + code: '+proj=ob_tran +o_proj=moll +o_lon_1=-11 +o_lon_2=6 +o_lat_1=-3 +o_lat_2=3 +x_0=10000 +y_0=50000 +R=6400000', + ll: [-90, 85], + xy: [3738567.7284, -7740351.1488], + acc: { + ll: 3, + xy: 3 + } + }, + { + code: '+proj=ob_tran +o_proj=moll +R=6378137.0 +o_lon_p=0 +o_lat_p=0 +lon_0=180', + ll: [10, 20], + xy: [-1384841.18787, 7581707.88240], + acc: { + ll: 3, + xy: 3 + } + }, + { + code: '+proj=ob_tran +o_proj=moll +o_lon_p=0 +o_lat_p=0 +lon_0=180 +R=6400000 +ellps=clrk80ign +pm=paris', + ll: [10, 20], + xy: [-1068593.9375, 7685891.0261], + acc: { + ll: 3, + xy: 3 + } + }, { code: 'PROJCS["Mount Dillon / Tobago Grid",GEOGCS["Mount Dillon",DATUM["Mount_Dillon",SPHEROID["Clarke 1858",6378293.645208759,294.2606763692654,AUTHORITY["EPSG","7007"]],AUTHORITY["EPSG","6157"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4157"]],UNIT["Clarke\'s link",0.201166195164,AUTHORITY["EPSG","9039"]],PROJECTION["Cassini_Soldner"],PARAMETER["latitude_of_origin",11.25217861111111],PARAMETER["central_meridian",-60.68600888888889],PARAMETER["false_easting",187500],PARAMETER["false_northing",180000],AUTHORITY["EPSG","2066"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]', ll: [-60.676753018, 11.2487234308],