Převod WGS 84 do S-JTSK
Co se týče zeměpisných souřadnic, běžný smrtelník ví, že existuje nějaká GPSka. Programátor pracující na českých projektech, které vyžadují práci se zeměpisnými souřadnicemi se určitě setkal se dvěma standardy a to WGS 84 a S-JTSK. Proto některé služby jako například Webdispečink poskytují data v obou standardech. Nicméně někdy potřebujete mezi sebou souřadnice převést.
Je za tím spoustu matematiky různých projekcí atd., čemuž nehodlám porozumět. Převod z S-JTSK do WGS 84 není sranda, ale existuje na to knihovna PROJ.4 (zvládá více světových systémů). Návod pro český systém naleznete wiki ČVUT.
Převod z WGS 84 do S-JTSK je výpočetně jednodušší. Reverse engineeringem jsem z javascriptu na stránce http://pecina.cz/krovak.html vytvořil následujíc java kód.
public class JTSK { | |
private double x; | |
private double y; | |
public JTSK(double x, double y) { | |
this.x = x; | |
this.y = y; | |
} | |
public double getX() { | |
return x; | |
} | |
public void setX(double x) { | |
this.x = x; | |
} | |
public double getY() { | |
return y; | |
} | |
public void setY(double y) { | |
this.y = y; | |
} | |
/** | |
* převádí WGS-84 do S-JTSK | |
* | |
* @param WGS84 | |
* @return JTSK | |
*/ | |
public static JTSK convert(Wgs84 wgs84) { | |
final double altitude = wgs84.getAltitude(); | |
double d2r = Math.PI / 180; | |
double a = 6378137.0; | |
double f1 = 298.257223563; | |
double dx = -570.69; | |
double dy = -85.69; | |
double dz = -462.84; | |
double wx = 4.99821 / 3600 * Math.PI / 180; | |
double wy = 1.58676 / 3600 * Math.PI / 180; | |
double wz = 5.2611 / 3600 * Math.PI / 180; | |
double m = -3.543e-6; | |
double latitude = wgs84.getLatitude(); | |
double longtitude = wgs84.getLongitude(); | |
double B = latitude * d2r; | |
double L = longtitude * d2r; | |
double H = altitude; | |
double e2 = 1 - Math.pow(1 - 1 / f1, 2); | |
double rho = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(B), 2)); | |
double x1 = (rho + H) * Math.cos(B) * Math.cos(L); | |
double y1 = (rho + H) * Math.cos(B) * Math.sin(L); | |
double z1 = ((1 - e2) * rho + H) * Math.sin(B); | |
double x2 = dx + (1 + m) * (x1 + wz * y1 - wy * z1); | |
double y2 = dy + (1 + m) * (-wz * x1 + y1 + wx * z1); | |
double z2 = dz + (1 + m) * (wy * x1 - wx * y1 + z1); | |
a = 6377397.15508; | |
f1 = 299.152812853; | |
double ab = f1 / (f1 - 1); | |
double p = Math.sqrt(Math.pow(x2, 2) + Math.pow(y2, 2)); | |
e2 = 1 - Math.pow(1 - 1 / f1, 2); | |
double th = Math.atan(z2 * ab / p); | |
double st = Math.sin(th); | |
double ct = Math.cos(th); | |
double t = (z2 + e2 * ab * a * (st * st * st)) / (p - e2 * a * (ct * ct * ct)); | |
B = Math.atan(t); | |
H = Math.sqrt(1 + t * t) * (p - a / Math.sqrt(1 + (1 - e2) * t * t)); | |
L = 2 * Math.atan(y2 / (p + x2)); | |
a = 6377397.15508; | |
double e = 0.081696831215303; | |
double n = 0.97992470462083; | |
double rho0 = 12310230.12797036; | |
double sinUQ = 0.863499969506341; | |
double cosUQ = 0.504348889819882; | |
double sinVQ = 0.420215144586493; | |
double cosVQ = 0.907424504992097; | |
double alpha = 1.000597498371542; | |
double k2 = 1.00685001861538; | |
double sinB = Math.sin(B); | |
t = (1 - e * sinB) / (1 + e * sinB); | |
t = Math.pow(1 + sinB, 2) / (1 - Math.pow(sinB, 2)) * Math.exp(e * Math.log(t)); | |
t = k2 * Math.exp(alpha * Math.log(t)); | |
double sinU = (t - 1) / (t + 1); | |
double cosU = Math.sqrt(1 - sinU * sinU); | |
double V = alpha * L; | |
double sinV = Math.sin(V); | |
double cosV = Math.cos(V); | |
double cosDV = cosVQ * cosV + sinVQ * sinV; | |
double sinDV = sinVQ * cosV - cosVQ * sinV; | |
double sinS = sinUQ * sinU + cosUQ * cosU * cosDV; | |
double cosS = Math.sqrt(1 - sinS * sinS); | |
double sinD = sinDV * cosU / cosS; | |
double cosD = Math.sqrt(1 - sinD * sinD); | |
double eps = n * Math.atan(sinD / cosD); | |
rho = rho0 * Math.exp(-n * Math.log((1 + sinS) / cosS)); | |
double CX = rho * Math.sin(eps); | |
double CY = rho * Math.cos(eps); | |
return new JTSK(-CX, -CY); | |
} | |
} |
public class Wgs84 { | |
private double latitude; | |
private double longitude; | |
/** | |
* nadmořská výška v metrech | |
*/ | |
private double altitude = 200; | |
public Wgs84(double latitude, double longitude) { | |
this.latitude = latitude; | |
this.longitude = longitude; | |
} | |
public Wgs84(double latitude, double longitude, double altitude) { | |
this(latitude, longitude); | |
this.altitude = altitude; | |
} | |
public double getLatitude() { | |
return latitude; | |
} | |
public void setLatitude(double latitude) { | |
this.latitude = latitude; | |
} | |
public double getLongitude() { | |
return longitude; | |
} | |
public void setLongitude(double longitude) { | |
this.longitude = longitude; | |
} | |
public double getAltitude() { | |
return altitude; | |
} | |
public void setAltitude(double altitude) { | |
this.altitude = altitude; | |
} | |
} |