Ticket #2540: swissgrid.patch

File swissgrid.patch, 6.5 KB (added by Gubaer, 4 years ago)
  • Projection.java

     
    3434        new Epsg4326(), 
    3535        new Mercator(), 
    3636        new Lambert(), 
    37         new LambertEST() 
     37        new LambertEST(), 
     38        new SwissGrid() 
    3839    }; 
    3940 
    4041    /** 
  • SwissGrid.java

     
     1//License: GPL. For details, see LICENSE file. 
     2 
     3package org.openstreetmap.josm.data.projection; 
     4 
     5import static org.openstreetmap.josm.tools.I18n.tr; 
     6 
     7import java.util.logging.Logger; 
     8 
     9import javax.swing.JOptionPane; 
     10 
     11import org.openstreetmap.josm.Main; 
     12import org.openstreetmap.josm.data.coor.EastNorth; 
     13import org.openstreetmap.josm.data.coor.LatLon; 
     14 
     15 
     16/** 
     17 * Projection for the SwissGrid, see http://de.wikipedia.org/wiki/Swiss_Grid. 
     18 *  
     19 * Calculations are based on formula from  
     20 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf 
     21 * 
     22 */ 
     23public class SwissGrid implements Projection { 
     24 
     25    static private Logger logger = Logger.getLogger(SwissGrid.class.getName()); 
     26    private boolean doAlertOnCoordinatesOufOfRange = true;  
     27     
     28     
     29    /** 
     30     * replies true if if wgs is in or reasonably close to Switzerland. False otherwise. 
     31     *  
     32     * @param wgs  lat/lon in WGS89 
     33     * @return 
     34     */ 
     35    protected boolean latlonInAcceptableRange(LatLon wgs) { 
     36         
     37        // coordinate transformation is invoked for boundary values regardless 
     38        // of current data set.   
     39        // 
     40        if (Math.abs(wgs.lon()) == Projection.MAX_LON && Math.abs(wgs.lat()) == Projection.MAX_LAT) { 
     41            return true; 
     42        } 
     43        return   wgs.lon() >= 5.7 && wgs.lon() <= 10.6 
     44               && wgs.lat() >= 45.7 && wgs.lat() <= 47.9; 
     45    } 
     46     
     47    /** 
     48     * displays an alert if lat/lon are not reasonably close to Switzerland. 
     49     */ 
     50    protected void alertCoordinatesOutOfRange() { 
     51        JOptionPane.showMessageDialog(Main.parent, 
     52                tr("The projection \"{0}\" is designed for\n" 
     53                + "latitudes between 45.7' and 47.9'\n" 
     54                + "and longitutes between 5.7' and 10.6' only.\n" 
     55                + "Use another projection system if you are not working\n" 
     56                + "on a data set of Switzerland or Lichtenstein.\n" 
     57                + "Do not upload any data after this message.", this.toString()), 
     58                "Current projection not suitable", 
     59                JOptionPane.WARNING_MESSAGE                 
     60        ); 
     61        doAlertOnCoordinatesOufOfRange = false;         
     62    } 
     63     
     64    /** 
     65     * @param wgs  WGS84 lat/lon (ellipsoid GRS80) (in degree) 
     66     * @return eastnorth projection in Swiss National Grid (ellipsoid Bessel) 
     67     */ 
     68    public EastNorth latlon2eastNorth(LatLon wgs) { 
     69 
     70 
     71            if (!latlonInAcceptableRange(wgs)) { 
     72                if (doAlertOnCoordinatesOufOfRange) { 
     73                    alertCoordinatesOutOfRange(); 
     74                } 
     75                logger.warning("lat/lon out of range for SwissGrid: lat/lon=" + wgs); 
     76            } 
     77             
     78            double phi = 3600d * wgs.lat(); 
     79            double lambda = 3600d * wgs.lon(); 
     80     
     81            double phiprime = (phi - 169028.66d) / 10000d; 
     82            double lambdaprime = (lambda - 26782.5d) / 10000d; 
     83 
     84            // precompute squares for lambdaprime and phiprime 
     85            // 
     86            double lambdaprime_2 = Math.pow(lambdaprime,2); 
     87            double phiprime_2 = Math.pow(phiprime,2); 
     88     
     89     
     90            double north = 
     91                  200147.07d 
     92                + 308807.95d                           * phiprime 
     93                +   3745.25d    * lambdaprime_2 
     94                +     76.63d                           * phiprime_2 
     95                -    194.56d    * lambdaprime_2        * phiprime                
     96                +    119.79d                           * Math.pow(phiprime,3); 
     97          
     98           
     99            double east = 
     100                  600072.37d 
     101                + 211455.93d  * lambdaprime 
     102                - 10938.51d   * lambdaprime             * phiprime 
     103                - 0.36d       * lambdaprime             * phiprime_2 
     104                - 44.54d      * Math.pow(lambdaprime,3); 
     105                     
     106            return new EastNorth(east, north); 
     107    } 
     108     
     109     
     110    /** 
     111     * @param xy SwissGrid east/north (in meters) 
     112     * @return LatLon WGS84 (in degree) 
     113     */ 
     114      
     115    public LatLon eastNorth2latlon(EastNorth xy) { 
     116    
     117        double yp = (xy.east() - 600000d) / 1000000d; 
     118        double xp = (xy.north() - 200000d) / 1000000d; 
     119         
     120        // precompute squares of xp and yp  
     121        // 
     122        double xp_2 = Math.pow(xp,2); 
     123        double yp_2 = Math.pow(yp,2); 
     124        
     125     
     126        // lambda = latitude, phi = longitude 
     127        double lmbp  =    2.6779094d             
     128                        + 4.728982d      * yp 
     129                        + 0.791484d      * yp               * xp  
     130                        + 0.1306d        * yp               * xp_2 
     131                        - 0.0436d        * Math.pow(yp,3); 
     132          
     133        double phip =     16.9023892d          
     134                        +  3.238272d                        * xp          
     135                        -  0.270978d    * yp_2       
     136                        -  0.002528d                        * xp_2 
     137                        -  0.0447d      * yp_2              * xp    
     138                        -  0.0140d                          * Math.pow(xp,3); 
     139         
     140 
     141          
     142        double lmb = lmbp * 100.0d / 36.0d; 
     143        double phi = phip * 100.0d / 36.0d; 
     144   
     145           
     146        return new LatLon(phi,lmb); 
     147         
     148    } 
     149 
     150    public double scaleFactor() { 
     151        return 1.0; 
     152    } 
     153     
     154    @Override public String toString() { 
     155        return tr("Swiss Grid (Switzerland)"); 
     156    } 
     157 
     158    public String toCode() { 
     159        return "EPSG:21781"; 
     160    } 
     161 
     162    public String getCacheDirectoryName() { 
     163        return "swissgrid"; 
     164    } 
     165 
     166 
     167    @Override 
     168    public boolean equals(Object o) { 
     169        return o instanceof SwissGrid; 
     170    } 
     171 
     172    @Override 
     173    public int hashCode() { 
     174        return SwissGrid.class.hashCode(); 
     175    } 
     176 
     177}