diff --git a/src/ol/ellipsoid/bessel1841ellipsoid.js b/src/ol/ellipsoid/bessel1841ellipsoid.js new file mode 100644 index 0000000000..198a8eb76b --- /dev/null +++ b/src/ol/ellipsoid/bessel1841ellipsoid.js @@ -0,0 +1,10 @@ +goog.provide('ol.ellipsoid.BESSEL1841'); + +goog.require('ol.Ellipsoid'); + + +/** + * @const + * @type {ol.Ellipsoid} + */ +ol.ellipsoid.BESSEL1841 = new ol.Ellipsoid(6377397.155, 1 / 299.15281285); diff --git a/src/ol/proj/chprojection.js b/src/ol/proj/chprojection.js new file mode 100644 index 0000000000..3302674635 --- /dev/null +++ b/src/ol/proj/chprojection.js @@ -0,0 +1,465 @@ +goog.provide('ol.proj.CH'); +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'); + + + +/** + * Internal base class for Swiss grid projections. + * @constructor + * @extends {ol.Projection} + * @param {{code: string, extent: ol.Extent}} options Options. + */ +ol.proj.CH = function(options) { + goog.base(this, { + code: options.code, + extent: options.extent, + global: false, + units: ol.ProjectionUnits.METERS + }); +}; +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. + */ +ol.proj.CH.add = function() { + ol.proj.EPSG2056.add(); + ol.proj.EPSG21781.add(); + var epsg2056 = ol.proj.get('EPSG:2056'); + var epsg21781 = ol.proj.get('EPSG:21781'); + ol.proj.addTransform(epsg2056, epsg21781, + goog.partial(ol.proj.CH.translate_, -2000000, -1000000)); + ol.proj.addTransform(epsg21781, epsg2056, + goog.partial(ol.proj.CH.translate_, 2000000, 1000000)); +}; + + +/** + * Transformation from EPSG:4326 to EPSG:2056/EPSG:21781. + * + * This uses an approximation that is accurate to about 1m. + * + * @see http://www.swisstopo.admin.ch/internet/swisstopo/en/home/products/software/products/skripts.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.fromEPSG4326Approximate_ = + 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 auxLat, auxLon, i; + for (i = 0; i < n; i += dimension) { + auxLat = 36 * input[i + 1] / 100 - 16.902866; + auxLon = 36 * input[i] / 100 - 2.67825; + output[i] = offsetY + 72.37 + + 211455.93 * auxLon - + 10938.51 * auxLon * auxLat - + 0.36 * auxLon * auxLat * auxLat - + 44.54 * auxLon * auxLon * auxLon; + output[i + 1] = offsetX + 147.07 + + 308807.95 * auxLat + + 3745.25 * auxLon * auxLon + + 76.63 * auxLat * auxLat - + 194.56 * auxLon * auxLon * auxLat + + 119.79 * auxLat * auxLat * auxLat; + } + return output; +}; + + +/** + * 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. + * + * This uses an approximation that is accurate to about 1m. + * + * @see http://www.swisstopo.admin.ch/internet/swisstopo/en/home/products/software/products/skripts.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.toEPSG4326Approximate_ = + 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 auxX, auxY, i; + for (i = 0; i < n; i += dimension) { + auxY = (input[i] - offsetY) / 1000000; + auxX = (input[i + 1] - offsetX) / 1000000; + output[i] = 100 * (2.6779094 + + 4.728982 * auxY + + 0.791484 * auxY * auxX + + 0.1306 * auxY * auxX * auxX - + 0.0436 * auxY * auxY * auxY) / 36; + output[i + 1] = 100 * (16.9023892 + + 3.238272 * auxX - + 0.270978 * auxY * auxY - + 0.002528 * auxX * auxX - + 0.0447 * auxY * auxY * auxX - + 0.014 * auxX * auxX * auxX) / 36; + } + return output; +}; + + +/** + * 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. + * + * Currently a simple offset is used. This is accurate to within 3m. + * + * @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.translate_ = + 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 i; + for (i = 0; i < n; i += dimension) { + output[i] = input[i] + offsetY; + output[i + 1] = input[i + 1] + offsetX; + } + return output; +}; + + +/** + * @inheritDoc + */ +ol.proj.CH.prototype.getPointResolution = function(resolution, point) { + return resolution; +}; + + + +/** + * The EPSG:2056 projection, also known as LV95 (CH1903+). + * @constructor + * @extends {ol.proj.CH} + */ +ol.proj.EPSG2056 = function() { + goog.base(this, { + code: 'EPSG:2056', + extent: ol.proj.EPSG2056.EXTENT + }); +}; +goog.inherits(ol.proj.EPSG2056, ol.proj.CH); + + +/** + * @const + * @type {ol.Extent} + */ +ol.proj.EPSG2056.EXTENT = + [2485869.5728, 2837076.5648, 1076443.1884, 1299941.7864]; + + +/** + * Add the EPSG:2056 projection and transformations to and from EPSG:4326. + */ +ol.proj.EPSG2056.add = function() { + ol.proj.addEquivalentProjections(ol.proj.EPSG4326.PROJECTIONS); + var epsg2056 = new ol.proj.EPSG2056(); + ol.proj.addProjection(epsg2056); + ol.proj.addEquivalentTransforms( + ol.proj.EPSG4326.PROJECTIONS, + [epsg2056], + goog.partial(ol.proj.CH.fromEPSG4326Rigorous_, 2600000, 1200000), + goog.partial(ol.proj.CH.toEPSG4326Rigorous_, 2600000, 1200000)); +}; + + + +/** + * The EPSG:21781 projection, also known as LV03 (CH1903). + * @constructor + * @extends {ol.proj.CH} + */ +ol.proj.EPSG21781 = function() { + goog.base(this, { + code: 'EPSG:21781', + extent: ol.proj.EPSG21781.EXTENT + }); +}; +goog.inherits(ol.proj.EPSG21781, ol.proj.CH); + + +/** + * @const + * @type {ol.Extent} + */ +ol.proj.EPSG21781.EXTENT = [485869.5728, 837076.5648, 76443.1884, 299941.7864]; + + +/** + * Add the EPSG:21781 projection and transformations to and from EPSG:4326. + */ +ol.proj.EPSG21781.add = function() { + ol.proj.addEquivalentProjections(ol.proj.EPSG4326.PROJECTIONS); + var epsg21781 = new ol.proj.EPSG21781(); + ol.proj.addProjection(epsg21781); + ol.proj.addEquivalentTransforms( + ol.proj.EPSG4326.PROJECTIONS, + [epsg21781], + 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 new file mode 100644 index 0000000000..a4a814c595 --- /dev/null +++ b/test/spec/ol/proj/chprojection.test.js @@ -0,0 +1,134 @@ +goog.provide('ol.test.proj.CH'); +goog.provide('ol.test.proj.EPSG2056'); +goog.provide('ol.test.proj.EPSG21781'); + + +describe('ol.proj.CH', function() { + + beforeEach(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'); + expect(output).to.be.an(Array); + expect(output).to.have.length(2); + expect(output[0]).to.roughlyEqual(660389.515487, 1e-9); + expect(output[1]).to.roughlyEqual(185731.630396, 1e-9); + }); + + it('can transform from EPSG:21781 to EPSG:2056', function() { + var output = ol.proj.transform( + [660389.515487, 185731.630396], 'EPSG:21781', 'EPSG:2056'); + expect(output).to.be.an(Array); + expect(output).to.have.length(2); + expect(output[0]).to.roughlyEqual(2660389.515487, 1e-10); + expect(output[1]).to.roughlyEqual(1185731.630396, 1e-10); + }); + +}); + + +describe('ol.proj.EPSG2056', function() { + + var epsg2056; + beforeEach(function() { + ol.proj.EPSG2056.add(); + epsg2056 = ol.proj.get('EPSG:2056'); + expect(epsg2056).to.be.an(ol.Projection); + }); + + it('transforms from EPSG:2056 to EPSG:4326', function() { + var wgs84 = ol.proj.transform( + [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( + (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 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); + }); + +}); + + + +describe('ol.proj.EPSG21781', function() { + + var epsg21781; + beforeEach(function() { + ol.proj.EPSG21781.add(); + epsg21781 = ol.proj.get('EPSG:21781'); + expect(epsg21781).to.be.an(ol.Projection); + }); + + 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'); + var roundTripped, x, y; + for (x = extent[0]; x < extent[1]; x += 50000) { + for (y = extent[2]; y < extent[3]; y += 50000) { + roundTripped = fromEPSG4326(toEPSG4326([x, y])); + expect(roundTripped).to.be.an(Array); + expect(roundTripped).to.have.length(2); + 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( + [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( + (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 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'); +goog.require('ol.proj.EPSG2056'); +goog.require('ol.proj.EPSG21781');