diff --git a/src/ol/proj/chprojection.js b/src/ol/proj/chprojection.js index 28fccc6614..3302674635 100644 --- a/src/ol/proj/chprojection.js +++ b/src/ol/proj/chprojection.js @@ -3,8 +3,10 @@ goog.provide('ol.proj.EPSG2056'); goog.provide('ol.proj.EPSG21781'); goog.require('goog.asserts'); +goog.require('goog.math'); goog.require('ol.Projection'); goog.require('ol.ProjectionUnits'); +goog.require('ol.ellipsoid.BESSEL1841'); goog.require('ol.proj'); goog.require('ol.proj.EPSG4326'); @@ -27,6 +29,92 @@ ol.proj.CH = function(options) { goog.inherits(ol.proj.CH, ol.Projection); +/** + * @const + * @type {number} + */ +ol.proj.CH.PHI0 = goog.math.toRadians((3600 * 46 + 60 * 57 + 8.66) / 3600); + + +/** + * @const + * @type {number} + */ +ol.proj.CH.LAMBDA0 = goog.math.toRadians((3600 * 7 + 60 * 26 + 22.5) / 3600); + + +/** + * @const + * @type {ol.Ellipsoid} + */ +ol.proj.CH.ELLIPSOID = ol.ellipsoid.BESSEL1841; + + +/** + * @const + * @type {number} + */ +ol.proj.CH.COS_PHI0 = Math.cos(ol.proj.CH.PHI0); + + +/** + * @const + * @type {number} + */ +ol.proj.CH.SIN_PHI0 = Math.sin(ol.proj.CH.PHI0); + + +/** + * @const + * @type {number} + */ +ol.proj.CH.R = ol.proj.CH.ELLIPSOID.a * Math.sqrt(1 - + ol.proj.CH.ELLIPSOID.eSquared) / (1 - ol.proj.CH.ELLIPSOID.eSquared * + ol.proj.CH.SIN_PHI0 * ol.proj.CH.SIN_PHI0); + + +/** + * @const + * @type {number} + */ +ol.proj.CH.ALPHA = Math.sqrt(1 + + ol.proj.CH.ELLIPSOID.eSquared * Math.pow(ol.proj.CH.COS_PHI0, 4) / + (1 - ol.proj.CH.ELLIPSOID.eSquared)); + + +/** + * @const + * @type {number} + */ +ol.proj.CH.SIN_B0 = ol.proj.CH.SIN_PHI0 / ol.proj.CH.ALPHA; + + +/** + * @const + * @type {number} + */ +ol.proj.CH.B0 = Math.asin(ol.proj.CH.SIN_B0); + + +/** + * @const + * @type {number} + */ +ol.proj.CH.COS_B0 = Math.cos(ol.proj.CH.B0); +// FIXME should we use Math.sqrt(1 - ol.proj.CH.SIN_B0 * ol.proj.CH.SIN_B0) ? + + +/** + * @const + * @type {number} + */ +ol.proj.CH.K = Math.log(Math.tan(Math.PI / 4 + ol.proj.CH.B0 / 2)) - + ol.proj.CH.ALPHA * Math.log(Math.tan(Math.PI / 4 + ol.proj.CH.PHI0 / 2)) + + ol.proj.CH.ALPHA * ol.proj.CH.ELLIPSOID.e * Math.log( + (1 + ol.proj.CH.ELLIPSOID.e * ol.proj.CH.SIN_PHI0) / + (1 - ol.proj.CH.ELLIPSOID.e * ol.proj.CH.SIN_PHI0)) / 2; + + /** * Add EPSG:2056 and EPSG:21781 projections, and transformations between them. */ @@ -57,7 +145,7 @@ ol.proj.CH.add = function() { * @private * @return {Array.} Output array of coordinate values. */ -ol.proj.CH.fromEPSG4326_ = +ol.proj.CH.fromEPSG4326Approximate_ = function(offsetY, offsetX, input, opt_output, opt_dimension) { var n = input.length; var dimension = goog.isDef(opt_dimension) ? opt_dimension : 2; @@ -93,6 +181,57 @@ ol.proj.CH.fromEPSG4326_ = }; +/** + * Transformation from EPSG:4326 to EPSG:2056/EPSG:21781. + * + * @see http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/projections.html + * + * @param {number} offsetY Y offset. + * @param {number} offsetX X offset. + * @param {Array.} input Input array of coordinate values. + * @param {Array.=} opt_output Output array of coordinate values. + * @param {number=} opt_dimension Dimension (default is 2). + * @private + * @return {Array.} Output array of coordinate values. + */ +ol.proj.CH.fromEPSG4326Rigorous_ = + function(offsetY, offsetX, input, opt_output, opt_dimension) { + var n = input.length; + var dimension = goog.isDef(opt_dimension) ? opt_dimension : 2; + var output; + if (goog.isDef(opt_output)) { + output = opt_output; + } else { + if (dimension > 2) { + output = input.slice(); + } else { + output = new Array(n); + } + } + goog.asserts.assert(dimension >= 2); + goog.asserts.assert(output.length % dimension === 0); + var b, bBar, eSinPhi, i, l, lambda, lBar, phi, s; + for (i = 0; i < n; i += dimension) { + lambda = goog.math.toRadians(input[i]); + phi = goog.math.toRadians(input[i + 1]); + eSinPhi = ol.proj.CH.ELLIPSOID.e * Math.sin(phi); + s = ol.proj.CH.ALPHA * Math.log(Math.tan(Math.PI / 4 + phi / 2)) - + ol.proj.CH.ALPHA * ol.proj.CH.ELLIPSOID.e * Math.log( + (1 + eSinPhi) / (1 - eSinPhi)) / 2 + ol.proj.CH.K; + b = 2 * (Math.atan(Math.exp(s)) - Math.PI / 4); + l = ol.proj.CH.ALPHA * (lambda - ol.proj.CH.LAMBDA0); + lBar = Math.atan2(Math.sin(l), + ol.proj.CH.SIN_B0 * Math.tan(b) + ol.proj.CH.COS_B0 * Math.cos(l)); + bBar = Math.asin(ol.proj.CH.COS_B0 * Math.sin(b) - + ol.proj.CH.SIN_B0 * Math.cos(b) * Math.cos(l)); + output[i] = offsetY + ol.proj.CH.R * lBar; + output[i + 1] = offsetX + ol.proj.CH.R * + Math.log((1 + Math.sin(bBar)) / (1 - Math.sin(bBar))) / 2; + } + return output; +}; + + /** * Transformation from EPSG:2056/EPSG:21781 to EPSG:4326. * @@ -108,7 +247,7 @@ ol.proj.CH.fromEPSG4326_ = * @private * @return {Array.} Output array of coordinate values. */ -ol.proj.CH.toEPSG4326_ = +ol.proj.CH.toEPSG4326Approximate_ = function(offsetY, offsetX, input, opt_output, opt_dimension) { var n = input.length; var dimension = goog.isDef(opt_dimension) ? opt_dimension : 2; @@ -144,6 +283,67 @@ ol.proj.CH.toEPSG4326_ = }; +/** + * Transformation from EPSG:2056/EPSG:21781 to EPSG:4326. + * + * @see http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/projections.html + * + * @param {number} offsetY Y offset. + * @param {number} offsetX X offset. + * @param {Array.} input Input array of coordinate values. + * @param {Array.=} opt_output Output array of coordinate values. + * @param {number=} opt_dimension Dimension (default is 2). + * @private + * @return {Array.} Output array of coordinate values. + */ +ol.proj.CH.toEPSG4326Rigorous_ = + function(offsetY, offsetX, input, opt_output, opt_dimension) { + var n = input.length; + var dimension = goog.isDef(opt_dimension) ? opt_dimension : 2; + var output; + if (goog.isDef(opt_output)) { + output = opt_output; + } else { + if (dimension > 2) { + output = input.slice(); + } else { + output = new Array(n); + } + } + goog.asserts.assert(dimension >= 2); + goog.asserts.assert(output.length % dimension === 0); + var b, bBar, eSinPhi, i, iterations, l, lambda, lastPhi, lBar, phi, s, x, y; + for (i = 0; i < n; i += dimension) { + y = input[i] - offsetY; + x = input[i + 1] - offsetX; + lBar = y / ol.proj.CH.R; + bBar = 2 * (Math.atan(Math.exp(x / ol.proj.CH.R)) - Math.PI / 4); + b = Math.asin(ol.proj.CH.COS_B0 * Math.sin(bBar) + + ol.proj.CH.SIN_B0 * Math.cos(bBar) * Math.cos(lBar)); + l = Math.atan2(Math.sin(lBar), ol.proj.CH.COS_B0 * Math.cos(lBar) - + ol.proj.CH.SIN_B0 * Math.tan(bBar)); + lambda = ol.proj.CH.LAMBDA0 + l / ol.proj.CH.ALPHA; + lastPhi = phi = b; + // Empirically, about 18 iterations are required for 1e-7 radian accuracy + for (iterations = 20; iterations > 0; --iterations) { + s = (Math.log(Math.tan(Math.PI / 4 + b / 2)) - + ol.proj.CH.K) / ol.proj.CH.ALPHA + + ol.proj.CH.ELLIPSOID.e * Math.log(Math.tan(Math.PI / 4 + + Math.asin(ol.proj.CH.ELLIPSOID.e * Math.sin(phi)) / 2)); + phi = 2 * Math.atan(Math.exp(s)) - Math.PI / 2; + if (Math.abs(phi - lastPhi) < 1e-7) { + break; + } + lastPhi = phi; + } + goog.asserts.assert(iterations !== 0); + output[i] = goog.math.toDegrees(lambda); + output[i + 1] = goog.math.toDegrees(phi); + } + return output; +}; + + /** * Transformation between EPSG:2056 and EPSG:21781. * @@ -223,8 +423,8 @@ ol.proj.EPSG2056.add = function() { ol.proj.addEquivalentTransforms( ol.proj.EPSG4326.PROJECTIONS, [epsg2056], - goog.partial(ol.proj.CH.fromEPSG4326_, 2600000, 1200000), - goog.partial(ol.proj.CH.toEPSG4326_, 2600000, 1200000)); + goog.partial(ol.proj.CH.fromEPSG4326Rigorous_, 2600000, 1200000), + goog.partial(ol.proj.CH.toEPSG4326Rigorous_, 2600000, 1200000)); }; @@ -260,6 +460,6 @@ ol.proj.EPSG21781.add = function() { ol.proj.addEquivalentTransforms( ol.proj.EPSG4326.PROJECTIONS, [epsg21781], - goog.partial(ol.proj.CH.fromEPSG4326_, 600000, 200000), - goog.partial(ol.proj.CH.toEPSG4326_, 600000, 200000)); + goog.partial(ol.proj.CH.fromEPSG4326Rigorous_, 600000, 200000), + goog.partial(ol.proj.CH.toEPSG4326Rigorous_, 600000, 200000)); }; diff --git a/test/spec/ol/proj/chprojection.test.js b/test/spec/ol/proj/chprojection.test.js index 23f9ea1a22..a4a814c595 100644 --- a/test/spec/ol/proj/chprojection.test.js +++ b/test/spec/ol/proj/chprojection.test.js @@ -9,6 +9,17 @@ describe('ol.proj.CH', function() { ol.proj.CH.add(); }); + it('has the correct constants', function() { + expect(ol.proj.CH.ELLIPSOID.eSquared).to.roughlyEqual( + 0.006674372230614, 1e-13); + expect(ol.proj.CH.R).to.roughlyEqual(6378815.90365, 1e-5); + expect(ol.proj.CH.ALPHA).to.roughlyEqual(1.00072913843038, 1e-14); + expect(ol.proj.CH.B0).to.roughlyEqual( + goog.math.toRadians((3600 * 46 + 60 * 54 + 27.83324844) / 3600), + 1e-13); + expect(ol.proj.CH.K).to.roughlyEqual(0.0030667323772751, 1e-13); + }); + it('can transform from EPSG:2056 to EPSG:21781', function() { var output = ol.proj.transform( [2660389.515487, 1185731.630396], 'EPSG:2056', 'EPSG:21781'); @@ -41,19 +52,24 @@ describe('ol.proj.EPSG2056', function() { it('transforms from EPSG:2056 to EPSG:4326', function() { var wgs84 = ol.proj.transform( - [2660389.515487, 1185731.630396], 'EPSG:2056', 'EPSG:4326'); + [2679520.05, 1212273.44], 'EPSG:2056', 'EPSG:4326'); expect(wgs84).to.be.an(Array); expect(wgs84).to.have.length(2); - expect(wgs84[0]).to.roughlyEqual(8.23, 1e-3); - expect(wgs84[1]).to.roughlyEqual(46.82, 1e-3); + expect(wgs84[0]).to.roughlyEqual( + (3600 * 8 + 60 * 29 + 11.111272) / 3600, 1e-8); + expect(wgs84[1]).to.roughlyEqual( + (3600 * 47 + 60 * 3 + 28.956592) / 3600, 1e-8); }); it('transforms from EPSG:4326 to EPSG:2056', function() { - var ch1903 = ol.proj.transform([8.23, 46.82], 'EPSG:4326', 'EPSG:2056'); - expect(ch1903).to.be.an(Array); - expect(ch1903).to.have.length(2); - expect(ch1903[0]).to.roughlyEqual(2660389.515487, 1); - expect(ch1903[1]).to.roughlyEqual(1185731.630396, 1); + var lv95 = ol.proj.transform([ + (3600 * 8 + 60 * 29 + 11.11127154) / 3600, + (3600 * 47 + 60 * 3 + 28.95659233) / 3600 + ], 'EPSG:4326', 'EPSG:2056'); + expect(lv95).to.be.an(Array); + expect(lv95).to.have.length(2); + expect(lv95[0]).to.roughlyEqual(2679520.05, 1e-3); + expect(lv95[1]).to.roughlyEqual(1212273.44, 1e-3); }); }); @@ -69,7 +85,7 @@ describe('ol.proj.EPSG21781', function() { expect(epsg21781).to.be.an(ol.Projection); }); - it('does not lose too much accuracy when round-tripping', function() { + it('maintains accuracy when round-tripping', function() { var extent = epsg21781.getExtent(); var fromEPSG4326 = ol.proj.getTransform('EPSG:4326', 'EPSG:21781'); var toEPSG4326 = ol.proj.getTransform('EPSG:21781', 'EPSG:4326'); @@ -79,32 +95,38 @@ describe('ol.proj.EPSG21781', function() { roundTripped = fromEPSG4326(toEPSG4326([x, y])); expect(roundTripped).to.be.an(Array); expect(roundTripped).to.have.length(2); - expect(roundTripped[0]).to.roughlyEqual(x, 1e1); - expect(roundTripped[1]).to.roughlyEqual(y, 1e1); + expect(roundTripped[0]).to.roughlyEqual(x, 1e-3); + expect(roundTripped[1]).to.roughlyEqual(y, 1e-3); } } }); it('transforms from EPSG:21781 to EPSG:4326', function() { var wgs84 = ol.proj.transform( - [660389.515487, 185731.630396], 'EPSG:21781', 'EPSG:4326'); + [679520.05, 212273.44], 'EPSG:21781', 'EPSG:4326'); expect(wgs84).to.be.an(Array); expect(wgs84).to.have.length(2); - expect(wgs84[0]).to.roughlyEqual(8.23, 1e-3); - expect(wgs84[1]).to.roughlyEqual(46.82, 1e-3); + expect(wgs84[0]).to.roughlyEqual( + (3600 * 8 + 60 * 29 + 11.111272) / 3600, 1e-8); + expect(wgs84[1]).to.roughlyEqual( + (3600 * 47 + 60 * 3 + 28.956592) / 3600, 1e-8); }); it('transforms from EPSG:4326 to EPSG:21781', function() { - var ch1903 = ol.proj.transform([8.23, 46.82], 'EPSG:4326', 'EPSG:21781'); - expect(ch1903).to.be.an(Array); - expect(ch1903).to.have.length(2); - expect(ch1903[0]).to.roughlyEqual(660389.515487, 1); - expect(ch1903[1]).to.roughlyEqual(185731.630396, 1); + var lv03 = ol.proj.transform([ + (3600 * 8 + 60 * 29 + 11.11127154) / 3600, + (3600 * 47 + 60 * 3 + 28.95659233) / 3600 + ], 'EPSG:4326', 'EPSG:21781'); + expect(lv03).to.be.an(Array); + expect(lv03).to.have.length(2); + expect(lv03[0]).to.roughlyEqual(679520.05, 1e-3); + expect(lv03[1]).to.roughlyEqual(212273.44, 1e-3); }); }); +goog.require('goog.math'); goog.require('ol.Projection'); goog.require('ol.proj'); goog.require('ol.proj.CH');