How to convert S-JTSK coordinates to WSG84

drndos March 12, 2013 1 Comment

Here in Slovakia and also in Czech republic our government uses special spatial geographic system called S-JTSK. This system uses x, y and z axis. However global standard is WSG84 (GPS position with latitude and longitude). Problem occurs while we want to convert this S-JTSK to WSG84. There are some programs which can do this but there are no good lightweight libraries which could be used and would be fast. So I analyzed one of these convertors in excel, extracted formulas and implemented it in Java.

First thing we need to do is to convert S-JTSK x, y to S-JTSK latitude and longitude.

[java collapse=”true”]
public double[] toSJTSKLatLong(double x, double y) {

double[] sjtskLatLong = new double[2];

double fi0;

double la0deg;

double fi0rad;

double la0;

double betakv;

double m;

double gamma;

double lambv;

double ro1;

double rGauss;

double c;

double betav;

double fi;

double la;

double a;

double b;

double e1;

double e2;

double finRad;

double vn;

double fin;

double n;

double k;

double laRad;

double fik;

double tga;

double fi1Rad;

double tga2;

double fi2Rad;

double tga3;

double fi3Rad;

double tga4;

double fi4Rad;

double tga5;

double firad;

double sjtskLat;

double laFerro;

double sjtskLong;

 

finRad = Math.PI * 49.5 / 180;

la0deg = 42 + 31.00 / 60 + 31.41725 / 3600;

ro1 = Math.sqrt(x * x + y * y);

rGauss = 6380065.5402;

betakv = Math.PI * 11.5 / 180;

fi0 = 59 + (42.00 / 60) + (42.69689 / 3600);

a = 6377397.155;

b = 6356078.96325;

 

 

e2 = Math.sqrt((a * a – b * b) / (b * b));

m = Math.cos(betakv);

vn = Math.sqrt(1 + e2 * e2 * Math.cos(finRad) * Math.cos(finRad));

gamma = Math.atan2(y, x);

c = rGauss * Math.sin(betakv) / m / Math.pow(Math.tan(betakv / 2), m);

 

fi0rad = fi0 * Math.PI / 180;

fin = Math.atan2(Math.tan(finRad), vn);

lambv = gamma / m;

betav = 2 * Math.atan(Math.pow(ro1 / c, 1 / m));

fi = Math.asin(Math.cos(betav) * Math.sin(fi0rad) – Math.sin(betav) * Math.cos(fi0rad) * Math.cos(lambv));

 

n = Math.sin(finRad) / Math.sin(fin);

la = Math.asin(Math.sin(betav) * Math.sin(lambv) / Math.cos(fi));

la0 = la0deg * Math.PI / 180;

 

laRad = (la0 – la) / n;

laFerro = laRad * 180 / Math.PI;

sjtskLong = laFerro – fgshift;

e1 = Math.sqrt((a * a – b * b) / (a * a));

 

k = Math.tan(Math.PI / 4 + fin / 2) / (Math.pow(Math.tan(Math.PI / 4 + finRad / 2), n) * Math.pow((1 – e1 * Math.sin(finRad)) / (1 + e1 * Math.sin(finRad)), n * e1 / 2));

fik = 2 * Math.atan(Math.pow(1 / k * Math.tan(Math.PI / 4 + fi / 2), 1 / n)) – Math.PI / 2;

tga = Math.pow((1 – e1 * Math.sin(fik)) / (1 + e1 * Math.sin(fik)), n * e1 / 2);

fi1Rad = 2 * Math.atan(Math.pow(1 / k / tga * Math.tan(Math.PI / 4 + fi / 2), 1 / n)) – Math.PI / 2;

tga2 = Math.pow((1 – e1 * Math.sin(fi1Rad)) / (1 + e1 * Math.sin(fi1Rad)), n * e1 / 2);

fi2Rad = 2 * Math.atan(Math.pow(1 / k / tga2 * Math.tan(Math.PI / 4 + fi / 2), 1 / n)) – Math.PI / 2;

tga3 = Math.pow((1 – e1 * Math.sin(fi2Rad)) / (1 + e1 * Math.sin(fi2Rad)), n * e1 / 2);

fi3Rad = 2 * Math.atan(Math.pow(1 / k / tga3 * Math.tan(Math.PI / 4 + fi / 2), 1 / n)) – Math.PI / 2;

tga4 = Math.pow((1 – e1 * Math.sin(fi3Rad)) / (1 + e1 * Math.sin(fi3Rad)), n * e1 / 2);

fi4Rad = 2 * Math.atan(Math.pow(1 / k / tga4 * Math.tan(Math.PI / 4 + fi / 2), 1 / n)) – Math.PI / 2;

tga5 = Math.pow((1 – e1 * Math.sin(fi4Rad)) / (1 + e1 * Math.sin(fi4Rad)), n * e1 / 2);

firad = 2 * Math.atan(Math.pow(1 / k / tga5 * Math.tan(Math.PI / 4 + fi / 2), 1 / n)) – Math.PI / 2;

sjtskLat = firad * 180 / Math.PI;

 

sjtskLatLong[0] = sjtskLat;

sjtskLatLong[1] = sjtskLong;

return sjtskLatLong;

}[/java]

When we have S-JTSK long and lat we can convert these into WSG84 lat long like this:

[java collapse=”true”]
public double[] toWSG84LatLong(double SJTSKlat, double SJTSKlong) {

double[] wsgLatLong = new double[2];

double fi1Rad;

double la1Rad;

double a1;

double f1;

double a2;

double f2;

double e1;

double M;

double N;

double dfisec;

double dlasec;

double wsgLat;

double wsgLong;

fi1Rad = SJTSKlat * Math.PI / 180;

la1Rad = SJTSKlong * Math.PI / 180;

a1 = 6377397.155;

f1 = 0.00334277318217481;

a2 = 6378137.00;

f2 = 0.00335281066474748;

e1 = Math.sqrt(2 * f1 – f1 * f1);

M = a1 * (1 – e1 * e1) / Math.pow((1 – e1 * e1 * Math.sin(fi1Rad) * Math.sin(fi1Rad)), 1.5);

N = a1 / Math.sqrt(1 – e1 * e1 * Math.sin(fi1Rad) * Math.sin(fi1Rad));

dlasec = (-dXshift * Math.sin(la1Rad) + dYshift * Math.cos(la1Rad)) / (N * Math.cos(fi1Rad) * Math.sin(Math.PI / 180 / 3600));

dfisec = (-dXshift * Math.sin(fi1Rad) * Math.cos(la1Rad) – dYshift * Math.sin(fi1Rad) * Math.sin(la1Rad) + dZshift * Math.cos(fi1Rad) + (a1 * (f2 – f1) + f1 * (a2 – a1)) * Math.sin(2 * fi1Rad)) / (M * Math.sin(Math.PI / 180 / 3600));

wsgLong = SJTSKlong + dlasec / 3600;

wsgLat = SJTSKlat + dfisec / 3600;

wsgLatLong[0] = wsgLat;

wsgLatLong[1] = wsgLong;

return wsgLatLong;

}[/java]

We have to define constants used in these methods

[java collapse=”true”]
private final double fgshift = 17.66666666666666;

private final double dXshift = 589;

private final double dYshift = 76;

private final double dZshift = 480;
[/java]

And we are done. Complete source code is available here SJTSKtoWSG84.java