WGS84 Projections

WGS84 Projections

For our WebGL Augmented Reality App, we needed to preprocess latitude and longitude values into localised coordinates. Lat/lng values are great, but they are not very useful if you need to represent them inside a 3D space. It is far better to use [x,y,z] coordinates which are stored as floats natively and can be used directly without projecting them into another coordinate space. Of course, you are better off using lat,lng values if you are looking at very large distances on a globe. For most other purposes (designing a game or 3d space where distances rendered are of smaller values) where a flat plane for a ground can be assumed, [x,y,z] coordinates are best.

Lets start by looking at what the formula looks like:

function toLocal(lat1, lng1, iLat, iLng, altitude) {

    //i for intial
    var R = 6371000; // radius of earth in metres
    var φ1 = degInRad(lat1);
    var φ2 = degInRad(iLat);
    var Δφ = degInRad(iLat - lat1);
    var Δλ = degInRad(iLng - lng1);
    var altitude = altitude || 0;

    var a = Math.sin(Δφ / 2) * Math.sin(Δφ / 2) +
        Math.cos(φ1) * Math.cos(φ2) *
        Math.sin(Δλ / 2) * Math.sin(Δλ / 2);
    var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

    var d = R * c;
    // see http://mathforum.org/library/drmath/view/55417.html
    
    var y = Math.sin(Δλ) * Math.cos(φ2);
    var x = Math.cos(φ1) * Math.sin(φ2) - Math.sin(φ1) * Math.cos(φ2) * Math.cos(Δλ);
    var θ = Math.atan2(y, x);

    var x = d * Math.cos(θ);
    var y = d * Math.sin(θ);
    return [x, altitude, y];
    // I use threejs THREE Vector for most of my projects as below
    //return new THREE.Vector(x,y,z)
    }

Since WGS84 coordinates have no way of representing height, this information needs to be added separately from another source. We used SRTM (Shuttle Radar Topography Mission) data open-sourced by NASA for our height maps (another post on how to use that later). I will not delve into the trigonometry of the formulas, but will explain the expected input and output.

The function takes in five values:

• lat1, lng1: the coordinates you want to convert
• iLat, ilng: the coordinates you want to convert with reference to
• altitude: the altitude of the lat1,lng1 point

The formula takes in these values and returns an [x,y,z] coordinate with:

• x is the x-axis distance in metres from point iLat, iLng point
• y is the height in metres above the ground
• z is the z-axis distance in metres from the iLat, iLng point

Note that our coordinate system was using height in the y-axis. If your system uses the altitude in z-axis, then you will need to change the return order accordingly.

If you don’t want to pass in iLat and iLng everytime, you can hardcode a reference value inside the function. It doesn’t matter if the reference point is a very large distance from the point your are projecting, as long as all the points to you want projected to local coordinates are clustered together. This means if you want to project the height map of Mount Everest and take iLat,iLng as [0,0], you will still get the right relative coordinates. The absilute distances will be very large, but relative to one another they will still be accurate.

I will take a look at using SRTM data in the next post. Also, a few helper functions for help in other GIS data conversions are below.

// returns the triangular x and y distance between two coordinates.
// absolute distance can be calculated through pythagoras: d = sqrt(x^2 + y^2)

function distanceBetween(lat1, lng1, lat2, lng2) {
    var x = distance(lat1, lng1, lat1, lng2);
    var y = distance(lat1, lng2, lat2, lng2);

    function distance(lat1, lng1, lat2, lng2) {
        var R = 6371000; // metres
        var φ1 = degInRad(lat1);
        var φ2 = degInRad(lat2);
        var Δφ = degInRad(lat2 - lat1);
        var Δλ = degInRad(lng2 - lng1);

        var a = Math.sin(Δφ / 2) * Math.sin(Δφ / 2) +
            Math.cos(φ1) * Math.cos(φ2) *
            Math.sin(Δλ / 2) * Math.sin(Δλ / 2);
        var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

        var d = R * c;
        return d;
    }

    return {
        v: x,
        h: y,
    };
};

// returns the midpoint wgs84 coordinate between two wgs84 coordinates
function midpoint(lat1, lng1, lat2, lng2) {
    var Δλ = degInRad(lng2 - lng1);

    //convert to radians
    lat1 = degInRad(lat1);
    lat2 = degInRad(lat2);
    lng1 = degInRad(lng1);

    var Bx = Math.cos(lat2) * Math.cos(Δλ);
    var By = Math.cos(lat2) * Math.sin(Δλ);
    var lat3 = Math.atan2(Math.sin(lat1) + Math.sin(lat2), Math.sqrt((Math.cos(lat1) + Bx) * (Math.cos(lat1) + Bx) + By * By));
    var lng3 = lng1 + Math.atan2(By, Math.cos(lat1) + Bx);

    return {
        lat: lat3 * 180 / Math.PI,
        lng: lng3 * 180 / Math.PI
    };

function degInRad(deg) {
    return deg * Math.PI / 180;
}