Index: test/unit/org/openstreetmap/josm/data/projection/SwissGridTest.java =================================================================== --- test/unit/org/openstreetmap/josm/data/projection/SwissGridTest.java (révision 0) +++ test/unit/org/openstreetmap/josm/data/projection/SwissGridTest.java (révision 0) @@ -0,0 +1,157 @@ +// License: GPL. For details, see LICENSE file. +package org.openstreetmap.josm.data.projection; + +import static org.junit.Assert.assertTrue; + +import org.junit.BeforeClass; +import org.junit.Test; +import org.openstreetmap.josm.Main; +import org.openstreetmap.josm.data.coor.EastNorth; +import org.openstreetmap.josm.data.coor.LatLon; + +public class SwissGridTest { + @BeforeClass + public static void setUp() { + Main.proj = new SwissGrid(); + } + + @Test + public void a_latlon2eastNorth_test() { + { + LatLon ll = new LatLon(46.518, 6.567); + EastNorth en = Main.proj.latlon2eastNorth(ll); + System.out.println(en); + assertTrue("Lausanne", Math.abs(en.east() - 533111.69) < 0.1); + assertTrue("Lausanne", Math.abs(en.north() - 152227.85) < 0.1); + } + + { + LatLon ll = new LatLon(47.78, 8.58); + EastNorth en = Main.proj.latlon2eastNorth(ll); + System.out.println(en); + assertTrue("Schafouse", Math.abs(en.east() - 685544.16) < 0.1); + assertTrue("Schafouse", Math.abs(en.north() - 292782.91) < 0.1); + } + + { + LatLon ll = new LatLon(46.58, 10.48); + EastNorth en = Main.proj.latlon2eastNorth(ll); + System.out.println(en); + assertTrue("Grinson", Math.abs(en.east() - 833068.04) < 0.1); + assertTrue("Grinson", Math.abs(en.north() - 163265.39) < 0.1); + } + { + LatLon ll = new LatLon(46.0 + 57.0 / 60 + 3.89813884505 / 3600, 7.0 + 26.0 / 60 + 19.076595154147 / 3600); + EastNorth en = Main.proj.latlon2eastNorth(ll); + System.out.println(en); + assertTrue("Berne", Math.abs(en.east() - 600000.0) < 0.1); + assertTrue("Berne", Math.abs(en.north() - 200000.0) < 0.1); + } + { + LatLon ll = new LatLon(46.044130, 8.730497); + EastNorth en = Main.proj.latlon2eastNorth(ll); + System.out.println(en); + assertTrue("Ref", Math.abs(en.east() - 700000.0) < 0.1); + assertTrue("Ref", Math.abs(en.north() - 100000.0) < 0.1); + } + } + + @Test + public void b_eastNorth2latlon_test() { + { + EastNorth en = new EastNorth(533111.69, 152227.85); + LatLon ll = Main.proj.eastNorth2latlon(en); + System.out.println(ll); + assertTrue("Lausanne", Math.abs(ll.lat() - 46.518) < 0.00001); + assertTrue("Lausanne", Math.abs(ll.lon() - 6.567) < 0.00001); + } + + { + EastNorth en = new EastNorth(685544.16, 292782.91); + LatLon ll = Main.proj.eastNorth2latlon(en); + System.out.println(ll); + assertTrue("Schafouse", Math.abs(ll.lat() - 47.78) < 0.00001); + assertTrue("Schafouse", Math.abs(ll.lon() - 8.58) < 0.00001); + } + + { + EastNorth en = new EastNorth(833068.04, 163265.39); + LatLon ll = Main.proj.eastNorth2latlon(en); + System.out.println(ll); + assertTrue("Ref", Math.abs(ll.lat() - 46.58) < 0.00001); + assertTrue("Ref", Math.abs(ll.lon() -10.48) < 0.00001); + } + + { + EastNorth en = new EastNorth(600000.0, 200000.0); + LatLon ll = Main.proj.eastNorth2latlon(en); + System.out.println(ll); + assertTrue("Berne", Math.abs(ll.lat() - (46.0 + 57.0 / 60 + 3.89813884505 / 3600)) < 0.00001); + assertTrue("Berne", Math.abs(ll.lon() - (7.0 + 26.0 / 60 + 19.076595154147 / 3600)) < 0.00001); + } + + { + EastNorth en = new EastNorth(700000.0, 100000.0); + LatLon ll = Main.proj.eastNorth2latlon(en); + System.out.println(ll); + assertTrue("Ref", Math.abs(ll.lat() - 46.044130) < 0.00001); + assertTrue("Ref", Math.abs(ll.lon() - 8.730497) < 0.00001); + } + } + + /** + * Send and return should have less than 1mm of difference. + */ + @Test + public void c_sendandreturn_test() { + { + EastNorth en = new EastNorth(533111.69, 152227.85); + LatLon ll = Main.proj.eastNorth2latlon(en); + EastNorth en2 = Main.proj.latlon2eastNorth(ll); + System.out.println(en.east() - en2.east()); + System.out.println(en.north() - en2.north()); + assertTrue("Lausanne", Math.abs(en.east() - en2.east()) < 0.002); + assertTrue("Lausanne", Math.abs(en.north() - en2.north()) < 0.002); + } + + { + EastNorth en = new EastNorth(685544.16, 292782.91); + LatLon ll = Main.proj.eastNorth2latlon(en); + EastNorth en2 = Main.proj.latlon2eastNorth(ll); + System.out.println(en.east() - en2.east()); + System.out.println(en.north() - en2.north()); + assertTrue("Schafouse", Math.abs(en.east() - en2.east()) < 0.002); + assertTrue("Schafouse", Math.abs(en.north() - en2.north()) < 0.002); + } + + { + EastNorth en = new EastNorth(833068.04, 163265.39); + LatLon ll = Main.proj.eastNorth2latlon(en); + EastNorth en2 = Main.proj.latlon2eastNorth(ll); + System.out.println(en.east() - en2.east()); + System.out.println(en.north() - en2.north()); + assertTrue("Ref", Math.abs(en.east() - en2.east()) < 0.002); + assertTrue("Ref", Math.abs(en.north() - en2.north()) < 0.002); + } + + { + EastNorth en = new EastNorth(600000.0, 200000.0); + LatLon ll = Main.proj.eastNorth2latlon(en); + EastNorth en2 = Main.proj.latlon2eastNorth(ll); + System.out.println(en.east() - en2.east()); + System.out.println(en.north() - en2.north()); + assertTrue("Berne", Math.abs(en.east() - en2.east()) < 0.002); + assertTrue("Berne", Math.abs(en.north() - en2.north()) < 0.002); + } + + { + EastNorth en = new EastNorth(700000.0, 100000.0); + LatLon ll = Main.proj.eastNorth2latlon(en); + EastNorth en2 = Main.proj.latlon2eastNorth(ll); + System.out.println(en.east() - en2.east()); + System.out.println(en.north() - en2.north()); + assertTrue("Ref", Math.abs(en.east() - en2.east()) < 0.002); + assertTrue("Ref", Math.abs(en.north() - en2.north()) < 0.002); + } + } +} Index: src/org/openstreetmap/josm/data/projection/SwissGrid.java =================================================================== --- src/org/openstreetmap/josm/data/projection/SwissGrid.java (révision 3470) +++ src/org/openstreetmap/josm/data/projection/SwissGrid.java (copie de travail) @@ -14,77 +14,107 @@ * 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 DELTA_PHI = 1e-8; + private static final double a = 6377397.155; + private static final double b = 6356078.962822; + private static final Ellipsoid SWISS_ELLIPSOID = new Ellipsoid(a, b); + 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 - SWISS_ELLIPSOID.e2) / (1 - (SWISS_ELLIPSOID.e2 * Math.pow(Math.sin(phi0), 2))); + private static final double alpha = Math.sqrt(1 + (SWISS_ELLIPSOID.eb2 * 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 * SWISS_ELLIPSOID.e / 2 + * Math.log((1 + SWISS_ELLIPSOID.e * Math.sin(phi0)) / (1 - SWISS_ELLIPSOID.e * Math.sin(phi0))); + + private static final double xTrans = 200000; + private static final double yTrans = 600000; + + private LatLon correctETRS89toCH1903(LatLon coord) { + double[] XYZ = Ellipsoid.WGS84.latLon2Cart(coord); + XYZ[0] -= dX; + XYZ[1] -= dY; + XYZ[2] -= dZ; + return SWISS_ELLIPSOID.cart2LatLon(XYZ, DELTA_PHI); + } + + private LatLon correctCH1903toETRS89(LatLon coord) { + double[] XYZ = SWISS_ELLIPSOID.latLon2Cart(coord); + XYZ[0] += dX; + XYZ[1] += dY; + XYZ[2] += dZ; + return Ellipsoid.WGS84.cart2LatLon(XYZ, DELTA_PHI); + } + /** * @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(new LatLon(Math.toRadians(wgs.lat()), Math.toRadians(wgs.lon()))); + double phi = coord.lat(); + double lambda = 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 * SWISS_ELLIPSOID.e / 2 + * Math.log((1 + SWISS_ELLIPSOID.e * Math.sin(phi)) / (1 - SWISS_ELLIPSOID.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.atan2(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.atan2(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; + double prevPhi = -1000; + int iteration = 0; + // iteration to finds S and phi + while (Math.abs(phi - prevPhi) > DELTA_PHI) { + if (++iteration > 20) { + throw new RuntimeException("Two many iterations"); + } + prevPhi = phi; + S = 1 / alpha * (Math.log(Math.tan(Math.PI / 4 + b / 2)) - K) + SWISS_ELLIPSOID.e + * Math.log(Math.tan(Math.PI / 4 + Math.asin(SWISS_ELLIPSOID.e * Math.sin(phi)) / 2)); + phi = 2 * Math.atan(Math.exp(S)) - Math.PI / 2; + } - return new LatLon(phi,lmb); + LatLon coord = correctCH1903toETRS89(new LatLon(phi, lambda)); + return new LatLon(Math.toDegrees(coord.lat()), Math.toDegrees(coord.lon())); } - @Override public String toString() { + @Override + public String toString() { return tr("Swiss Grid (Switzerland)"); } @@ -101,11 +131,8 @@ return "swissgrid"; } - public Bounds getWorldBoundsLatLon() - { - return new Bounds( - new LatLon(45.7, 5.7), - new LatLon(47.9, 10.6)); + public Bounds getWorldBoundsLatLon() { + return new Bounds(new LatLon(45.7, 5.7), new LatLon(47.9, 10.6)); } public double getDefaultZoomInPPD() { Index: src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java =================================================================== --- src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java (révision 3470) +++ src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java (copie de travail) @@ -347,22 +347,8 @@ * @param ell reference ellipsoid */ private LatLon cart2LatLon(double X, double Y, double Z, Ellipsoid ell) { - double norm = Math.sqrt(X * X + Y * Y); - double lg = 2.0 * Math.atan(Y / (X + norm)); - double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z))))); - double delta = 1.0; - while (delta > epsilon) { - double s2 = Math.sin(lt); - s2 *= s2; - double l = Math.atan((Z / norm) - / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2))))); - delta = Math.abs(l - lt); - lt = l; - } - double s2 = Math.sin(lt); - s2 *= s2; - // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2); - return new LatLon(lt, lg); + double[] XYZ = {X, Y, Z}; + return ell.cart2LatLon(XYZ, epsilon); } /** Index: src/org/openstreetmap/josm/data/projection/Ellipsoid.java =================================================================== --- src/org/openstreetmap/josm/data/projection/Ellipsoid.java (révision 3470) +++ src/org/openstreetmap/josm/data/projection/Ellipsoid.java (copie de travail) @@ -7,6 +7,8 @@ package org.openstreetmap.josm.data.projection; +import org.openstreetmap.josm.data.coor.LatLon; + /** * the reference ellipsoids */ @@ -163,4 +165,47 @@ } return lati1; } + + + /** + * Initializes from cartesian coordinates + * + * @param XYZ the coordinates in meters (X, Y, Z) + * @param epsilon the maximum delta we want on the phi angle + * @return The corresponding latitude and longitude in radian + */ + public LatLon cart2LatLon(double[] XYZ, double epsilon) { + double norm = Math.sqrt(XYZ[0] * XYZ[0] + XYZ[1] * XYZ[1]); + double lg = 2.0 * Math.atan(XYZ[1] / (XYZ[0] + norm)); + double lt = Math.atan(XYZ[2] / (norm * (1.0 - (a * e2 / Math.sqrt(XYZ[0] * XYZ[0] + XYZ[1] * XYZ[1] + XYZ[2] * XYZ[2]))))); + double delta = 1.0; + while (delta > epsilon) { + double s2 = Math.sin(lt); + s2 *= s2; + double l = Math.atan((XYZ[2] / norm) + / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2))))); + delta = Math.abs(l - lt); + lt = l; + } + return new LatLon(lt, lg); + } + + /** + * Initializes to cartesian coordinates + * + * @param coord The Latitude and longitude in radian + * @return the corresponding (X, Y Z) cartesian coordinates in meters. + */ + public double[] latLon2Cart(LatLon coord) { + double phi = coord.lat(); + double lambda = coord.lon(); + + double Rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2)); + double[] XYZ = new double[3]; + XYZ[0] = Rn * Math.cos(phi) * Math.cos(lambda); + XYZ[1] = Rn * Math.cos(phi) * Math.sin(lambda); + XYZ[2] = Rn * (1 - e2) * Math.sin(phi); + + return XYZ; + } }