Merge pull request #531 from potion-cellar/ob-tran

feat: add ob_tran projection
This commit is contained in:
Andreas Hocevar 2025-11-08 18:57:36 +01:00 committed by GitHub
commit 33f8709947
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
10 changed files with 561 additions and 3 deletions

View File

@ -29,7 +29,8 @@ var projs = [
'tpers',
'geos',
'eqearth',
'bonne'
'bonne',
'ob_tran'
];
module.exports = function (grunt) {
grunt.initConfig({

View File

@ -38,6 +38,7 @@ import wkt from 'wkt-parser';
* @property {number} [rectified_grid_angle]
* @property {boolean} [approx]
* @property {boolean} [over]
* @property {string} [projStr]
* @property {<T extends import('./core').TemplateCoordinates>(coordinates: T, enforceAxis?: boolean) => T} inverse
* @property {<T extends import('./core').TemplateCoordinates>(coordinates: T, enforceAxis?: boolean) => T} forward
*/

View File

@ -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) {

View File

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

View File

@ -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) {

View File

@ -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) {

347
lib/projections/ob_tran.js Normal file
View File

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

View File

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

View File

@ -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') {

View File

@ -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],