Last active
August 29, 2015 14:05
-
-
Save atlefren/95302a356854d1e05863 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/** | |
Simple script to transform coordinates from WGS84 Geographic to UTM Zone 33N. | |
Based on Holsen.js, which in turn is based on "Holsens småprogrammer" | |
License: MIT | |
Written by: Atle F. Sveen | |
**/ | |
var toUtm33 = (function () { | |
'use strict'; | |
//make these easier to type/read | |
var cos = Math.cos; | |
var sin = Math.sin; | |
var tan = Math.tan; | |
var pow = Math.pow; | |
var sqrt = Math.sqrt; | |
//Ellipsoid constants for WGS84-ellipsoid | |
var ellipsoid = { | |
a: 6378137.000, | |
f: 1 / 298.257223563 | |
}; | |
ellipsoid.b = -1 * ((ellipsoid.f * ellipsoid.a) - ellipsoid.a); | |
//Values specific for Zone 33N | |
var coordsystem = { | |
scaleFactor: 0.9996, | |
falseEasting: 500000, | |
falseNorthing: 0, | |
centralMeridian: 15, | |
latitudeOfOrigin: 0 | |
}; | |
function toRad(deg) { | |
return deg / (180 / Math.PI); | |
} | |
function ellipsoidParams(ellipsoid) { | |
var n = (ellipsoid.a - ellipsoid.b) / (ellipsoid.a + ellipsoid.b); | |
return { | |
n: n, | |
k: ellipsoid.a / (n + 1), | |
k1: 1 + n * (n / 4) + pow(n, 4) / 64, | |
k2: (n - n * n * n / 8) * 3, | |
k3: (n * n - pow(n, 4) / 4) * 15 / 8 | |
}; | |
} | |
function getGeodeticArch(lon1, lon2, ell) { | |
var db = lon1 - lon2; | |
var bm = lon1 - db / 2; | |
return (ell.k * (ell.k1 * db - ell.k2 * cos(2 * bm) * sin(db) + | |
ell.k3 * cos(4 * bm) * sin(2 * db)) - | |
ell.k * (pow(ell.n, 3) * cos(6 * bm) * sin(3 * db) * 35 / 24) + | |
ell.k * (pow(ell.n, 4) * cos(8 * bm) * sin(4 * db) * 315) / 256); | |
} | |
return function (lon, lat) { | |
lon = toRad(lon) - toRad(coordsystem.centralMeridian); | |
lat = toRad(lat); | |
var a = ellipsoid.a; | |
var b = ellipsoid.b; | |
var et = (pow(a, 2) - pow(b, 2)) / pow(b, 2) * pow(cos(lat), 2); | |
var n1 = pow(a, 2) / (sqrt(1 + et) * b); | |
var t = tan(lat); | |
var a1 = n1 * cos(lat); | |
var a2 = -(n1 * t * pow(cos(lat), 2)) / 2.0; | |
var a3 = -(n1 * pow(cos(lat), 3)) * (1 - pow(t, 2) + et) / 6; | |
var a4 = n1 * t * pow(cos(lat), 4) * (5 - pow(t, 2) + 9 * et + 4 * pow(et, 2)) / 24; | |
var a5 = n1 * pow(cos(lat), 5) * (5 - 18 * pow(t, 2) + pow(t, 4) + 14 * et - 58 * et * pow(t, 2)) / 120; | |
var a6 = n1 * t * pow(cos(lat), 6) * (61 - 58 * pow(t, 2) + pow(t, 4) + 270 * et - 330 * et * pow(t, 2)) / 720; | |
var north = -a2 * pow(lon, 2) + a4 * pow(lon, 4) + a6 * pow(lon, 6); | |
var east = a1 * lon - a3 * pow(lon, 3) + a5 * pow(lon, 5); | |
north = north + getGeodeticArch( | |
lat, | |
toRad(coordsystem.latitudeOfOrigin), | |
ellipsoidParams(ellipsoid) | |
); | |
//apply UTM "tricks" | |
return { | |
y: north * coordsystem.scaleFactor + coordsystem.falseNorthing, | |
x: east * coordsystem.scaleFactor + coordsystem.falseEasting | |
}; | |
}; | |
}()); |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment