ashatora/openlayers
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Use rigorous Swiss grid / EPSG:4326 transforms

Browse Source
This commit is contained in:
Tom Payne
2013-06-12 16:01:31 +02:00
parent 826556775c
commit 941432af75
2 changed files with 247 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files

212
src/ol/proj/chprojection.js
View File

@@ -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.<number>} 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.<number>} input Input array of coordinate values.
* @param {Array.<number>=} opt_output Output array of coordinate values.
* @param {number=} opt_dimension Dimension (default is 2).
* @private
* @return {Array.<number>} 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.<number>} 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.<number>} input Input array of coordinate values.
* @param {Array.<number>=} opt_output Output array of coordinate values.
* @param {number=} opt_dimension Dimension (default is 2).
* @private
* @return {Array.<number>} 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));
};