Ticket #5327: swissgrid.2.patch

File swissgrid.2.patch, 17.7 KB (added by sbrunner, 14 years ago)

Manly use Ellipsoide.

  • test/unit/org/openstreetmap/josm/data/projection/SwissGridTest.java

     
     1// License: GPL. For details, see LICENSE file.
     2package org.openstreetmap.josm.data.projection;
     3
     4import static org.junit.Assert.assertTrue;
     5
     6import org.junit.BeforeClass;
     7import org.junit.Test;
     8import org.openstreetmap.josm.Main;
     9import org.openstreetmap.josm.data.coor.EastNorth;
     10import org.openstreetmap.josm.data.coor.LatLon;
     11
     12public 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            LatLon ll = new LatLon(46.0 + 57.0 / 60 + 3.89813884505 / 3600, 7.0 + 26.0 / 60 + 19.076595154147 / 3600);
     45            EastNorth en = Main.proj.latlon2eastNorth(ll);
     46            System.out.println(en);
     47            assertTrue("Berne", Math.abs(en.east() - 600000.0) < 0.1);
     48            assertTrue("Berne", Math.abs(en.north() - 200000.0) < 0.1);
     49        }
     50        {
     51            LatLon ll = new LatLon(46.044130, 8.730497);
     52            EastNorth en = Main.proj.latlon2eastNorth(ll);
     53            System.out.println(en);
     54            assertTrue("Ref", Math.abs(en.east() - 700000.0) < 0.1);
     55            assertTrue("Ref", Math.abs(en.north() - 100000.0) < 0.1);
     56        }
     57    }
     58
     59    @Test
     60    public void b_eastNorth2latlon_test() {
     61        {
     62            EastNorth en = new EastNorth(533111.69, 152227.85);
     63            LatLon ll = Main.proj.eastNorth2latlon(en);
     64            System.out.println(ll);
     65            assertTrue("Lausanne", Math.abs(ll.lat() - 46.518) < 0.00001);
     66            assertTrue("Lausanne", Math.abs(ll.lon() - 6.567) < 0.00001);
     67        }
     68
     69        {
     70            EastNorth en = new EastNorth(685544.16, 292782.91);
     71            LatLon ll = Main.proj.eastNorth2latlon(en);
     72            System.out.println(ll);
     73            assertTrue("Schafouse", Math.abs(ll.lat() - 47.78) < 0.00001);
     74            assertTrue("Schafouse", Math.abs(ll.lon() - 8.58) < 0.00001);
     75        }
     76
     77        {
     78            EastNorth en = new EastNorth(833068.04, 163265.39);
     79            LatLon ll = Main.proj.eastNorth2latlon(en);
     80            System.out.println(ll);
     81            assertTrue("Ref", Math.abs(ll.lat() - 46.58) < 0.00001);
     82            assertTrue("Ref", Math.abs(ll.lon() -10.48) < 0.00001);
     83        }
     84
     85        {
     86            EastNorth en = new EastNorth(600000.0, 200000.0);
     87            LatLon ll = Main.proj.eastNorth2latlon(en);
     88            System.out.println(ll);
     89            assertTrue("Berne", Math.abs(ll.lat() - (46.0 + 57.0 / 60 + 3.89813884505 / 3600)) < 0.00001);
     90            assertTrue("Berne", Math.abs(ll.lon() - (7.0 + 26.0 / 60 + 19.076595154147 / 3600)) < 0.00001);
     91        }
     92       
     93        {
     94            EastNorth en = new EastNorth(700000.0, 100000.0);
     95            LatLon ll = Main.proj.eastNorth2latlon(en);
     96            System.out.println(ll);
     97            assertTrue("Ref", Math.abs(ll.lat() - 46.044130) < 0.00001);
     98            assertTrue("Ref", Math.abs(ll.lon() - 8.730497) < 0.00001);
     99        }
     100    }
     101
     102    /**
     103     * Send and return should have less than 1mm of difference.
     104     */
     105    @Test
     106    public void c_sendandreturn_test() {
     107        {
     108            EastNorth en = new EastNorth(533111.69, 152227.85);
     109            LatLon ll = Main.proj.eastNorth2latlon(en);
     110            EastNorth en2 = Main.proj.latlon2eastNorth(ll);
     111            System.out.println(en.east() - en2.east());
     112            System.out.println(en.north() - en2.north());
     113            assertTrue("Lausanne", Math.abs(en.east() - en2.east()) < 0.002);
     114            assertTrue("Lausanne", Math.abs(en.north() - en2.north()) < 0.002);
     115        }
     116
     117        {
     118            EastNorth en = new EastNorth(685544.16, 292782.91);
     119            LatLon ll = Main.proj.eastNorth2latlon(en);
     120            EastNorth en2 = Main.proj.latlon2eastNorth(ll);
     121            System.out.println(en.east() - en2.east());
     122            System.out.println(en.north() - en2.north());
     123            assertTrue("Schafouse", Math.abs(en.east() - en2.east()) < 0.002);
     124            assertTrue("Schafouse", Math.abs(en.north() - en2.north()) < 0.002);
     125        }
     126
     127        {
     128            EastNorth en = new EastNorth(833068.04, 163265.39);
     129            LatLon ll = Main.proj.eastNorth2latlon(en);
     130            EastNorth en2 = Main.proj.latlon2eastNorth(ll);
     131            System.out.println(en.east() - en2.east());
     132            System.out.println(en.north() - en2.north());
     133            assertTrue("Ref", Math.abs(en.east() - en2.east()) < 0.002);
     134            assertTrue("Ref", Math.abs(en.north() - en2.north()) < 0.002);
     135        }
     136
     137        {
     138            EastNorth en = new EastNorth(600000.0, 200000.0);
     139            LatLon ll = Main.proj.eastNorth2latlon(en);
     140            EastNorth en2 = Main.proj.latlon2eastNorth(ll);
     141            System.out.println(en.east() - en2.east());
     142            System.out.println(en.north() - en2.north());
     143            assertTrue("Berne", Math.abs(en.east() - en2.east()) < 0.002);
     144            assertTrue("Berne", Math.abs(en.north() - en2.north()) < 0.002);
     145        }
     146       
     147        {
     148            EastNorth en = new EastNorth(700000.0, 100000.0);
     149            LatLon ll = Main.proj.eastNorth2latlon(en);
     150            EastNorth en2 = Main.proj.latlon2eastNorth(ll);
     151            System.out.println(en.east() - en2.east());
     152            System.out.println(en.north() - en2.north());
     153            assertTrue("Ref", Math.abs(en.east() - en2.east()) < 0.002);
     154            assertTrue("Ref", Math.abs(en.north() - en2.north()) < 0.002);
     155        }
     156    }
     157}

  • src/org/openstreetmap/josm/data/projection/SwissGrid.java

     
    1414 * Calculations are based on formula from
    1515 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf
    1616 *
     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
    1719 */
    1820public 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 SWISS_ELLIPSOID = 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 - SWISS_ELLIPSOID.e2) / (1 - (SWISS_ELLIPSOID.e2 * Math.pow(Math.sin(phi0), 2)));
     33    private static final double alpha = Math.sqrt(1 + (SWISS_ELLIPSOID.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 * SWISS_ELLIPSOID.e / 2
     37            * Math.log((1 + SWISS_ELLIPSOID.e * Math.sin(phi0)) / (1 - SWISS_ELLIPSOID.e * Math.sin(phi0)));
     38
     39    private static final double xTrans = 200000;
     40    private static final double yTrans = 600000;
     41
     42    private LatLon correctETRS89toCH1903(LatLon coord) {
     43        double[] XYZ = Ellipsoid.WGS84.latLon2Cart(coord);
     44        XYZ[0] -= dX;
     45        XYZ[1] -= dY;
     46        XYZ[2] -= dZ;
     47        return SWISS_ELLIPSOID.cart2LatLon(XYZ, DELTA_PHI);
     48    }
     49
     50    private LatLon correctCH1903toETRS89(LatLon coord) {
     51        double[] XYZ = SWISS_ELLIPSOID.latLon2Cart(coord);
     52        XYZ[0] += dX;
     53        XYZ[1] += dY;
     54        XYZ[2] += dZ;
     55        return Ellipsoid.WGS84.cart2LatLon(XYZ, DELTA_PHI);
     56    }
     57
    1958    /**
    2059     * @param wgs  WGS84 lat/lon (ellipsoid GRS80) (in degree)
    2160     * @return eastnorth projection in Swiss National Grid (ellipsoid Bessel)
    2261     */
    2362    public EastNorth latlon2eastNorth(LatLon wgs) {
    24             double phi = 3600d * wgs.lat();
    25             double lambda = 3600d * wgs.lon();
     63        LatLon coord = correctETRS89toCH1903(new LatLon(Math.toRadians(wgs.lat()), Math.toRadians(wgs.lon())));
     64        double phi = coord.lat();
     65        double lambda = coord.lon();
    2666
    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 * SWISS_ELLIPSOID.e / 2
     68                * Math.log((1 + SWISS_ELLIPSOID.e * Math.sin(phi)) / (1 - SWISS_ELLIPSOID.e * Math.sin(phi))) + K;
     69        double b = 2 * (Math.atan(Math.exp(S)) - Math.PI / 4);
     70        double l = alpha * (lambda - lambda0);
    2971
    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));
    3474
    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)));
    4277
    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);
    5179    }
    5280
    5381    /**
    5482     * @param xy SwissGrid east/north (in meters)
    5583     * @return LatLon WGS84 (in degree)
    5684     */
    57 
    5885    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;
    6188
    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);
    6691
    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));
    7394
    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;
    8098
    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) + SWISS_ELLIPSOID.e
     108                    * Math.log(Math.tan(Math.PI / 4 + Math.asin(SWISS_ELLIPSOID.e * Math.sin(phi)) / 2));
     109            phi = 2 * Math.atan(Math.exp(S)) - Math.PI / 2;
     110        }
    83111
    84         return new LatLon(phi,lmb);
     112        LatLon coord = correctCH1903toETRS89(new LatLon(phi, lambda));
     113        return new LatLon(Math.toDegrees(coord.lat()), Math.toDegrees(coord.lon()));
    85114    }
    86115
    87     @Override public String toString() {
     116    @Override
     117    public String toString() {
    88118        return tr("Swiss Grid (Switzerland)");
    89119    }
    90120
    … …  
    101131        return "swissgrid";
    102132    }
    103133
    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));
    109136    }
    110137
    111138    public double getDefaultZoomInPPD() {

  • src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java

     
    347347     * @param ell reference ellipsoid
    348348     */
    349349    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        return ell.cart2LatLon(XYZ, epsilon);
    366352    }
    367353
    368354    /**

  • src/org/openstreetmap/josm/data/projection/Ellipsoid.java

     
    77
    88package org.openstreetmap.josm.data.projection;
    99
     10import org.openstreetmap.josm.data.coor.LatLon;
     11
    1012/**
    1113 * the reference ellipsoids
    1214 */
    … …  
    163165        }
    164166        return lati1;
    165167    }
     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 radian
     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(lt, lg);
     191    }
     192   
     193    /**
     194     * Initializes to cartesian coordinates
     195     *
     196     * @param coord The Latitude and longitude in radian
     197     * @return the corresponding (X, Y Z) cartesian coordinates in meters.
     198     */
     199    public double[] latLon2Cart(LatLon coord) {
     200        double phi = coord.lat();
     201        double lambda = 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    }
    166211}