Changeset 3473 in josm


Ignore:
Timestamp:
Aug 26, 2010 11:56:47 PM (3 years ago)
Author:
bastiK
Message:

see #5327 (patch by sbrunner) - Swissgrild uses approximate formulas => better use exact formulas

Location:
trunk
Files:
2 added
4 edited

Legend:

Unmodified
Added
Removed

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

    r3083 r3473  
    88package org.openstreetmap.josm.data.projection; 
    99 
     10import org.openstreetmap.josm.data.coor.LatLon; 
     11 
    1012/** 
    1113 * the reference ellipsoids 
    … …  
    3032     */ 
    3133    public static final Ellipsoid WGS84 = new Ellipsoid(6378137.0, 6356752.3142); 
     34 
     35    /** 
     36     * Bessel 1841 ellipsoid 
     37     */ 
     38    public static final Ellipsoid Bessel1841 = new Ellipsoid(6377397.155, 6356078.962822); 
    3239 
    3340    /** 
    … …  
    164171        return lati1; 
    165172    } 
     173 
     174    /** 
     175     * convert cartesian coordinates to ellipsoidal coordinates 
     176     * 
     177     * @param XYZ the coordinates in meters (X, Y, Z) 
     178     * @return The corresponding latitude and longitude in degrees 
     179     */ 
     180    public LatLon cart2LatLon(double[] XYZ) { 
     181        return cart2LatLon(XYZ, 1e-11); 
     182    } 
     183    public LatLon cart2LatLon(double[] XYZ, double epsilon) { 
     184        double norm = Math.sqrt(XYZ[0] * XYZ[0] + XYZ[1] * XYZ[1]); 
     185        double lg = 2.0 * Math.atan(XYZ[1] / (XYZ[0] + norm)); 
     186        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]))))); 
     187        double delta = 1.0; 
     188        while (delta > epsilon) { 
     189            double s2 = Math.sin(lt); 
     190            s2 *= s2; 
     191            double l = Math.atan((XYZ[2] / norm) 
     192                    / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2))))); 
     193            delta = Math.abs(l - lt); 
     194            lt = l; 
     195        } 
     196        return new LatLon(Math.toDegrees(lt), Math.toDegrees(lg)); 
     197    } 
     198     
     199    /** 
     200     * convert ellipsoidal coordinates to cartesian coordinates 
     201     * 
     202     * @param coord The Latitude and longitude in degrees 
     203     * @return the corresponding (X, Y Z) cartesian coordinates in meters. 
     204     */ 
     205    public double[] latLon2Cart(LatLon coord) { 
     206        double phi = Math.toRadians(coord.lat()); 
     207        double lambda = Math.toRadians(coord.lon()); 
     208 
     209        double Rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2)); 
     210        double[] XYZ = new double[3]; 
     211        XYZ[0] = Rn * Math.cos(phi) * Math.cos(lambda); 
     212        XYZ[1] = Rn * Math.cos(phi) * Math.sin(lambda); 
     213        XYZ[2] = Rn * (1 - e2) * Math.sin(phi); 
     214         
     215        return XYZ; 
     216    } 
    166217} 

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

    r3083 r3473  
    11//License: GPL. For details, see LICENSE file. 
    2  
    32package org.openstreetmap.josm.data.projection; 
    43 
    54import static org.openstreetmap.josm.tools.I18n.tr; 
    65 
     6import java.util.Collection; 
     7import java.util.Collections; 
     8import javax.swing.Box; 
     9import javax.swing.JPanel; 
     10 
    711import org.openstreetmap.josm.data.Bounds; 
    812import org.openstreetmap.josm.data.coor.EastNorth; 
    913import org.openstreetmap.josm.data.coor.LatLon; 
     14import org.openstreetmap.josm.gui.widgets.HtmlPanel; 
     15import org.openstreetmap.josm.tools.GBC; 
    1016 
    1117/** 
    12  * Projection for the SwissGrid, see http://de.wikipedia.org/wiki/Swiss_Grid. 
     18 * Projection for the SwissGrid CH1903 / L03, see http://de.wikipedia.org/wiki/Swiss_Grid. 
    1319 * 
    1420 * Calculations are based on formula from 
    1521 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf 
    1622 * 
     23 * August 2010 update to this formula (rigorous formulas) 
     24 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf 
    1725 */ 
    18 public class SwissGrid implements Projection { 
     26public class SwissGrid implements Projection, ProjectionSubPrefs { 
     27 
     28    private static final double dX = 674.374; 
     29    private static final double dY = 15.056; 
     30    private static final double dZ = 405.346; 
     31 
     32    private static final double phi0 = Math.toRadians(46.0 + 57.0 / 60 + 8.66 / 3600); 
     33    private static final double lambda0 = Math.toRadians(7.0 + 26.0 / 60 + 22.50 / 3600); 
     34    private static final double R = Ellipsoid.Bessel1841.a * Math.sqrt(1 - Ellipsoid.Bessel1841.e2) / (1 - (Ellipsoid.Bessel1841.e2 * Math.pow(Math.sin(phi0), 2))); 
     35    private static final double alpha = Math.sqrt(1 + (Ellipsoid.Bessel1841.eb2 * Math.pow(Math.cos(phi0), 4))); 
     36    private static final double b0 = Math.asin(Math.sin(phi0) / alpha); 
     37    private static final double K = Math.log(Math.tan(Math.PI / 4 + b0 / 2)) - alpha 
     38            * Math.log(Math.tan(Math.PI / 4 + phi0 / 2)) + alpha * Ellipsoid.Bessel1841.e / 2 
     39            * Math.log((1 + Ellipsoid.Bessel1841.e * Math.sin(phi0)) / (1 - Ellipsoid.Bessel1841.e * Math.sin(phi0))); 
     40 
     41    private static final double xTrans = 200000; 
     42    private static final double yTrans = 600000; 
     43 
     44    private static final double DELTA_PHI = 1e-11; 
     45 
     46    private LatLon correctEllipoideGSR80toBressel1841(LatLon coord) { 
     47        double[] XYZ = Ellipsoid.WGS84.latLon2Cart(coord); 
     48        XYZ[0] -= dX; 
     49        XYZ[1] -= dY; 
     50        XYZ[2] -= dZ; 
     51        return Ellipsoid.Bessel1841.cart2LatLon(XYZ); 
     52    } 
     53 
     54    private LatLon correctEllipoideBressel1841toGRS80(LatLon coord) { 
     55        double[] XYZ = Ellipsoid.Bessel1841.latLon2Cart(coord); 
     56        XYZ[0] += dX; 
     57        XYZ[1] += dY; 
     58        XYZ[2] += dZ; 
     59        return Ellipsoid.WGS84.cart2LatLon(XYZ); 
     60    } 
     61 
    1962    /** 
    2063     * @param wgs  WGS84 lat/lon (ellipsoid GRS80) (in degree) 
    2164     * @return eastnorth projection in Swiss National Grid (ellipsoid Bessel) 
    2265     */ 
     66    @Override 
    2367    public EastNorth latlon2eastNorth(LatLon wgs) { 
    24             double phi = 3600d * wgs.lat(); 
    25             double lambda = 3600d * wgs.lon(); 
     68        LatLon coord = correctEllipoideGSR80toBressel1841(wgs); 
     69        double phi = Math.toRadians(coord.lat()); 
     70        double lambda = Math.toRadians(coord.lon()); 
    2671 
    27             double phiprime = (phi - 169028.66d) / 10000d; 
    28             double lambdaprime = (lambda - 26782.5d) / 10000d; 
     72        double S = alpha * Math.log(Math.tan(Math.PI / 4 + phi / 2)) - alpha * Ellipsoid.Bessel1841.e / 2 
     73                * Math.log((1 + Ellipsoid.Bessel1841.e * Math.sin(phi)) / (1 - Ellipsoid.Bessel1841.e * Math.sin(phi))) + K; 
     74        double b = 2 * (Math.atan(Math.exp(S)) - Math.PI / 4); 
     75        double l = alpha * (lambda - lambda0); 
    2976 
    30             // precompute squares for lambdaprime and phiprime 
    31             // 
    32             double lambdaprime_2 = Math.pow(lambdaprime,2); 
    33             double phiprime_2 = Math.pow(phiprime,2); 
     77        double lb = Math.atan2(Math.sin(l), Math.sin(b0) * Math.tan(b) + Math.cos(b0) * Math.cos(l)); 
     78        double bb = Math.asin(Math.cos(b0) * Math.sin(b) - Math.sin(b0) * Math.cos(b) * Math.cos(l)); 
    3479 
    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); 
     80        double y = R * lb; 
     81        double x = R / 2 * Math.log((1 + Math.sin(bb)) / (1 - Math.sin(bb))); 
    4282 
    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); 
     83        return new EastNorth(y + yTrans, x + xTrans); 
    5184    } 
    5285 
    … …  
    5588     * @return LatLon WGS84 (in degree) 
    5689     */ 
     90    @Override 
     91    public LatLon eastNorth2latlon(EastNorth xy) { 
     92        double x = xy.north() - xTrans; 
     93        double y = xy.east() - yTrans; 
    5794 
    58     public LatLon eastNorth2latlon(EastNorth xy) { 
    59         double yp = (xy.east() - 600000d) / 1000000d; 
    60         double xp = (xy.north() - 200000d) / 1000000d; 
     95        double lb = y / R; 
     96        double bb = 2 * (Math.atan(Math.exp(x / R)) - Math.PI / 4); 
    6197 
    62         // precompute squares of xp and yp 
    63         // 
    64         double xp_2 = Math.pow(xp,2); 
    65         double yp_2 = Math.pow(yp,2); 
     98        double b = Math.asin(Math.cos(b0) * Math.sin(bb) + Math.sin(b0) * Math.cos(bb) * Math.cos(lb)); 
     99        double l = Math.atan2(Math.sin(lb), Math.cos(b0) * Math.cos(lb) - Math.sin(b0) * Math.tan(bb)); 
    66100 
    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); 
     101        double lambda = lambda0 + l / alpha; 
     102        double phi = b; 
     103        double S = 0; 
    73104 
    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); 
     105        double prevPhi = -1000; 
     106        int iteration = 0; 
     107        // iteration to finds S and phi 
     108        while (Math.abs(phi - prevPhi) > DELTA_PHI) { 
     109            if (++iteration > 30) { 
     110                throw new RuntimeException("Two many iterations"); 
     111            } 
     112            prevPhi = phi; 
     113            S = 1 / alpha * (Math.log(Math.tan(Math.PI / 4 + b / 2)) - K) + Ellipsoid.Bessel1841.e 
     114                    * Math.log(Math.tan(Math.PI / 4 + Math.asin(Ellipsoid.Bessel1841.e * Math.sin(phi)) / 2)); 
     115            phi = 2 * Math.atan(Math.exp(S)) - Math.PI / 2; 
     116        } 
    80117 
    81         double lmb = lmbp * 100.0d / 36.0d; 
    82         double phi = phip * 100.0d / 36.0d; 
    83  
    84         return new LatLon(phi,lmb); 
     118        LatLon coord = correctEllipoideBressel1841toGRS80(new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda))); 
     119        return coord; 
    85120    } 
    86121 
    87     @Override public String toString() { 
     122    @Override 
     123    public String toString() { 
    88124        return tr("Swiss Grid (Switzerland)"); 
    89125    } 
    90126 
     127    @Override 
    91128    public String toCode() { 
    92129        return "EPSG:21781"; 
    … …  
    98135    } 
    99136 
     137    @Override 
    100138    public String getCacheDirectoryName() { 
    101139        return "swissgrid"; 
    102140    } 
    103141 
    104     public Bounds getWorldBoundsLatLon() 
    105     { 
    106         return new Bounds( 
    107                 new LatLon(45.7, 5.7), 
    108                 new LatLon(47.9, 10.6)); 
     142    @Override 
     143    public Bounds getWorldBoundsLatLon() { 
     144        return new Bounds(new LatLon(45.7, 5.7), new LatLon(47.9, 10.6)); 
    109145    } 
    110146 
     147    @Override 
    111148    public double getDefaultZoomInPPD() { 
    112149        // This will set the scale bar to about 100 m 
    113150        return 1.01; 
    114151    } 
     152 
     153    @Override 
     154    public void setupPreferencePanel(JPanel p) { 
     155        p.add(new HtmlPanel("<i>CH1903 / LV03 (without local corrections)</i>"), GBC.eol().fill(GBC.HORIZONTAL)); 
     156        p.add(Box.createVerticalGlue(), GBC.eol().fill(GBC.BOTH)); 
     157    } 
     158 
     159    @Override 
     160    public void setPreferences(Collection<String> args) { 
     161    } 
     162 
     163    @Override 
     164    public Collection<String> getPreferences(JPanel p) { 
     165        return Collections.singletonList("CH1903"); 
     166    } 
     167 
     168    @Override 
     169    public Collection<String> getPreferencesFromCode(String code) { 
     170        if ("EPSG:21781".equals(code)) 
     171            return Collections.singletonList("CH1903"); 
     172        return null; 
     173    } 
    115174} 

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

    r3343 r3473  
    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        LatLon coord = ell.cart2LatLon(XYZ, epsilon); 
     352        return new LatLon(Math.toRadians(coord.lat()), Math.toRadians(coord.lon())); 
    366353    } 
    367354 

  • trunk/src/org/openstreetmap/josm/gui/preferences/ProjectionPreference.java

    r3461 r3473  
    2929import org.openstreetmap.josm.data.projection.ProjectionSubPrefs; 
    3030import org.openstreetmap.josm.gui.NavigatableComponent; 
    31 import org.openstreetmap.josm.gui.widgets.VerticallyScrollablePanel; 
    3231import org.openstreetmap.josm.tools.GBC; 
    3332 
    … …  
    108107     * This is the panel holding all projection preferences 
    109108     */ 
    110     private JPanel projPanel = new VerticallyScrollablePanel(); 
     109    private JPanel projPanel = new JPanel(new GridBagLayout()); 
    111110 
    112111    /** 
Note: See TracChangeset for help on using the changeset viewer.