WGS-84 to OSGB36 Coords Conversion

OSGB.java

/**
 * Java based conversion of Chris Veness' scripts for Lat/Lon WG84 to OSGB36 coords
 *
 * http://www.movable-type.co.uk/scripts/latlong-convert-coords.html
 * (c) Chris Veness 2005-2010  Released under an LGPL license
 * http://www.fsf.org/licensing/licenses/lgpl.html
 * 
 */
public class OSGB {

    // WGS 84 ELLIPSOID
    private static double WGS84_A = 6378137;
    private static double WGS84_B = 6356752.314245;
    private static double WGS84_E2 = ((WGS84_A * WGS84_A) - (WGS84_B * WGS84_B)) / (WGS84_A * WGS84_A);


    // NatGrid scale factor on central meridian
    private static final double F0 = 0.9996012717;

    // Airy 1830 major & minor semi-axes - note .909 not .910
    // http://www.ordnancesurvey.co.uk/oswebsite/gps/docs/A_Guide_to_Coordinate_Systems_in_Great_Britain.pdf
    private static final double AIRY_A = 6377563.396;
    private static final double AIRY_B = 6356256.909;

    // NatGrid true origin
    private static final double LAT0 = Math.toRadians(49);  // Phi0
    private static final double LON0 = Math.toRadians(-2);  // Lamda0

    // northing & easting of true origin, metres
    private static final double N0 = -100000;
    private static final double E0 = 400000;

    // eccentricity squared
    private static final double E2 = ((AIRY_A * AIRY_A) - (AIRY_B * AIRY_B)) / (AIRY_A * AIRY_A);

    private static final double N = (AIRY_A - AIRY_B) / (AIRY_A + AIRY_B);
    private static final double N2 = N * N;
    private static final double N3 = N * N * N;


    // http://www.ordnancesurvey.co.uk/oswebsite/gps/docs/A_Guide_to_Coordinate_Systems_in_Great_Britain.pdf
    // section 6.6
    private static final double TX = -446.448;
    private static final double TY = 125.157;
    private static final double TZ = -542.060;
    private static final double RX = Math.toRadians(-0.1502 / 3600);
    private static final double RY = Math.toRadians(-0.2470 / 3600);
    private static final double RZ = Math.toRadians(-0.8421 / 3600);
    private static final double S = 20.4894 / 1e6 + 1;



    /**
     * Converts lat and lon in WGS84 datum to lat/lon in Airy
     *
     * @param lat   Latitude in radians
     * @param lon   Longitude in radians
     *
     * @return  Easting/Northing pair
     */
    public static double[] transform(double lat, double lon) {
        // -- 1: convert polar to cartesian coordinates (using ellipse 1)

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

        double nu = WGS84_A / Math.sqrt(1 - WGS84_E2 * sinPhi * sinPhi);

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

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

        double x2 = TX + x1 * S - y1 * RZ + z1 * RY;
        double y2 = TY + x1 * RZ + y1 * S - z1 * RX;
        double z2 = TZ - x1 * RY + y1 * RX + z1 * S;

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

        double precision = 4 / AIRY_A;

        double p = Math.sqrt(x2 * x2 + y2 * y2);
        double phi = Math.atan2(z2, p * (1 - E2)), phiP = 2 * Math.PI;
        while (Math.abs(phi - phiP) > precision) {
            nu = AIRY_A / Math.sqrt(1 - E2 * Math.sin(phi) * Math.sin(phi));
            phiP = phi;
            phi = Math.atan2(z2 + E2 * nu * Math.sin(phi), p);
        }
        double lambda = Math.atan2(y2, x2);

        // -- 4: now we're in OSGB, get the EN coords
        double cosLat = Math.cos(phi), sinLat = Math.sin(phi);
        nu = AIRY_A * F0 / Math.sqrt(1 - E2 * sinLat * sinLat);              // transverse radius of curvature
        double rho = AIRY_A * F0 * (1 - E2) / Math.pow(1 - E2 * sinLat * sinLat, 1.5);  // meridional radius of curvature
        double eta2 = nu / rho - 1;

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

        double cos3lat = cosLat * cosLat * cosLat;
        double cos5lat = cos3lat * cosLat * cosLat;
        double tan2lat = Math.tan(phi) * Math.tan(phi);
        double tan4lat = tan2lat * tan2lat;

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

        double dLon = lambda - LON0;
        double dLon2 = dLon * dLon, dLon3 = dLon2 * dLon, dLon4 = dLon3 * dLon, dLon5 = dLon4 * dLon, dLon6 = dLon5 * dLon;

        double N = I + II * dLon2 + III * dLon4 + IIIA * dLon6;
        double E = E0 + IV * dLon + V * dLon3 + VI * dLon5;

        return new double[] {E, N};

    }


    public static void main(String[] args) {
        double[] en = transform(Math.toRadians(50.687874), Math.toRadians(-1.940884));
        System.out.println(String.format("Easting = %f, Northing = %f", en[0], en[1]));
    }


}

Hi,

Thanks for your great website;
I am incorporating this code in simulation with MATSim using Java. But I am getting an error on Line 61 .... transform(double lat, double lon)
It is saying, syntax error on token "()" and "," . I changes "," to ";" and the error for this part has gone, but I don't know what to do with the error on "()".

I wonder if you have any suggestion.

Thank you,
Somayeh

Hey Somayeh,

Without seeing how you're integrating it, I can't be of much help. You should be able to take the class verbatim and put it into whatever package you need and then call the static method as I have done in the main method.

Alex