MongoDB scripts for WGS84 -> OSGB conversions on the server side

monos.js

/*
 --------------------------------------------------------------------------------------------------------------------
 EXTENSIONS
 --------------------------------------------------------------------------------------------------------------------
 */

/** Converts numeric degrees to radians */
if (typeof Number.prototype.toRad == 'undefined') {
    Number.prototype.toRad = function () {
        return this * Math.PI / 180;
    };
}

/** Converts radians to numeric (signed) degrees */
if (typeof Number.prototype.toDeg == 'undefined') {
    Number.prototype.toDeg = function () {
        return this * 180 / Math.PI;
    };
}

/*
 --------------------------------------------------------------------------------------------------------------------
 GLOBAL NAMESPACE
 --------------------------------------------------------------------------------------------------------------------
 */
var os = (function(){
    var self = {};

    // -------------------------------------------------------------------
    // PRIVATES
    // -------------------------------------------------------------------
    var ellipse = {
        WGS84:      { a:6378137,     b:6356752.3142,    f:1 / 298.257223563  },
        Airy1830:   { a:6377563.396, b:6356256.909,     f:1 / 299.3249646    }
    };

    var datumTransform = {
        toOSGB36:{  tx:-446.448,    ty:125.157,     tz:-542.060,    // m
                    rx:-0.1502,     ry:-0.2470,     rz:-0.8421,     // sec
                    s:20.4894                                       // ppm
        },
        toED50:{    tx:89.5,        ty:93.8,        tz:123.1,       // m
                    rx:0.0,         ry:0.0,         rz:0.156,       // sec
                    s:-1.2                                          // ppm
        },
        toIrl1975:{ tx:-482.530,    ty:130.596,     tz:-564.557,    // m
                    rx:-1.042,      ry:-0.214,      rz:-0.631,      // sec
                    s:-8.150                                        // ppm
        }
    };


    // Airy 1830 major & minor semi-axes
    var A = ellipse.Airy1830.a, B = ellipse.Airy1830.b;

    // NatGrid scale factor on central meridian
    var F0 = 0.9996012717;

    // nat grid true origin
    var LAT0 = (49).toRad(), LON0 = (-2).toRad();

    // northing & easting of true origin, metres
    var N0 = -100000, E0 = 400000;

    // eccentricity squared
    var E2 = 1 - (B * B) / (A * A);


    var N = (A - B) / (A + B);
    var N2 = N * N;
    var N3 = N2 * N;

    // -------------------------------------------------------------------
    // PUBLICS
    // -------------------------------------------------------------------

    /**
     * Converts a lat/lon in WGS84 to lat/lon in OSGB36.
     *
     * @param lat   Latitude in decimal degrees
     * @param lon   Longitude in decimal degrees
     */
    self.toOSGB36 = function (lat, lon) {
        // -- 1: convert polar to cartesian coordinates (using ellipse 1)

        // source ellipsoid
        var e1 = ellipse.WGS84;
        var e2 = ellipse.Airy1830;
        var t = datumTransform.toOSGB36;

        lat = lat.toRad();
        lon = lon.toRad();

        var a = e1.a, b = e1.b;

        var sinPhi = Math.sin(lat);
        var cosPhi = Math.cos(lat);
        var sinLambda = Math.sin(lon);
        var cosLambda = Math.cos(lon);
        var H = 24.7;  // for the moment

        var eSq = (a * a - b * b) / (a * a);
        var nu = a / Math.sqrt(1 - eSq * sinPhi * sinPhi);

        var x1 = (nu + H) * cosPhi * cosLambda;
        var y1 = (nu + H) * cosPhi * sinLambda;
        var z1 = ((1 - eSq) * nu + H) * sinPhi;


        // -- 2: apply helmert transform using appropriate params

        var tx = t.tx, ty = t.ty, tz = t.tz;
        var rx = (t.rx / 3600).toRad();  // normalise seconds to radians
        var ry = (t.ry / 3600).toRad();
        var rz = (t.rz / 3600).toRad();
        var s1 = t.s / 1e6 + 1;          // normalise ppm to (s+1)

        // apply transform
        var x2 = tx + x1 * s1 - y1 * rz + z1 * ry;
        var y2 = ty + x1 * rz + y1 * s1 - z1 * rx;
        var z2 = tz - x1 * ry + y1 * rx + z1 * s1;


        // -- 3: convert cartesian to polar coordinates (using ellipse 2)

        a = e2.a;
        b = e2.b;
        var precision = 4 / a;  // results accurate to around 4 metres

        eSq = (a * a - b * b) / (a * a);
        var p = Math.sqrt(x2 * x2 + y2 * y2);
        var phi = Math.atan2(z2, p * (1 - eSq)), phiP = 2 * Math.PI;
        while (Math.abs(phi - phiP) > precision) {
            nu = a / Math.sqrt(1 - eSq * Math.sin(phi) * Math.sin(phi));
            phiP = phi;
            phi = Math.atan2(z2 + eSq * nu * Math.sin(phi), p);
        }
        var lambda = Math.atan2(y2, x2);
        H = p / Math.cos(phi) - nu;

        return {lat:phi.toDeg(), lon:lambda.toDeg(), height:H};
    };

    /**
     * Converts a lat/lon (in OSGB36) to E/N
     *
     * @param lat   Latitude in decimal degress
     * @param lon   Longitude in decimal degrees
     */
    self.toEN = function (lat, lon) {
        lat = lat.toRad();
        lon = lon.toRad();

        var cosLat = Math.cos(lat), sinLat = Math.sin(lat);
        var nu = A * F0 / Math.sqrt(1 - E2 * sinLat * sinLat);              // transverse radius of curvature
        var rho = A * F0 * (1 - E2) / Math.pow(1 - E2 * sinLat * sinLat, 1.5);  // meridional radius of curvature
        var eta2 = nu / rho - 1;

        var Ma = (1 + N + (5 / 4) * N2 + (5 / 4) * N3) * (lat - LAT0);
        var Mb = (3 * N + 3 * N * N + (21 / 8) * N3) * Math.sin(lat - LAT0) * Math.cos(lat + LAT0);
        var Mc = ((15 / 8) * N2 + (15 / 8) * N3) * Math.sin(2 * (lat - LAT0)) * Math.cos(2 * (lat + LAT0));
        var Md = (35 / 24) * N3 * Math.sin(3 * (lat - LAT0)) * Math.cos(3 * (lat + LAT0));
        var M = B * F0 * (Ma - Mb + Mc - Md);              // meridional arc

        var cos3lat = cosLat * cosLat * cosLat;
        var cos5lat = cos3lat * cosLat * cosLat;
        var tan2lat = Math.tan(lat) * Math.tan(lat);
        var tan4lat = tan2lat * tan2lat;

        var I = M + N0;
        var II = (nu / 2) * sinLat * cosLat;
        var III = (nu / 24) * sinLat * cos3lat * (5 - tan2lat + 9 * eta2);
        var IIIA = (nu / 720) * sinLat * cos5lat * (61 - 58 * tan2lat + tan4lat);
        var IV = nu * cosLat;
        var V = (nu / 6) * cos3lat * (nu / rho - tan2lat);
        var VI = (nu / 120) * cos5lat * (5 - 18 * tan2lat + tan4lat + 14 * eta2 - 58 * tan2lat * eta2);

        var dLon = lon - LON0;
        var dLon2 = dLon * dLon, dLon3 = dLon2 * dLon, dLon4 = dLon3 * dLon, dLon5 = dLon4 * dLon, dLon6 = dLon5 * dLon;

        var northing = I + II * dLon2 + III * dLon4 + IIIA * dLon6;
        var easting = E0 + IV * dLon + V * dLon3 + VI * dLon5;

        return {E:easting, N:northing};
    };

    self.toLatLon = function(easting, northing) {

        var lat = LAT0, M = 0;
        do {
            lat = (northing - N0 - M) / (A * F0) + lat;

            var Ma = (1 + N + (5 / 4) * N2 + (5 / 4) * N3) * (lat - LAT0);
            var Mb = (3 * N + 3 * N * N + (21 / 8) * N3) * Math.sin(lat - LAT0) * Math.cos(lat + LAT0);
            var Mc = ((15 / 8) * N2 + (15 / 8) * N3) * Math.sin(2 * (lat - LAT0)) * Math.cos(2 * (lat + LAT0));
            var Md = (35 / 24) * N3 * Math.sin(3 * (lat - LAT0)) * Math.cos(3 * (lat + LAT0));
            M = B * F0 * (Ma - Mb + Mc - Md);                // meridional arc

        } while (northing - N0 - M >= 0.00001);  // ie until < 0.01mm

        var cosLat = Math.cos(lat), sinLat = Math.sin(lat);
        var nu = A * F0 / Math.sqrt(1 - E2 * sinLat * sinLat);              // transverse radius of curvature
        var rho = A * F0 * (1 - E2) / Math.pow(1 - E2 * sinLat * sinLat, 1.5);  // meridional radius of curvature
        var eta2 = nu / rho - 1;

        var tanLat = Math.tan(lat);
        var tan2lat = tanLat * tanLat, tan4lat = tan2lat * tan2lat, tan6lat = tan4lat * tan2lat;
        var secLat = 1 / cosLat;
        var nu3 = nu * nu * nu, nu5 = nu3 * nu * nu, nu7 = nu5 * nu * nu;
        var VII = tanLat / (2 * rho * nu);
        var VIII = tanLat / (24 * rho * nu3) * (5 + 3 * tan2lat + eta2 - 9 * tan2lat * eta2);
        var IX = tanLat / (720 * rho * nu5) * (61 + 90 * tan2lat + 45 * tan4lat);
        var X = secLat / nu;
        var XI = secLat / (6 * nu3) * (nu / rho + 2 * tan2lat);
        var XII = secLat / (120 * nu5) * (5 + 28 * tan2lat + 24 * tan4lat);
        var XIIA = secLat / (5040 * nu7) * (61 + 662 * tan2lat + 1320 * tan4lat + 720 * tan6lat);

        var dE = (easting - E0), dE2 = dE * dE, dE3 = dE2 * dE, dE4 = dE2 * dE2, dE5 = dE3 * dE2, dE6 = dE4 * dE2, dE7 = dE5 * dE2;
        lat = lat - VII * dE2 + VIII * dE4 - IX * dE6;
        var lon = LON0 + X * dE - XI * dE3 + XII * dE5 - XIIA * dE7;

        return {lat:lat.toDeg(), lon:lon.toDeg()};
    };

    return self;

})();

db.system.js.save({"_id":"os", value:os});

test.tolatlon.js

db.pos.remove();

db.pos.insert({pos: {lat: 50.687767, lon: -1.941008}, name: "Sandbanks"});
db.pos.insert({pos: {lat: 50.937927, lon: -1.470626}, name: "Ordnance Survey"});
db.pos.insert({pos: {lat: 51.50338, lon: -0.11969}, name: "London Eye"});
db.pos.insert({pos: {lat: 52.62177, lon: -1.12444}, name: "Leicester University"});

db.pos.find().forEach(function(elm) {
    var ll = os.toOSGB36(elm.pos.lat, elm.pos.lon);
    var en = os.toEN(ll.lat, ll.lon);
    print(elm.name + ": (" + elm.pos.lon + "," + elm.pos.lat + ") => (" + en.E + "," + en.N + ")");
});

Quick test script for converting Lat/Lon in WGS84 into OSGB36 E/N using MongoDB. The code is modified from the excellent http://www.movable-type.co.uk/ javascript resource.

To run:

mongo monos.js test.tolatlon.js

You should see the following:

MongoDB shell version: 2.2.0
connecting to: test
loading file: monos.js
loading file: test.tolatlon.js
Sandbanks: (-1.941008,50.687767) => (404263.5874244009,87569.9587131209)
Ordnance Survey: (-1.470626,50.937927) => (437292.88340633264,115520.25779962823)
London Eye: (-0.11969,51.50338) => (530603.0252780108,179947.48600626152)
Leicester University: (-1.12444,52.62177) => (459369.5806165703,303026.8833740335)

To continue playing, add the --shell option when running the commands...