Index: org/openstreetmap/josm/data/projection/SwissGrid.java =================================================================== --- org/openstreetmap/josm/data/projection/SwissGrid.java (révision 3429) +++ org/openstreetmap/josm/data/projection/SwissGrid.java (copie de travail) @@ -3,7 +3,6 @@ package org.openstreetmap.josm.data.projection; import static org.openstreetmap.josm.tools.I18n.tr; - import org.openstreetmap.josm.data.Bounds; import org.openstreetmap.josm.data.coor.EastNorth; import org.openstreetmap.josm.data.coor.LatLon; @@ -14,77 +13,143 @@ * Calculations are based on formula from * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf * + * August 2010 update to this formula (rigorous formulas) + * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf */ public class SwissGrid implements Projection { + + private static final double dX = 674.374; + private static final double dY = 15.056; + private static final double dZ = 405.346; + + private static final double wgs84_a = 6378137.000; + private static final double wgs84_b = 6356752.314245; +// private static final double wgs84_f = (wgs84_a-wgs84_b)/wgs84_a; + private static final double wgs84_E2 = (Math.pow(wgs84_a, 2)-Math.pow(wgs84_b, 2))/Math.pow(wgs84_a, 2); +// private static final double wgs84_E = Math.sqrt(wgs84_E2); + + + private static final double a = 6377397.155; + private static final double b = 6356078.962822; +// private static final double f = (a-b)/a; + private static final double E2 = (Math.pow(a, 2)-Math.pow(b, 2))/Math.pow(a, 2); + private static final double E = Math.sqrt(E2); + private static final double phi0 = Math.toRadians(46.0 + 57.0 / 60 + 8.66 / 3600); + private static final double lambda0 = Math.toRadians(7.0 + 26.0 / 60 + 22.50 / 3600); + private static final double R = a * Math.sqrt(1 - E2) / (1 - (E2 * Math.pow(Math.sin(phi0), 2))); + private static final double alpha = Math.sqrt(1 + (E2 / (1 - E2) * Math.pow(Math.cos(phi0), 4))); + private static final double b0 = Math.asin(Math.sin(phi0) / alpha); + private static final double K = Math.log(Math.tan(Math.PI / 4 + b0 / 2)) - alpha + * Math.log(Math.tan(Math.PI / 4 + phi0 / 2)) + alpha * E / 2 + * Math.log((1 + E * Math.sin(phi0)) / (1 - E * Math.sin(phi0))); + + private static final double xTrans = 200000; + private static final double yTrans = 600000; + + private LatLon correctETRS89toCH1903(LatLon coord) { + double phi = Math.toRadians(coord.lat()); + double lambda = Math.toRadians(coord.lon()); + + double Rn = wgs84_a/Math.sqrt(1-wgs84_E2*Math.pow(Math.sin(phi), 2)); + double X = Rn*Math.cos(phi)*Math.cos(lambda); + double Y = Rn*Math.cos(phi)*Math.sin(lambda); + double Z = Rn*(1-wgs84_E2)*Math.sin(phi); + X -= dX; + Y -= dY; + Z -= dZ; + lambda = Math.atan(Y/X); + double sX2Ys = Math.sqrt(X*X+Y*Y); + phi = Math.atan((Z/sX2Ys)); + double h = 0; + + for (int i = 0 ; i < 6 ; i++) { + Rn = a/Math.sqrt(1-E2*Math.pow(Math.sin(phi), 2)); + h = sX2Ys/Math.cos(phi)-Rn; + phi = Math.atan((Z/sX2Ys)/(1-Rn*E2/(Rn+h))); + } + + return new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda)); + } + + private LatLon correctCH1903toETRS89(LatLon coord) { + double phi = Math.toRadians(coord.lat()); + double lambda = Math.toRadians(coord.lon()); + + double Rn = a/Math.sqrt(1-E2*Math.pow(Math.sin(phi), 2)); + double X = Rn*Math.cos(phi)*Math.cos(lambda); + double Y = Rn*Math.cos(phi)*Math.sin(lambda); + double Z = Rn*(1-E2)*Math.sin(phi); + X += dX; + Y += dY; + Z += dZ; + lambda = Math.atan(Y/X); + double sX2Ys = Math.sqrt(X*X+Y*Y); + phi = Math.atan((Z/sX2Ys)); + double h = 0; + + for (int i = 0 ; i < 5 ; i++) { + Rn = wgs84_a/Math.sqrt(1-wgs84_E2*Math.pow(Math.sin(phi), 2)); + h = sX2Ys/Math.cos(phi)-Rn; + phi = Math.atan((Z/sX2Ys)/(1-Rn*wgs84_E2/(Rn+h))); + } + + return new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda)); + } + /** * @param wgs WGS84 lat/lon (ellipsoid GRS80) (in degree) * @return eastnorth projection in Swiss National Grid (ellipsoid Bessel) */ public EastNorth latlon2eastNorth(LatLon wgs) { - double phi = 3600d * wgs.lat(); - double lambda = 3600d * wgs.lon(); + LatLon coord = correctETRS89toCH1903(wgs); + double phi = Math.toRadians(coord.lat()); + double lambda = Math.toRadians(coord.lon()); - double phiprime = (phi - 169028.66d) / 10000d; - double lambdaprime = (lambda - 26782.5d) / 10000d; + double S = alpha * Math.log(Math.tan(Math.PI / 4 + phi / 2)) - alpha * E / 2 + * Math.log((1 + E * Math.sin(phi)) / (1 - E * Math.sin(phi))) + K; + double b = 2 * (Math.atan(Math.exp(S)) - Math.PI / 4); + double l = alpha * (lambda - lambda0); - // precompute squares for lambdaprime and phiprime - // - double lambdaprime_2 = Math.pow(lambdaprime,2); - double phiprime_2 = Math.pow(phiprime,2); + double lb = Math.atan(Math.sin(l) / (Math.sin(b0) * Math.tan(b) + Math.cos(b0) * Math.cos(l))); + double bb = Math.asin(Math.cos(b0) * Math.sin(b) - Math.sin(b0) * Math.cos(b) * Math.cos(l)); - double north = - 200147.07d - + 308807.95d * phiprime - + 3745.25d * lambdaprime_2 - + 76.63d * phiprime_2 - - 194.56d * lambdaprime_2 * phiprime - + 119.79d * Math.pow(phiprime,3); + double y = R * lb; + double x = R / 2 * Math.log((1 + Math.sin(bb)) / (1 - Math.sin(bb))); - double east = - 600072.37d - + 211455.93d * lambdaprime - - 10938.51d * lambdaprime * phiprime - - 0.36d * lambdaprime * phiprime_2 - - 44.54d * Math.pow(lambdaprime,3); - - return new EastNorth(east, north); + return new EastNorth(y + yTrans, x + xTrans); } /** * @param xy SwissGrid east/north (in meters) * @return LatLon WGS84 (in degree) */ - public LatLon eastNorth2latlon(EastNorth xy) { - double yp = (xy.east() - 600000d) / 1000000d; - double xp = (xy.north() - 200000d) / 1000000d; + double x = xy.north() - xTrans; + double y = xy.east() - yTrans; - // precompute squares of xp and yp - // - double xp_2 = Math.pow(xp,2); - double yp_2 = Math.pow(yp,2); + double lb = y / R; + double bb = 2 * (Math.atan(Math.exp(x / R)) - Math.PI / 4); - // lambda = latitude, phi = longitude - double lmbp = 2.6779094d - + 4.728982d * yp - + 0.791484d * yp * xp - + 0.1306d * yp * xp_2 - - 0.0436d * Math.pow(yp,3); + double b = Math.asin(Math.cos(b0) * Math.sin(bb) + Math.sin(b0) * Math.cos(bb) * Math.cos(lb)); + double l = Math.atan(Math.sin(lb) / (Math.cos(b0) * Math.cos(lb) - Math.sin(b0) * Math.tan(bb))); - double phip = 16.9023892d - + 3.238272d * xp - - 0.270978d * yp_2 - - 0.002528d * xp_2 - - 0.0447d * yp_2 * xp - - 0.0140d * Math.pow(xp,3); + double lambda = lambda0 + l / alpha; + double phi = b; + double S = 0; - double lmb = lmbp * 100.0d / 36.0d; - double phi = phip * 100.0d / 36.0d; + // 5 iteration to fins S and phi + for (int i = 0; i < 6; i++) { + S = 1 / alpha * (Math.log(Math.tan(Math.PI / 4 + b / 2)) - K) + E + * Math.log(Math.tan(Math.PI / 4 + Math.asin(E * Math.sin(phi)) / 2)); + phi = 2 * Math.atan(Math.exp(S)) - Math.PI / 2; + } - return new LatLon(phi,lmb); + LatLon coord = new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda)); + return correctCH1903toETRS89(coord); } - @Override public String toString() { + @Override + public String toString() { return tr("Swiss Grid (Switzerland)"); }