Ticket #5327: swissgrid.3.patch
File swissgrid.3.patch, 17.9 KB (added by , 14 years ago) |
---|
-
test/unit/org/openstreetmap/josm/data/projection/SwissGridTest.java
1 // License: GPL. For details, see LICENSE file. 2 package org.openstreetmap.josm.data.projection; 3 4 import static org.junit.Assert.assertTrue; 5 6 import org.junit.BeforeClass; 7 import org.junit.Test; 8 import org.openstreetmap.josm.Main; 9 import org.openstreetmap.josm.data.coor.EastNorth; 10 import org.openstreetmap.josm.data.coor.LatLon; 11 12 public class SwissGridTest { 13 @BeforeClass 14 public static void setUp() { 15 Main.proj = new SwissGrid(); 16 } 17 18 @Test 19 public void a_latlon2eastNorth_test() { 20 { 21 LatLon ll = new LatLon(46.518, 6.567); 22 EastNorth en = Main.proj.latlon2eastNorth(ll); 23 System.out.println(en); 24 assertTrue("Lausanne", Math.abs(en.east() - 533111.69) < 0.1); 25 assertTrue("Lausanne", Math.abs(en.north() - 152227.85) < 0.1); 26 } 27 28 { 29 LatLon ll = new LatLon(47.78, 8.58); 30 EastNorth en = Main.proj.latlon2eastNorth(ll); 31 System.out.println(en); 32 assertTrue("Schafouse", Math.abs(en.east() - 685544.16) < 0.1); 33 assertTrue("Schafouse", Math.abs(en.north() - 292782.91) < 0.1); 34 } 35 36 { 37 LatLon ll = new LatLon(46.58, 10.48); 38 EastNorth en = Main.proj.latlon2eastNorth(ll); 39 System.out.println(en); 40 assertTrue("Grinson", Math.abs(en.east() - 833068.04) < 0.1); 41 assertTrue("Grinson", Math.abs(en.north() - 163265.39) < 0.1); 42 } 43 44 { 45 LatLon ll = new LatLon(46.0 + 57.0 / 60 + 3.89813884505 / 3600, 7.0 + 26.0 / 60 + 19.076595154147 / 3600); 46 EastNorth en = Main.proj.latlon2eastNorth(ll); 47 System.out.println(en); 48 assertTrue("Berne", Math.abs(en.east() - 600000.0) < 0.1); 49 assertTrue("Berne", Math.abs(en.north() - 200000.0) < 0.1); 50 } 51 { 52 LatLon ll = new LatLon(46.0 + 2.0 / 60 + 38.87 / 3600, 8.0 + 43.0 / 60 + 49.79 / 3600); 53 EastNorth en = Main.proj.latlon2eastNorth(ll); 54 System.out.println(en); 55 assertTrue("Ref", Math.abs(en.east() - 700000.0) < 0.1); 56 assertTrue("Ref", Math.abs(en.north() - 100000.0) < 0.1); 57 } 58 } 59 60 @Test 61 public void b_eastNorth2latlon_test() { 62 { 63 EastNorth en = new EastNorth(533111.69, 152227.85); 64 LatLon ll = Main.proj.eastNorth2latlon(en); 65 System.out.println(ll); 66 assertTrue("Lausanne", Math.abs(ll.lat() - 46.518) < 0.00001); 67 assertTrue("Lausanne", Math.abs(ll.lon() - 6.567) < 0.00001); 68 } 69 70 { 71 EastNorth en = new EastNorth(685544.16, 292782.91); 72 LatLon ll = Main.proj.eastNorth2latlon(en); 73 System.out.println(ll); 74 assertTrue("Schafouse", Math.abs(ll.lat() - 47.78) < 0.00001); 75 assertTrue("Schafouse", Math.abs(ll.lon() - 8.58) < 0.00001); 76 } 77 78 { 79 EastNorth en = new EastNorth(833068.04, 163265.39); 80 LatLon ll = Main.proj.eastNorth2latlon(en); 81 System.out.println(ll); 82 assertTrue("Grinson", Math.abs(ll.lat() - 46.58) < 0.00001); 83 assertTrue("Grinson", Math.abs(ll.lon() - 10.48) < 0.00001); 84 } 85 86 { 87 EastNorth en = new EastNorth(600000.0, 200000.0); 88 LatLon ll = Main.proj.eastNorth2latlon(en); 89 System.out.println(ll); 90 assertTrue("Berne", Math.abs(ll.lat() - (46.0 + 57.0 / 60 + 3.89813884505 / 3600)) < 0.00001); 91 assertTrue("Berne", Math.abs(ll.lon() - (7.0 + 26.0 / 60 + 19.076595154147 / 3600)) < 0.00001); 92 } 93 94 { 95 EastNorth en = new EastNorth(700000.0, 100000.0); 96 LatLon ll = Main.proj.eastNorth2latlon(en); 97 System.out.println(ll); 98 assertTrue("Ref", Math.abs(ll.lat() - 46.0 + 2.0 / 60 + 38.87 / 3600) < 0.00001); 99 assertTrue("Ref", Math.abs(ll.lon() - 8.0 + 43.0 / 60 + 49.79 / 3600) < 0.00001); 100 } 101 } 102 103 /** 104 * Send and return should have less than 2mm of difference. 105 */ 106 @Test 107 public void c_sendandreturn_test() { 108 { 109 EastNorth en = new EastNorth(533111.69, 152227.85); 110 LatLon ll = Main.proj.eastNorth2latlon(en); 111 EastNorth en2 = Main.proj.latlon2eastNorth(ll); 112 System.out.println(en.east() - en2.east()); 113 System.out.println(en.north() - en2.north()); 114 assertTrue("Lausanne", Math.abs(en.east() - en2.east()) < 0.002); 115 assertTrue("Lausanne", Math.abs(en.north() - en2.north()) < 0.002); 116 } 117 118 { 119 EastNorth en = new EastNorth(685544.16, 292782.91); 120 LatLon ll = Main.proj.eastNorth2latlon(en); 121 EastNorth en2 = Main.proj.latlon2eastNorth(ll); 122 System.out.println(en.east() - en2.east()); 123 System.out.println(en.north() - en2.north()); 124 assertTrue("Schafouse", Math.abs(en.east() - en2.east()) < 0.002); 125 assertTrue("Schafouse", Math.abs(en.north() - en2.north()) < 0.002); 126 } 127 128 { 129 EastNorth en = new EastNorth(833068.04, 163265.39); 130 LatLon ll = Main.proj.eastNorth2latlon(en); 131 EastNorth en2 = Main.proj.latlon2eastNorth(ll); 132 System.out.println(en.east() - en2.east()); 133 System.out.println(en.north() - en2.north()); 134 assertTrue("Grinson", Math.abs(en.east() - en2.east()) < 0.002); 135 assertTrue("Grinson", Math.abs(en.north() - en2.north()) < 0.002); 136 } 137 138 { 139 EastNorth en = new EastNorth(600000.0, 200000.0); 140 LatLon ll = Main.proj.eastNorth2latlon(en); 141 EastNorth en2 = Main.proj.latlon2eastNorth(ll); 142 System.out.println(en.east() - en2.east()); 143 System.out.println(en.north() - en2.north()); 144 assertTrue("Berne", Math.abs(en.east() - en2.east()) < 0.002); 145 assertTrue("Berne", Math.abs(en.north() - en2.north()) < 0.002); 146 } 147 148 { 149 EastNorth en = new EastNorth(700000.0, 100000.0); 150 LatLon ll = Main.proj.eastNorth2latlon(en); 151 EastNorth en2 = Main.proj.latlon2eastNorth(ll); 152 System.out.println(en.east() - en2.east()); 153 System.out.println(en.north() - en2.north()); 154 assertTrue("Ref", Math.abs(en.east() - en2.east()) < 0.002); 155 assertTrue("Ref", Math.abs(en.north() - en2.north()) < 0.002); 156 } 157 } 158 } -
src/org/openstreetmap/josm/data/projection/SwissGrid.java
14 14 * Calculations are based on formula from 15 15 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf 16 16 * 17 * August 2010 update to this formula (rigorous formulas) 18 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf 17 19 */ 18 20 public class SwissGrid implements Projection { 21 22 private static final double dX = 674.374; 23 private static final double dY = 15.056; 24 private static final double dZ = 405.346; 25 26 private static final double DELTA_PHI = 1e-8; 27 private static final double a = 6377397.155; 28 private static final double b = 6356078.962822; 29 private static final Ellipsoid BRESSEL1841 = new Ellipsoid(a, b); 30 private static final double phi0 = Math.toRadians(46.0 + 57.0 / 60 + 8.66 / 3600); 31 private static final double lambda0 = Math.toRadians(7.0 + 26.0 / 60 + 22.50 / 3600); 32 private static final double R = a * Math.sqrt(1 - BRESSEL1841.e2) / (1 - (BRESSEL1841.e2 * Math.pow(Math.sin(phi0), 2))); 33 private static final double alpha = Math.sqrt(1 + (BRESSEL1841.eb2 * Math.pow(Math.cos(phi0), 4))); 34 private static final double b0 = Math.asin(Math.sin(phi0) / alpha); 35 private static final double K = Math.log(Math.tan(Math.PI / 4 + b0 / 2)) - alpha 36 * Math.log(Math.tan(Math.PI / 4 + phi0 / 2)) + alpha * BRESSEL1841.e / 2 37 * Math.log((1 + BRESSEL1841.e * Math.sin(phi0)) / (1 - BRESSEL1841.e * Math.sin(phi0))); 38 39 private static final double xTrans = 200000; 40 private static final double yTrans = 600000; 41 42 private LatLon correctEllipoideGSR80toBressel1841(LatLon coord) { 43 double[] XYZ = Ellipsoid.WGS84.latLon2Cart(coord); 44 XYZ[0] -= dX; 45 XYZ[1] -= dY; 46 XYZ[2] -= dZ; 47 return BRESSEL1841.cart2LatLon(XYZ, DELTA_PHI); 48 } 49 50 private LatLon correctEllipoideBressel1841toGRS80(LatLon coord) { 51 double[] XYZ = BRESSEL1841.latLon2Cart(coord); 52 XYZ[0] += dX; 53 XYZ[1] += dY; 54 XYZ[2] += dZ; 55 return Ellipsoid.WGS84.cart2LatLon(XYZ, DELTA_PHI); 56 } 57 19 58 /** 20 59 * @param wgs WGS84 lat/lon (ellipsoid GRS80) (in degree) 21 60 * @return eastnorth projection in Swiss National Grid (ellipsoid Bessel) 22 61 */ 23 62 public EastNorth latlon2eastNorth(LatLon wgs) { 24 double phi = 3600d * wgs.lat(); 25 double lambda = 3600d * wgs.lon(); 63 LatLon coord = correctEllipoideGSR80toBressel1841(wgs); 64 double phi = Math.toRadians(coord.lat()); 65 double lambda = Math.toRadians(coord.lon()); 26 66 27 double phiprime = (phi - 169028.66d) / 10000d; 28 double lambdaprime = (lambda - 26782.5d) / 10000d; 67 double S = alpha * Math.log(Math.tan(Math.PI / 4 + phi / 2)) - alpha * BRESSEL1841.e / 2 68 * Math.log((1 + BRESSEL1841.e * Math.sin(phi)) / (1 - BRESSEL1841.e * Math.sin(phi))) + K; 69 double b = 2 * (Math.atan(Math.exp(S)) - Math.PI / 4); 70 double l = alpha * (lambda - lambda0); 29 71 30 // precompute squares for lambdaprime and phiprime 31 // 32 double lambdaprime_2 = Math.pow(lambdaprime,2); 33 double phiprime_2 = Math.pow(phiprime,2); 72 double lb = Math.atan2(Math.sin(l), Math.sin(b0) * Math.tan(b) + Math.cos(b0) * Math.cos(l)); 73 double bb = Math.asin(Math.cos(b0) * Math.sin(b) - Math.sin(b0) * Math.cos(b) * Math.cos(l)); 34 74 35 double north = 36 200147.07d 37 + 308807.95d * phiprime 38 + 3745.25d * lambdaprime_2 39 + 76.63d * phiprime_2 40 - 194.56d * lambdaprime_2 * phiprime 41 + 119.79d * Math.pow(phiprime,3); 75 double y = R * lb; 76 double x = R / 2 * Math.log((1 + Math.sin(bb)) / (1 - Math.sin(bb))); 42 77 43 double east = 44 600072.37d 45 + 211455.93d * lambdaprime 46 - 10938.51d * lambdaprime * phiprime 47 - 0.36d * lambdaprime * phiprime_2 48 - 44.54d * Math.pow(lambdaprime,3); 49 50 return new EastNorth(east, north); 78 return new EastNorth(y + yTrans, x + xTrans); 51 79 } 52 80 53 81 /** 54 82 * @param xy SwissGrid east/north (in meters) 55 83 * @return LatLon WGS84 (in degree) 56 84 */ 57 58 85 public LatLon eastNorth2latlon(EastNorth xy) { 59 double yp = (xy.east() - 600000d) / 1000000d;60 double xp = (xy.north() - 200000d) / 1000000d;86 double x = xy.north() - xTrans; 87 double y = xy.east() - yTrans; 61 88 62 // precompute squares of xp and yp 63 // 64 double xp_2 = Math.pow(xp,2); 65 double yp_2 = Math.pow(yp,2); 89 double lb = y / R; 90 double bb = 2 * (Math.atan(Math.exp(x / R)) - Math.PI / 4); 66 91 67 // lambda = latitude, phi = longitude 68 double lmbp = 2.6779094d 69 + 4.728982d * yp 70 + 0.791484d * yp * xp 71 + 0.1306d * yp * xp_2 72 - 0.0436d * Math.pow(yp,3); 92 double b = Math.asin(Math.cos(b0) * Math.sin(bb) + Math.sin(b0) * Math.cos(bb) * Math.cos(lb)); 93 double l = Math.atan2(Math.sin(lb), Math.cos(b0) * Math.cos(lb) - Math.sin(b0) * Math.tan(bb)); 73 94 74 double phip = 16.9023892d 75 + 3.238272d * xp 76 - 0.270978d * yp_2 77 - 0.002528d * xp_2 78 - 0.0447d * yp_2 * xp 79 - 0.0140d * Math.pow(xp,3); 95 double lambda = lambda0 + l / alpha; 96 double phi = b; 97 double S = 0; 80 98 81 double lmb = lmbp * 100.0d / 36.0d; 82 double phi = phip * 100.0d / 36.0d; 99 double prevPhi = -1000; 100 int iteration = 0; 101 // iteration to finds S and phi 102 while (Math.abs(phi - prevPhi) > DELTA_PHI) { 103 if (++iteration > 20) { 104 throw new RuntimeException("Two many iterations"); 105 } 106 prevPhi = phi; 107 S = 1 / alpha * (Math.log(Math.tan(Math.PI / 4 + b / 2)) - K) + BRESSEL1841.e 108 * Math.log(Math.tan(Math.PI / 4 + Math.asin(BRESSEL1841.e * Math.sin(phi)) / 2)); 109 phi = 2 * Math.atan(Math.exp(S)) - Math.PI / 2; 110 } 83 111 84 return new LatLon(phi,lmb); 112 LatLon coord = correctEllipoideBressel1841toGRS80(new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda))); 113 return coord; 85 114 } 86 115 87 @Override public String toString() { 116 @Override 117 public String toString() { 88 118 return tr("Swiss Grid (Switzerland)"); 89 119 } 90 120 … … 101 131 return "swissgrid"; 102 132 } 103 133 104 public Bounds getWorldBoundsLatLon() 105 { 106 return new Bounds( 107 new LatLon(45.7, 5.7), 108 new LatLon(47.9, 10.6)); 134 public Bounds getWorldBoundsLatLon() { 135 return new Bounds(new LatLon(45.7, 5.7), new LatLon(47.9, 10.6)); 109 136 } 110 137 111 138 public double getDefaultZoomInPPD() { -
src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java
347 347 * @param ell reference ellipsoid 348 348 */ 349 349 private LatLon cart2LatLon(double X, double Y, double Z, Ellipsoid ell) { 350 double norm = Math.sqrt(X * X + Y * Y); 351 double lg = 2.0 * Math.atan(Y / (X + norm)); 352 double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z))))); 353 double delta = 1.0; 354 while (delta > epsilon) { 355 double s2 = Math.sin(lt); 356 s2 *= s2; 357 double l = Math.atan((Z / norm) 358 / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2))))); 359 delta = Math.abs(l - lt); 360 lt = l; 361 } 362 double s2 = Math.sin(lt); 363 s2 *= s2; 364 // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2); 365 return new LatLon(lt, lg); 350 double[] XYZ = {X, Y, Z}; 351 LatLon coord = ell.cart2LatLon(XYZ, epsilon); 352 return new LatLon(Math.toRadians(coord.lat()), Math.toRadians(coord.lon())); 366 353 } 367 354 368 355 /** -
src/org/openstreetmap/josm/data/projection/Ellipsoid.java
7 7 8 8 package org.openstreetmap.josm.data.projection; 9 9 10 import org.openstreetmap.josm.data.coor.LatLon; 11 10 12 /** 11 13 * the reference ellipsoids 12 14 */ … … 163 165 } 164 166 return lati1; 165 167 } 168 169 170 /** 171 * Initializes from cartesian coordinates 172 * 173 * @param XYZ the coordinates in meters (X, Y, Z) 174 * @param epsilon the maximum delta we want on the phi angle 175 * @return The corresponding latitude and longitude in degrees 176 */ 177 public LatLon cart2LatLon(double[] XYZ, double epsilon) { 178 double norm = Math.sqrt(XYZ[0] * XYZ[0] + XYZ[1] * XYZ[1]); 179 double lg = 2.0 * Math.atan(XYZ[1] / (XYZ[0] + norm)); 180 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]))))); 181 double delta = 1.0; 182 while (delta > epsilon) { 183 double s2 = Math.sin(lt); 184 s2 *= s2; 185 double l = Math.atan((XYZ[2] / norm) 186 / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2))))); 187 delta = Math.abs(l - lt); 188 lt = l; 189 } 190 return new LatLon(Math.toDegrees(lt), Math.toDegrees(lg)); 191 } 192 193 /** 194 * Initializes to cartesian coordinates 195 * 196 * @param coord The Latitude and longitude in degrees 197 * @return the corresponding (X, Y Z) cartesian coordinates in meters. 198 */ 199 public double[] latLon2Cart(LatLon coord) { 200 double phi = Math.toRadians(coord.lat()); 201 double lambda = Math.toRadians(coord.lon()); 202 203 double Rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2)); 204 double[] XYZ = new double[3]; 205 XYZ[0] = Rn * Math.cos(phi) * Math.cos(lambda); 206 XYZ[1] = Rn * Math.cos(phi) * Math.sin(lambda); 207 XYZ[2] = Rn * (1 - e2) * Math.sin(phi); 208 209 return XYZ; 210 } 166 211 }