From 738bee536e4d1c08bd8b9919998c2091563e510f Mon Sep 17 00:00:00 2001 From: Fu Zhen Date: Fri, 20 Oct 2023 18:13:39 +0800 Subject: [PATCH] add support for UTM and Traverse Mercator projection (#2108) * add support for UTM and Traverse Mercator projection * fix specs * fix specs * fix specs --- src/geo/projection/Projection.EPSG4326.js | 1 + src/geo/projection/Projection.EPSG4490.js | 22 -- src/geo/projection/Projection.EPSG9807.js | 89 ++++++ src/geo/projection/Projection.UTM.js | 37 +++ src/geo/projection/etmerc.js | 256 ++++++++++++++++++ src/geo/projection/index.js | 3 +- src/map/spatial-reference/SpatialReference.js | 47 +++- test/map/ProjectionSpec.js | 49 +++- 8 files changed, 467 insertions(+), 37 deletions(-) delete mode 100644 src/geo/projection/Projection.EPSG4490.js create mode 100644 src/geo/projection/Projection.EPSG9807.js create mode 100644 src/geo/projection/Projection.UTM.js create mode 100644 src/geo/projection/etmerc.js diff --git a/src/geo/projection/Projection.EPSG4326.js b/src/geo/projection/Projection.EPSG4326.js index f48ecea14..2ecbe99b8 100644 --- a/src/geo/projection/Projection.EPSG4326.js +++ b/src/geo/projection/Projection.EPSG4326.js @@ -21,6 +21,7 @@ export default extend({}, Common, /** @lends projection.EPSG4326 */ { * @constant */ code: 'EPSG:4326', + aliases: ['EPSG:4490'], project: function (p, out) { if (out) { out.x = p.x; diff --git a/src/geo/projection/Projection.EPSG4490.js b/src/geo/projection/Projection.EPSG4490.js deleted file mode 100644 index 4c8bffbfb..000000000 --- a/src/geo/projection/Projection.EPSG4490.js +++ /dev/null @@ -1,22 +0,0 @@ -import { extend } from '../../core/util'; -import PROJ4326 from './Projection.EPSG4326'; - -/** - * For CGCS2000 - * - * @class - * @category geo - * @protected - * @memberOf projection - * @name EPSG4490 - * @mixes projection.EPSG4326 - * @mixes measurer.WGS84Sphere - */ -export default extend({}, PROJ4326, /** @lends projection.EPSG4490 */ { - /** - * "EPSG:4490", Code of the projection - * @type {String} - * @constant - */ - code: 'EPSG:4490' -}); diff --git a/src/geo/projection/Projection.EPSG9807.js b/src/geo/projection/Projection.EPSG9807.js new file mode 100644 index 000000000..d57ea6f7d --- /dev/null +++ b/src/geo/projection/Projection.EPSG9807.js @@ -0,0 +1,89 @@ +import { isNil, extend } from '../../core/util'; +import Common from './Projection'; +import Coordinate from '../Coordinate'; +import { WGS84Sphere } from '../measurer'; +import etmerc from './etmerc.js'; + +// from proj_api.h +const RAD_TO_DEG = 57.295779513082321, + DEG_TO_RAD = 0.017453292519943296; + +// from pj_transform.c +const SRS_WGS84_SEMIMAJOR = 6378137; +const SRS_WGS84_ESQUARED = 0.0066943799901413165; + +/** + * Traverse Mercator Projection + * + * @class + * @category geo + * @protected + * @memberOf projection + * @name EPSG9807 + * @mixes projection.Common + * @mixes measurer.WGS84Sphere + */ +export default { + code: 'EPSG:9807', + aliases: ['Traverse_Mercator'], + create(params) { + const P = { + a: SRS_WGS84_SEMIMAJOR, + es: SRS_WGS84_ESQUARED, + x0: isNil(params.falseEasting) ? 500000 : params.falseEasting, + y0: isNil(params.falseNorthing) ? 0 : params.falseNorthing, + k0: params.scaleFactor || 0.9996, + lam0: (params.centralMeridian || 0) * DEG_TO_RAD, + phi0: (params.latitudeOfOrigin || 0) * DEG_TO_RAD, + originLam0: params.startLongtitude || 0, + originPhi0: params.startLatitude || 0 + }; + etmerc(P); + const lp = { lam: 0, phi: 0 }; + const xy = {}; + let originX = 0; + let originY = 0; + if (P.originLam0 || P.originPhi0) { + lp.lam = P.originLam0 * DEG_TO_RAD - P.lam0; + lp.phi = P.originPhi0 * DEG_TO_RAD; + P.fwd(lp, xy); + originX = P.a * xy.x + P.x0; + originY = P.a * xy.y + P.y0; + } + + return extend({}, Common, { + /** + * "EPSG:9807", Code of the projection + * @type {String} + * @constant + */ + code: 'EPSG:9807', + project: function (p, out) { + lp.lam = p.x * DEG_TO_RAD - P.lam0; + lp.phi = p.y * DEG_TO_RAD; + P.fwd(lp, xy); + const x = P.a * xy.x + P.x0 - originX; + const y = P.a * xy.y + P.y0 - originY; + if (out) { + out.x = x; + out.y = y; + return out; + } + return new Coordinate(x, y); + }, + unproject: function (p, out) { + xy.x = (p.x - P.x0 + originX) / P.a; + xy.y = (p.y - P.y0 + originY) / P.a; + P.inv(xy, lp); + const x = (lp.lam + P.lam0) * RAD_TO_DEG; + const y = lp.phi * RAD_TO_DEG; + if (out) { + out.x = x; + out.y = y; + return out; + } + return new Coordinate(x, y); + } + }, WGS84Sphere); + } +}; diff --git a/src/geo/projection/Projection.UTM.js b/src/geo/projection/Projection.UTM.js new file mode 100644 index 000000000..c13f86d1f --- /dev/null +++ b/src/geo/projection/Projection.UTM.js @@ -0,0 +1,37 @@ +import { extend } from '../../core/util'; +import PROJ9807 from './Projection.EPSG9807'; + +/** + * Universal Traverse Mercator projection + * + * @class + * @category geo + * @protected + * @memberOf projection + * @name EPSG4490 + * @mixes projection.EPSG4326 + * @mixes measurer.WGS84Sphere + */ +export default extend({}, PROJ9807, /** @lends projection.EPSG4490 */ { + /** + * "EPSG:4490", Code of the projection + * @type {String} + * @constant + */ + code: 'utm', + aliases: [], + create(params) { + const P = {}; + let zone = parseInt(params.zone); + P.falseNorthing = params.south ? 10000000 : 0; + P.falseEasting = 500000; + if (zone > 0 && zone <= 60) { + zone--; + } else { + throw new Error('zone must be > 0 and <= 60.'); + } + P.centralMeridian = (zone + 0.5) * 6 - 180; + P.scaleFactor = 0.9996; + return PROJ9807.create(P); + } +}); diff --git a/src/geo/projection/etmerc.js b/src/geo/projection/etmerc.js new file mode 100644 index 000000000..ca0d8839a --- /dev/null +++ b/src/geo/projection/etmerc.js @@ -0,0 +1,256 @@ +/* eslint-disable */ +// https://github.com/mbloch/mapshaper-proj/blob/master/src/projections/etmerc.js + +// add math.h functions to library scope +// (to make porting projection functions simpler) +var fabs = Math.abs, + floor = Math.floor, + sin = Math.sin, + cos = Math.cos, + tan = Math.tan, + asin = Math.asin, + acos = Math.acos, + atan = Math.atan, + atan2 = Math.atan2, + sqrt = Math.sqrt, + pow = Math.pow, + exp = Math.exp, + log = Math.log, + hypot = Math.hypot, + sinh = Math.sinh, + cosh = Math.cosh, + MIN = Math.min, + MAX = Math.max; + +// constants from math.h +var HUGE_VAL = Infinity, + M_PI = Math.PI; + +// from proj_api.h +var RAD_TO_DEG = 57.295779513082321, + DEG_TO_RAD = 0.017453292519943296; + +// from pj_transform.c +var SRS_WGS84_SEMIMAJOR = 6378137; +var SRS_WGS84_ESQUARED = 0.0066943799901413165; + +// math constants from project.h +var M_FORTPI = M_PI / 4, + M_HALFPI = M_PI / 2, + M_PI_HALFPI = 1.5 * M_PI, + M_TWOPI = 2 * M_PI, + M_TWO_D_PI = 2 / M_PI, + M_TWOPI_HALFPI = 2.5 * M_PI; + +// datum types +var PJD_UNKNOWN = 0, + PJD_3PARAM = 1, + PJD_7PARAM = 2, + PJD_GRIDSHIFT = 3, + PJD_WGS84 = 4; + +// named errors +var PJD_ERR_GEOCENTRIC = -45, + PJD_ERR_AXIS = -47, + PJD_ERR_GRID_AREA = -48, + PJD_ERR_CATALOG = -49; + +// common +var EPS10 = 1e-10; + +function e_error() {} + +function pj_etmerc(P) { + var cgb = [], + cbg = [], + utg = [], + gtu = [], + Qn, Zb, f, n, np, Z; + if (P.es <= 0) e_error(-34); + /* flattening */ + f = P.es / (1 + sqrt(1 - P.es)); /* Replaces: f = 1 - sqrt(1-P.es); */ + /* third flattening */ + np = n = f/(2 - f); + /* COEF. OF TRIG SERIES GEO <-> GAUSS */ + /* cgb := Gaussian -> Geodetic, KW p190 - 191 (61) - (62) */ + /* cbg := Geodetic -> Gaussian, KW p186 - 187 (51) - (52) */ + /* PROJ_ETMERC_ORDER = 6th degree : Engsager and Poder: ICC2007 */ + cgb[0] = n*(2 + n*(-2/3 + n * (-2 + n*(116/45 + n * (26/45 + n*(-2854/675 )))))); + cbg[0] = n*(-2 + n*( 2/3 + n*( 4/3 + n*(-82/45 + n*(32/45 + n*(4642/4725)))))); + np *= n; + cgb[1] = np*(7/3 + n*(-8/5 + n*(-227/45 + n*(2704/315 + n*(2323/945))))); + cbg[1] = np*(5/3 + n*(-16/15 + n*( -13/9 + n*(904/315 + n*(-1522/945))))); + np *= n; + /* n^5 coeff corrected from 1262/105 -> -1262/105 */ + cgb[2] = np*(56/15 + n*(-136/35 + n*(-1262/105 + n*(73814/2835)))); + cbg[2] = np*(-26/15 + n*(34/21 + n*(8/5 + n*(-12686/2835)))); + np *= n; + /* n^5 coeff corrected from 322/35 -> 332/35 */ + cgb[3] = np*(4279/630 + n*(-332/35 + n*(-399572/14175))); + cbg[3] = np*(1237/630 + n*(-12/5 + n*( -24832/14175))); + np *= n; + cgb[4] = np*(4174/315 + n*(-144838/6237)); + cbg[4] = np*(-734/315 + n*(109598/31185)); + np *= n; + cgb[5] = np*(601676/22275); + cbg[5] = np*(444337/155925); + + /* Constants of the projections */ + /* Transverse Mercator (UTM, ITM, etc) */ + np = n*n; + /* Norm. mer. quad, K&W p.50 (96), p.19 (38b), p.5 (2) */ + Qn = P.k0/(1 + n) * (1 + np*(1/4 + np*(1/64 + np/256))); + /* coef of trig series */ + /* utg := ell. N, E -> sph. N, E, KW p194 (65) */ + /* gtu := sph. N, E -> ell. N, E, KW p196 (69) */ + utg[0] = n*(-0.5 + n*( 2/3 + n*(-37/96 + n*( 1/360 + n*(81/512 + n*(-96199/604800)))))); + gtu[0] = n*(0.5 + n*(-2/3 + n*(5/16 + n*(41/180 + n*(-127/288 + n*(7891/37800)))))); + utg[1] = np*(-1/48 + n*(-1/15 + n*(437/1440 + n*(-46/105 + n*(1118711/3870720))))); + gtu[1] = np*(13/48 + n*(-3/5 + n*(557/1440 + n*(281/630 + n*(-1983433/1935360))))); + np *= n; + utg[2] = np*(-17/480 + n*(37/840 + n*(209/4480 + n*(-5569/90720 )))); + gtu[2] = np*(61/240 + n*(-103/140 + n*(15061/26880 + n*(167603/181440)))); + np *= n; + utg[3] = np*(-4397/161280 + n*(11/504 + n*(830251/7257600))); + gtu[3] = np*(49561/161280 + n*(-179/168 + n*(6601661/7257600))); + np *= n; + utg[4] = np*(-4583/161280 + n*(108847/3991680)); + gtu[4] = np*(34729/80640 + n*(-3418889/1995840)); + np *= n; + utg[5] = np*(-20648693/638668800); + gtu[5] = np*(212378941/319334400); + + /* Gaussian latitude value of the origin latitude */ + Z = gatg(cbg, P.phi0); + + /* Origin northing minus true northing at the origin latitude */ + /* i.e. true northing = N - P.Zb */ + Zb = -Qn*(Z + clens(gtu, 2*Z)); + P.fwd = e_fwd; + P.inv = e_inv; + + function e_fwd(lp, xy) { + var sin_Cn, cos_Cn, cos_Ce, sin_Ce, tmp; + var Cn = lp.phi, Ce = lp.lam; + + /* ell. LAT, LNG -> Gaussian LAT, LNG */ + Cn = gatg(cbg, Cn); + /* Gaussian LAT, LNG -> compl. sph. LAT */ + sin_Cn = sin(Cn); + cos_Cn = cos(Cn); + sin_Ce = sin(Ce); + cos_Ce = cos(Ce); + Cn = atan2(sin_Cn, cos_Ce*cos_Cn); + Ce = atan2(sin_Ce*cos_Cn, hypot(sin_Cn, cos_Cn*cos_Ce)); + /* compl. sph. N, E -> ell. norm. N, E */ + Ce = asinhy(tan(Ce)); + tmp = clenS(gtu, 2*Cn, 2*Ce); + Cn += tmp[0]; + Ce += tmp[1]; + if (fabs (Ce) <= 2.623395162778) { + xy.y = Qn * Cn + Zb; /* Northing */ + xy.x = Qn * Ce; /* Easting */ + } else { + xy.x = xy.y = HUGE_VAL; + } + } + + function e_inv(xy, lp) { + var sin_Cn, cos_Cn, cos_Ce, sin_Ce, tmp; + var Cn = xy.y, Ce = xy.x; + /* normalize N, E */ + Cn = (Cn - Zb)/Qn; + Ce = Ce/Qn; + if (fabs(Ce) <= 2.623395162778) { /* 150 degrees */ + /* norm. N, E -> compl. sph. LAT, LNG */ + tmp = clenS(utg, 2*Cn, 2*Ce); + Cn += tmp[0]; + Ce += tmp[1]; + Ce = atan(sinh(Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - M_FORTPI); */ + /* compl. sph. LAT -> Gaussian LAT, LNG */ + sin_Cn = sin(Cn); + cos_Cn = cos(Cn); + sin_Ce = sin(Ce); + cos_Ce = cos(Ce); + Ce = atan2(sin_Ce, cos_Ce*cos_Cn); + Cn = atan2(sin_Cn*cos_Ce, hypot(sin_Ce, cos_Ce*cos_Cn)); + /* Gaussian LAT, LNG -> ell. LAT, LNG */ + lp.phi = gatg (cgb, Cn); + lp.lam = Ce; + } + else { + lp.phi = lp.lam = HUGE_VAL; + } + } + + function log1py(x) { + var y = 1 + x, + z = y - 1; + return z === 0 ? x : x * log(y) / z; + } + + function asinhy(x) { + var y = fabs(x); + y = log1py(y * (1 + y/(hypot(1, y) + 1))); + return x < 0 ? -y : y; + } + + function gatg(pp, B) { + var cos_2B = 2 * cos(2 * B), + i = pp.length - 1, + h1 = pp[i], + h2 = 0, + h; + while (--i >= 0) { + h = -h2 + cos_2B * h1 + pp[i]; + h2 = h1; + h1 = h; + } + return (B + h * sin(2 * B)); + } + + function clens(pp, arg_r) { + var r = 2 * cos(arg_r), + i = pp.length - 1, + hr1 = pp[i], + hr2 = 0, + hr; + while (--i >= 0) { + hr = -hr2 + r * hr1 + pp[i]; + hr2 = hr1; + hr1 = hr; + } + return sin(arg_r) * hr; + } + + function clenS(pp, arg_r, arg_i) { + var sin_arg_r = sin(arg_r), + cos_arg_r = cos(arg_r), + sinh_arg_i = sinh(arg_i), + cosh_arg_i = cosh(arg_i), + r = 2 * cos_arg_r * cosh_arg_i, + i = -2 * sin_arg_r * sinh_arg_i, + j = pp.length - 1, + hr = pp[j], + hi1 = 0, + hr1 = 0, + hi = 0, + hr2, hi2; + while (--j >= 0) { + hr2 = hr1; + hi2 = hi1; + hr1 = hr; + hi1 = hi; + hr = -hr2 + r*hr1 - i * hi1 + pp[j]; + hi = -hi2 + i*hr1 + r * hi1; + } + r = sin_arg_r * cosh_arg_i; + i = cos_arg_r * sinh_arg_i; + return [r * hr - i * hi, r * hi + i * hr]; + } +} + + +export default pj_etmerc; +/* eslint-enable */ + diff --git a/src/geo/projection/index.js b/src/geo/projection/index.js index fca08b410..bc0c42329 100644 --- a/src/geo/projection/index.js +++ b/src/geo/projection/index.js @@ -3,7 +3,8 @@ import EPSG3857 from './Projection.EPSG3857'; export { default as EPSG4326 } from './Projection.EPSG4326'; -export { default as EPSG4490 } from './Projection.EPSG4490'; +export { default as EPSG9807 } from './Projection.EPSG9807'; +export { default as UTM } from './Projection.UTM'; export { default as BAIDU } from './Projection.Baidu'; export { default as IDENTITY } from './Projection.IDENTITY'; export { EPSG3857 }; diff --git a/src/map/spatial-reference/SpatialReference.js b/src/map/spatial-reference/SpatialReference.js index 2441455c1..f7dbb2aa1 100644 --- a/src/map/spatial-reference/SpatialReference.js +++ b/src/map/spatial-reference/SpatialReference.js @@ -1,4 +1,4 @@ -import { extend, isNil, isObject, hasOwn, sign, isString } from '../../core/util'; +import { extend, isNil, hasOwn, sign, isString } from '../../core/util'; import Coordinate from '../../geo/Coordinate'; import Extent from '../../geo/Extent'; import * as projections from '../../geo/projection'; @@ -121,27 +121,48 @@ export default class SpatialReference { return Object.keys(DefaultSpatialReference); } - static getProjectionInstance(prjName) { - if (!prjName) { + static getProjectionInstance(projection) { + if (!projection) { return null; } - if (isObject(prjName)) { - if (!prjName.locate) { - prjName = extend({}, prjName); - if (prjName.measure === 'identity') { - extend(prjName, Measurer.getInstance('IDENTITY')); + if (isString(projection)) { + projection = { + code: projection + }; + } + // a custom one + if (projection.project) { + if (!projection.locate) { + projection = extend({}, projection); + if (projection.measure === 'identity') { + extend(projection, Measurer.getInstance('IDENTITY')); } else { - extend(prjName, Measurer.getInstance('EPSG:4326')); + extend(projection, Measurer.getInstance('EPSG:4326')); } } - return prjName; + return projection; } - prjName = (prjName + '').toLowerCase(); + const prjName = (projection.code + '').toLowerCase(); for (const p in projections) { if (hasOwn(projections, p)) { + const names = projections[p].aliases || []; const code = projections[p]['code']; - if (code && code.toLowerCase() === prjName) { - return projections[p]; + names.push(code); + for (let i = 0; i < names.length; i++) { + if (names[i].toLowerCase() === prjName) { + if (projections[p].create) { + const instance = projections[p].create(projection); + instance.code = names[i]; + return instance; + } else { + if (projections[p].code === names[i]) { + return projections[p]; + } + const instance = extend({}, projections[p]); + instance.code = names[i]; + return instance; + } + } } } } diff --git a/test/map/ProjectionSpec.js b/test/map/ProjectionSpec.js index c1ef04f50..d49523243 100644 --- a/test/map/ProjectionSpec.js +++ b/test/map/ProjectionSpec.js @@ -205,5 +205,52 @@ describe('Map.Projection', function () { var b = baiduProj.project(new maptalks.Coordinate(-Infinity, -Infinity)); }); }); -}); + describe('utm projection', function () { + it('utm project and unproject', function () { + var spatialReference = { + projection: { + code: 'utm', + zone: 51 + }, + resolutions: (function () { + let res = Math.pow(2, 8); + const resolutions = []; + for (let i = 0; i < 23; i++) { + resolutions[i] = res; + res *= 0.5; + } + return resolutions; + })(), + fullExtent: { + 'top': 200000, + 'left': -200000, + 'bottom': -200000, + 'right': 200000 + } + }; + map.remove(); + map = new maptalks.Map(container, { + center: [121.50, 31.24], + zoom: 5, + centerCross: true, + spatialReference + }); + var center = map.getCenter(); + var projection = map.getProjection(); + var p0 = center.add(-0.018, -0.003); + var p1 = center.add(0.018, -0.003); + + var pj0 = projection.project(p0); + var pj1 = projection.project(p1); + + expect(pj0.toArray()).to.be.eql([355434.28874248517, 3456861.229232939]); + expect(pj1.toArray()).to.be.eql([358863.0849529234, 3456814.667859022]); + + var pp0 = projection.unproject(pj0); + var pp1 = projection.unproject(pj1); + expect(pp0.toArray()).to.be.eql([121.482, 31.23699999999999]); + expect(pp1.toArray()).to.be.eql([121.518, 31.237]); + }); + }); +});