Ticket #5327: swiss.patch

File swiss.patch, 8.8 KB (added by sbrunner, 14 years ago)

Patch to use the rigorous formulas

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

     
    33package org.openstreetmap.josm.data.projection;
    44
    55import static org.openstreetmap.josm.tools.I18n.tr;
    6 
    76import org.openstreetmap.josm.data.Bounds;
    87import org.openstreetmap.josm.data.coor.EastNorth;
    98import org.openstreetmap.josm.data.coor.LatLon;
    … …  
    1413 * Calculations are based on formula from
    1514 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf
    1615 *
     16 * August 2010 update to this formula (rigorous formulas)
     17 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf
    1718 */
    1819public class SwissGrid implements Projection {
     20   
     21    private static final double dX = 674.374;
     22    private static final double dY = 15.056;
     23    private static final double dZ = 405.346;
     24   
     25    private static final double wgs84_a = 6378137.000;
     26    private static final double wgs84_b = 6356752.314245;
     27//    private static final double wgs84_f = (wgs84_a-wgs84_b)/wgs84_a;
     28    private static final double wgs84_E2 = (Math.pow(wgs84_a, 2)-Math.pow(wgs84_b, 2))/Math.pow(wgs84_a, 2);
     29//    private static final double wgs84_E = Math.sqrt(wgs84_E2);
     30
     31   
     32    private static final double a = 6377397.155;
     33    private static final double b = 6356078.962822;
     34//    private static final double f = (a-b)/a;
     35    private static final double E2 = (Math.pow(a, 2)-Math.pow(b, 2))/Math.pow(a, 2);
     36    private static final double E = Math.sqrt(E2);
     37    private static final double phi0 = Math.toRadians(46.0 + 57.0 / 60 + 8.66 / 3600);
     38    private static final double lambda0 = Math.toRadians(7.0 + 26.0 / 60 + 22.50 / 3600);
     39    private static final double R = a * Math.sqrt(1 - E2) / (1 - (E2 * Math.pow(Math.sin(phi0), 2)));
     40    private static final double alpha = Math.sqrt(1 + (E2 / (1 - E2) * Math.pow(Math.cos(phi0), 4)));
     41    private static final double b0 = Math.asin(Math.sin(phi0) / alpha);
     42    private static final double K = Math.log(Math.tan(Math.PI / 4 + b0 / 2)) - alpha
     43            * Math.log(Math.tan(Math.PI / 4 + phi0 / 2)) + alpha * E / 2
     44            * Math.log((1 + E * Math.sin(phi0)) / (1 - E * Math.sin(phi0)));
     45
     46    private static final double xTrans = 200000;
     47    private static final double yTrans = 600000;
     48
     49    private LatLon correctETRS89toCH1903(LatLon coord) {
     50        double phi = Math.toRadians(coord.lat());
     51        double lambda = Math.toRadians(coord.lon());
     52       
     53        double Rn = wgs84_a/Math.sqrt(1-wgs84_E2*Math.pow(Math.sin(phi), 2));
     54        double X = Rn*Math.cos(phi)*Math.cos(lambda);
     55        double Y = Rn*Math.cos(phi)*Math.sin(lambda);
     56        double Z = Rn*(1-wgs84_E2)*Math.sin(phi);
     57        X -= dX;
     58        Y -= dY;
     59        Z -= dZ;
     60        lambda = Math.atan(Y/X);
     61        double sX2Ys = Math.sqrt(X*X+Y*Y);
     62        phi = Math.atan((Z/sX2Ys));
     63        double h = 0;
     64       
     65        for (int i = 0 ; i < 6 ; i++) {
     66            Rn = a/Math.sqrt(1-E2*Math.pow(Math.sin(phi), 2));
     67            h = sX2Ys/Math.cos(phi)-Rn;
     68            phi = Math.atan((Z/sX2Ys)/(1-Rn*E2/(Rn+h)));
     69        }
     70
     71        return new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda));
     72    }
     73   
     74    private LatLon correctCH1903toETRS89(LatLon coord) {
     75        double phi = Math.toRadians(coord.lat());
     76        double lambda = Math.toRadians(coord.lon());
     77       
     78        double Rn = a/Math.sqrt(1-E2*Math.pow(Math.sin(phi), 2));
     79        double X = Rn*Math.cos(phi)*Math.cos(lambda);
     80        double Y = Rn*Math.cos(phi)*Math.sin(lambda);
     81        double Z = Rn*(1-E2)*Math.sin(phi);
     82        X += dX;
     83        Y += dY;
     84        Z += dZ;
     85        lambda = Math.atan(Y/X);
     86        double sX2Ys = Math.sqrt(X*X+Y*Y);
     87        phi = Math.atan((Z/sX2Ys));
     88        double h = 0;
     89       
     90        for (int i = 0 ; i < 5 ; i++) {
     91            Rn = wgs84_a/Math.sqrt(1-wgs84_E2*Math.pow(Math.sin(phi), 2));
     92            h = sX2Ys/Math.cos(phi)-Rn;
     93            phi = Math.atan((Z/sX2Ys)/(1-Rn*wgs84_E2/(Rn+h)));
     94        }
     95
     96        return new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda));
     97    }
     98   
    1999    /**
    20100     * @param wgs  WGS84 lat/lon (ellipsoid GRS80) (in degree)
    21101     * @return eastnorth projection in Swiss National Grid (ellipsoid Bessel)
    22102     */
    23103    public EastNorth latlon2eastNorth(LatLon wgs) {
    24             double phi = 3600d * wgs.lat();
    25             double lambda = 3600d * wgs.lon();
     104        LatLon coord = correctETRS89toCH1903(wgs);
     105        double phi = Math.toRadians(coord.lat());
     106        double lambda = Math.toRadians(coord.lon());
    26107
    27             double phiprime = (phi - 169028.66d) / 10000d;
    28             double lambdaprime = (lambda - 26782.5d) / 10000d;
     108        double S = alpha * Math.log(Math.tan(Math.PI / 4 + phi / 2)) - alpha * E / 2
     109                * Math.log((1 + E * Math.sin(phi)) / (1 - E * Math.sin(phi))) + K;
     110        double b = 2 * (Math.atan(Math.exp(S)) - Math.PI / 4);
     111        double l = alpha * (lambda - lambda0);
    29112
    30             // precompute squares for lambdaprime and phiprime
    31             //
    32             double lambdaprime_2 = Math.pow(lambdaprime,2);
    33             double phiprime_2 = Math.pow(phiprime,2);
     113        double lb = Math.atan(Math.sin(l) / (Math.sin(b0) * Math.tan(b) + Math.cos(b0) * Math.cos(l)));
     114        double bb = Math.asin(Math.cos(b0) * Math.sin(b) - Math.sin(b0) * Math.cos(b) * Math.cos(l));
    34115
    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);
     116        double y = R * lb;
     117        double x = R / 2 * Math.log((1 + Math.sin(bb)) / (1 - Math.sin(bb)));
    42118
    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);
     119        return new EastNorth(y + yTrans, x + xTrans);
    51120    }
    52121
    53122    /**
    54123     * @param xy SwissGrid east/north (in meters)
    55124     * @return LatLon WGS84 (in degree)
    56125     */
    57 
    58126    public LatLon eastNorth2latlon(EastNorth xy) {
    59         double yp = (xy.east() - 600000d) / 1000000d;
    60         double xp = (xy.north() - 200000d) / 1000000d;
     127        double x = xy.north() - xTrans;
     128        double y = xy.east() - yTrans;
    61129
    62         // precompute squares of xp and yp
    63         //
    64         double xp_2 = Math.pow(xp,2);
    65         double yp_2 = Math.pow(yp,2);
     130        double lb = y / R;
     131        double bb = 2 * (Math.atan(Math.exp(x / R)) - Math.PI / 4);
    66132
    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);
     133        double b = Math.asin(Math.cos(b0) * Math.sin(bb) + Math.sin(b0) * Math.cos(bb) * Math.cos(lb));
     134        double l = Math.atan(Math.sin(lb) / (Math.cos(b0) * Math.cos(lb) - Math.sin(b0) * Math.tan(bb)));
    73135
    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);
     136        double lambda = lambda0 + l / alpha;
     137        double phi = b;
     138        double S = 0;
    80139
    81         double lmb = lmbp * 100.0d / 36.0d;
    82         double phi = phip * 100.0d / 36.0d;
     140        // 5 iteration to fins S and phi
     141        for (int i = 0; i < 6; i++) {
     142            S = 1 / alpha * (Math.log(Math.tan(Math.PI / 4 + b / 2)) - K) + E
     143                    * Math.log(Math.tan(Math.PI / 4 + Math.asin(E * Math.sin(phi)) / 2));
     144            phi = 2 * Math.atan(Math.exp(S)) - Math.PI / 2;
     145        }
    83146
    84         return new LatLon(phi,lmb);
     147        LatLon coord = new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda));
     148        return correctCH1903toETRS89(coord);
    85149    }
    86150
    87     @Override public String toString() {
     151    @Override
     152    public String toString() {
    88153        return tr("Swiss Grid (Switzerland)");
    89154    }