Skip to content

Instantly share code, notes, and snippets.

@atlefren atlefren/toUtm33.js
Last active Aug 29, 2015

Embed
What would you like to do?
/**
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
You can’t perform that action at this time.