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);
}
}