On 2005-07-31, Simon Brooke <

[email protected]> wrote:

> conversion to latitude/longitude. So I've got to recast this lot

><URL:http://www.gps.gov.uk/guideB.pdf> into Java. The really frustrating
I threw this together a while ago. Does it help?

Sue.

public class LatLongToOSGB

{

private static double scalechange = (20.4894 / 1000000.0);

private static double rz = (-0.8421);

private static double ry = (-0.247);

private static double rx = (-0.1502);

private static double r0c0 = (1 + scalechange);

private static double r0c1 = (-1 * rz * 0.000004848136811);

private static double r0c2 = (ry * 0.000004848136811);

private static double r1c0 = (rz * 0.000004848136811);

private static double r1c1 = r0c0;

private static double r1c2 = (-1 * rx * 0.000004848136811);

private static double r2c0 = (-1 * ry * 0.000004848136811);

private static double r2c1 = (rx * 0.000004848136811);

private static double r2c2 = r0c0;

private static double pi = 3.1415926535;

private static double NGa = 6377563.396;

private static double NGb = 6356256.910;

private static double WGS84a = 6378137.0;

private static double WGS84b = 6356752.3141;

private static double NGFo = 0.9996012717;

private static double NGphiOrigin = (49.0 * pi/180);

private static double NGlamdaOrigin = (-2.0 * pi/180);

private static double NGeOrigin = 400000.0;

private static double NGnOrigin = -100000.0;

private static double m_aSquared = (NGa * NGa);

private static double m_bSquared = (NGb * NGb);

private static double m_eSquared = ((m_aSquared - m_bSquared)/ m_aSquared);

private static double m_n = ((NGa - NGb)/ (NGa + NGb));

private static double m_nToTheTwo = (m_n * m_n);

private static double m_nToTheThree = (m_nToTheTwo * m_n);

private static double m_NGaNGFo = (NGa * NGFo);

private static double m_NGbNGFo = (NGb * NGFo);

private static double m_WGS84aSquared = (WGS84a * WGS84a);

private static double m_WGS84bSquared = (WGS84b * WGS84b);

private static double m_WGS84eSquared = ((m_WGS84aSquared - m_WGS84bSquared)/ m_WGS84aSquared);

private double northing;

private double easting;

public double phi;

public double lamda;

public double latitude;

public double longitude;

public LatLongToOSGB()

{

}

private void validateLatLong() {

// some validation

// if (latitude < 30.0)

// System.err.println("Too far South!");

// if (latitude > 70.0)

//System.err.println("Too far North!");

//if (longitude < -10.0)

//System.err.println("Too far West!");

//if (longitude > 3.0)

//System.err.println("Too far East!");

}

private void validateOSGB()

{

// some validation

if (northing < 0.0)

System.err.println("Too far South!");

if (northing > 1300000.0)

System.err.println("Too far North!");

if (easting < 0.0)

System.err.println("Too far West!");

if (easting > 700000.0)

System.err.println("Too far East!");

}

private void HelTrans(double ElipsoidHeight)

{

// Helmert Transformation

// lat long height

//This lat long to 3-d uses WGS84 m_eSquaredIN, NGaIN

double v = WGS84a/Math.sqrt(1 - m_WGS84eSquared * Math.pow(Math.sin(phi), 2) );

double inCoOrdX = (v + ElipsoidHeight)* Math.cos(phi) * Math.cos(lamda);

double inCoOrdY = (v + ElipsoidHeight)* Math.cos(phi) * Math.sin(lamda);

double inCoOrdZ = ((1 - m_WGS84eSquared)* v + ElipsoidHeight) * Math.sin(phi);

double prod1 = (1 + scalechange) * inCoOrdX +

(-1 * rz * 0.000004848136811) * inCoOrdY +

(ry * 0.000004848136811) * inCoOrdZ +

-446.448;

double prod2 = (rz * 0.000004848136811) * inCoOrdX +

(1 + scalechange) * inCoOrdY +

(-1 * rx * 0.000004848136811) * inCoOrdZ +

125.157;

double prod3 = (-1 * ry * 0.000004848136811) * inCoOrdX +

(rx * 0.000004848136811) * inCoOrdY +

(1 + scalechange) * inCoOrdZ +

-542.06;

// This 3-d to lat long uses OSGB m_eSquared, NGa

double p = Math.sqrt(prod1 * prod1 + prod2 * prod2);

double lat1 = Math.atan2(prod3, p * (1 - m_eSquared));

double v1 = NGa / Math.sqrt(1 - m_eSquared * Math.pow(Math.sin(lat1), 2));

double lat2 = Math.atan2(prod3 + m_eSquared * v1 * Math.sin(lat1), p);

double v2 = NGa / Math.sqrt(1 - m_eSquared* Math.pow(Math.sin(lat2), 2));

phi = Math.atan2(prod3 + m_eSquared * v2 * Math.sin(lat2), p);

lamda = Math.atan2(prod2, prod1);

ElipsoidHeight = p /Math.cos(phi) - v2;

}

void RevHelmert()

{

// Reverse Helmert Transformation

// lat long height

// This lat long to 3-d uses OSGB m_eSquared, NGa

double v = NGa/Math.sqrt(1.0 - m_eSquared * Math.pow(Math.sin(phi), 2) );

double inCoOrdX = (v )* Math.cos(phi) * Math.cos(lamda);

double inCoOrdY = (v )* Math.cos(phi) * Math.sin(lamda);

double inCoOrdZ = ((1.0 - m_eSquared)* v ) * Math.sin(phi);

// All parameter signs reversed.

double prod1 = (1.0 - scalechange) * inCoOrdX +

(rz * 0.000004848136811) * inCoOrdY +

(ry * -0.000004848136811) * inCoOrdZ +

446.448;

double prod2 = (rz * -0.000004848136811) * inCoOrdX +

(1.0 - scalechange) * inCoOrdY +

(rx * 0.000004848136811) * inCoOrdZ +

-125.157;

double prod3 = (ry * 0.000004848136811) * inCoOrdX +

(rx * -0.000004848136811) * inCoOrdY +

(1.0 - scalechange) * inCoOrdZ +

542.06;

// This 3-d to lat long uses WGS84 m_eSquaredIN, NGaIN

double p = Math.sqrt(prod1 * prod1 + prod2 * prod2);

double lat1 = Math.atan2(prod3, p * (1.0 - m_WGS84eSquared));

double v1 = WGS84a / Math.sqrt(1.0 - m_WGS84eSquared * Math.pow(Math.sin(lat1), 2));

double lat2 = Math.atan2(prod3 + m_WGS84eSquared * v1 * Math.sin(lat1), p);

double v2 = WGS84a /Math.sqrt(1.0 - m_WGS84eSquared * Math.pow(Math.sin(lat2), 2));

phi = Math.atan2(prod3 + m_WGS84eSquared * v2 * Math.sin(lat2), p);

v2 = WGS84a /Math.sqrt(1.0 - m_WGS84eSquared* Math.pow(Math.sin(phi), 2)); // we iterate once more.

phi = Math.atan2(prod3 + m_WGS84eSquared * v2 * Math.sin(phi), p);

lamda = Math.atan2(prod2, prod1);

}

private void doCalc(boolean bApplyOSGB1936Transform, double ElipsoidHeight)

{

if (bApplyOSGB1936Transform)

HelTrans(ElipsoidHeight);

double phisine = Math.sin(phi);

double phisineSquared = phisine * phisine;

double phicos = Math.cos(phi);

double phicosCubed = phicos * phicos * phicos;

double phicosToTheFifth = phicos * phicos * phicosCubed; // ****************************** FIXED

double phiTan = Math.tan(phi);

double phiTanSquared = phiTan * phiTan;

double phiTanToTheFourth = phiTanSquared * phiTanSquared;

double OneMinuseSquaredphisineSquared = 1.0 - m_eSquared * phisineSquared;

double phiMinusNGphiOrigin = phi - NGphiOrigin;

double phiPlusNGphiOrigin = phi + NGphiOrigin;

double v = m_NGaNGFo / Math.sqrt(OneMinuseSquaredphisineSquared);

double p = m_NGaNGFo * (1.0 - m_eSquared) * Math.pow((OneMinuseSquaredphisineSquared), -1.5);

double nnSquared = v/p - 1.0;

double frag1 = (1.0 + m_n + (5.0/4.0) * (m_nToTheTwo + m_nToTheThree)) * phiMinusNGphiOrigin;

double frag2 = (3.0 * (m_n + m_nToTheTwo) + (21.0/8.0) * m_nToTheThree) * Math.sin(phiMinusNGphiOrigin) * Math.cos(phiPlusNGphiOrigin);

double frag3 = ((15.0/8.0) * (m_nToTheTwo + m_nToTheThree)) * Math.sin(2.0 * (phiMinusNGphiOrigin)) * Math.cos(2.0 * (phiPlusNGphiOrigin));

double frag4 = ((35.0/24) * m_nToTheThree) * Math.sin(3.0 * (phiMinusNGphiOrigin)) * Math.cos(3.0 * (phiPlusNGphiOrigin));

double lamdaDiff = lamda - NGlamdaOrigin;

double lamdaDiffSquared = lamdaDiff * lamdaDiff;

double lamdaDiffToTheFourth = lamdaDiffSquared * lamdaDiffSquared;

{

double I = m_NGbNGFo * (frag1 - frag2 + frag3 - frag4) + NGnOrigin;

double II = (v/2.0) * phisine * phicos;

double III = (v/24.0) * phisine * phicosCubed * (5 - phiTanSquared + 9 * nnSquared);

double IIIA = (v/720.0) * phisine * phicosToTheFifth * (61 - 58 * phiTanSquared + phiTanToTheFourth);

northing = I + II * lamdaDiffSquared + III * lamdaDiffToTheFourth + IIIA * lamdaDiffToTheFourth * lamdaDiffSquared;

}

{

double IV = v * phicos;

double V = (v/6.0) * phicosCubed * (v/p - phiTanSquared);

double VI = (v/120.0) * phicosToTheFifth * (5 - 18 * phiTanSquared + phiTanToTheFourth + 14 * nnSquared - (58 * (phiTanSquared * nnSquared)) );

easting = NGeOrigin + IV * lamdaDiff + V * lamdaDiffSquared * lamdaDiff + VI * lamdaDiffToTheFourth * lamdaDiff;

}

}

void RevCalc()

{

// See page 41 of guide to OSGB coordinate system

// easting = 651409.903;

// northing = 313177.270;

double M = 0.0;

double phi_dash = ((northing - NGnOrigin - M)/ (NGa * NGFo)) + NGphiOrigin;

do

{

double term1 = (1 + m_n + (5.0 / 4.0) * m_nToTheTwo + (5.0 / 4.0) * m_nToTheThree) * (phi_dash - NGphiOrigin);

double term2 = (3.0 * m_n + 3.0 * m_nToTheTwo + (21.0 / 8.0) * m_nToTheThree) * Math.sin(phi_dash - NGphiOrigin) * Math.cos(phi_dash + NGphiOrigin);

double term3 = ((15.0 / 8.0) * m_nToTheTwo + (15.0 / 8.0) * m_nToTheThree) * Math.sin(2.0 *(phi_dash - NGphiOrigin))

* Math.cos(2.0 * (phi_dash + NGphiOrigin));

double term4 = ((35.0 / 24.0) * m_nToTheThree) * Math.sin(3.0 *(phi_dash - NGphiOrigin)) * Math.cos(3.0 * (phi_dash + NGphiOrigin));

M = NGb * NGFo * (term1 - term2 + term3 - term4);

phi_dash += ((northing - NGnOrigin - M)/ (NGa * NGFo));

}

while ((northing - NGnOrigin - M) >= 0.001);

double v = (NGa * NGFo) * Math.pow((1.0 - m_eSquared * Math.sin(phi_dash) * Math.sin(phi_dash)), -0.5);

double p = (NGa * NGFo) * (1.0 - m_eSquared) * Math.pow((1.0 - m_eSquared * Math.sin(phi_dash) * Math.sin(phi_dash)), -1.5);

double nsquared = v/p -1;

double tan_phi_dash = Math.tan(phi_dash);

double tan_phi_dash_squared = tan_phi_dash * tan_phi_dash;

double pv = p * v;

double VII = tan_phi_dash / (2.0 * pv);

double vsquared = v * v;

double VIII = (tan_phi_dash / (24.0 * pv * vsquared)) * (5.0 + 3.0 * tan_phi_dash_squared +

nsquared - 9.0 * tan_phi_dash_squared * nsquared);

double IX = (tan_phi_dash/ (720.0 * pv * vsquared * vsquared)) * (61.0 + 90.0 * tan_phi_dash_squared +

+ 45.0 * (tan_phi_dash_squared * tan_phi_dash_squared));

double secx = 1.0 / Math.cos(phi_dash);

double X = secx / v;

double XI = secx / (6.0 * v * vsquared) * ( v/p + 2 * tan_phi_dash_squared);

double vforth = vsquared * vsquared;

double XII = secx / (120.0 * v * vforth) * (5.0 + 28.0 * tan_phi_dash_squared + 24.0 * tan_phi_dash_squared * tan_phi_dash_squared);

double tan_phi_dash_sixth = tan_phi_dash_squared * tan_phi_dash_squared * tan_phi_dash_squared;

double XIIA = secx / (5040.0 * v * vforth * vsquared) * (61.0 + 662.0 * tan_phi_dash_squared +

1320.0 * tan_phi_dash_squared * tan_phi_dash_squared +

720.0 * tan_phi_dash_sixth);

double edeltasquared = (easting - NGeOrigin) * (easting - NGeOrigin);

phi = phi_dash - VII * edeltasquared +

VIII * edeltasquared * edeltasquared -

IX * edeltasquared * edeltasquared * edeltasquared;

lamda = NGlamdaOrigin + X * (easting - NGeOrigin) -

XI * (easting - NGeOrigin) * edeltasquared +

XII * (easting - NGeOrigin) * edeltasquared * edeltasquared -

XIIA * (easting - NGeOrigin) * edeltasquared * edeltasquared * edeltasquared;

RevHelmert();

latitude = phi * 180 /pi;

longitude = lamda * 180 /pi;

}

public String getOSGBString()

{

int Northing = (int)(northing + 0.5);

int Easting = (int)(easting + 0.5);

int osgbEasting = Easting / 100000;

int osgbNorthing = 14 - (Northing / 100000);

int iNorth = ((osgbNorthing * 5) /25) * 5 + 'H';

if (osgbEasting > 4)

iNorth++;

if (iNorth >= 'I')

iNorth++;

int iEast = (osgbEasting % 5) + 'A';

iEast += (osgbNorthing * 5) % 25;

if (iEast >= 'I')

iEast++;

String result_ = "";

result_ += (char)iNorth;

result_ += (char)iEast;

result_ += (Easting % 100000);

result_ += (Northing % 100000);

return result_;

}

public boolean setOSGBString(String sz1) {

String sz = "";

for (int xx = 0; xx < sz1.length(); xx++)

{

if (Character.isLetterOrDigit(sz1.charAt(xx)))

sz += sz1.charAt(xx);

}

int len = sz.length();

if (len == 0)

return false;

if (len % 2 == 1)

return false;

// pull apart the 2 initial letters - if present

if (Character.isLetter(sz.charAt(0)) != Character.isLetter(sz.charAt(1)))

return false;

String szNumericField;

easting = 0;

northing = 0;

if (Character.isLetter(sz.charAt(0)))

{

szNumericField = sz.substring(2);

if (!Convert(sz.charAt(0), sz.charAt(1))) {

return false;

}

} else {

szNumericField = sz;

}

szNumericField.trim();

// this field must have an even number of digits

len = szNumericField.length();

if (len % 2 == 1)

return false;

int factor = 10000;

for (int x = 0; x < len / 2; x++)

{

if (!Character.isDigit(szNumericField.charAt(x)))

return false;

if (!Character.isDigit(szNumericField.charAt(x + (len / 2))))

return false;

easting += factor * (szNumericField.charAt(x) & 0x0f);

northing += factor * (szNumericField.charAt(x + (len / 2)) & 0x0f);

factor = factor / 10;

}

RevCalc();

return true;

}

boolean Convert(char ch1, char ch2)

{

ch1 &= 0x1f;

if (ch1 > 9)

ch1--; // skip i

ch1--;

ch2 &= 0x1f;

if (ch2 > 9)

ch2--; // skip i

ch2--;

easting = (ch2 % 5) * 100000;

northing = (4 - (ch2 / 5)) * 100000;

easting += ((ch1 % 5) - 2) * 500000;

northing += (3 - (ch1 / 5)) * 500000;

if ((easting >= 700000) || (northing >= 1300000) || (easting < 0) || (northing < 0))

return false;

return true;

}

public void setLatLong(double latitude, double longitude) {

this.latitude = latitude;

this.longitude = longitude;

phi = latitude * pi / 180.0;

lamda = longitude * pi / 180.0;

doCalc(true, 0.0);

validateOSGB();

}

public double getLatitude() {

return latitude;

}

public double getLongitude() {

return longitude;

}

public String toString() {

String s = "OSGB " + getOSGBString() + " lat/long " + getLatitude() + " " + getLongitude();

return s;

}

public static void main(String[] argv) {

LatLongToOSGB latLong = new LatLongToOSGB();

if (argv.length == 1){

latLong.setOSGBString(argv[0]);

} else if (argv.length == 2) {

latLong.setLatLong(new Double(argv[0]).doubleValue(),

new Double(argv[1]).doubleValue());

} else {

System.out.println("Usage single arg OS string or 2 args for lat and long");

System.exit(1);

}

System.out.println(latLong);

}

}