Ignore:
Timestamp:
2011-08-07T12:17:39+02:00 (8 years ago)
Author:
bastiK
Message:

major projection rework

More modular structure, inspired by Proj.4.

There are almost no semantic changes to the projection algorithms. Mostly factors of 'a' and 180/PI have been moved from one place to the other. In UTM_France_DOM, the ellipsoid conversion for the first 3 projections has been changed from hayford <-> GRS80 to hayford <-> WGS84.

Some redundant algorithms have been removed. In particular:

  • UTM_France_DOM used to have its own Transverse Mercator implementation. It is different from the implementation in TransverseMercator.java as it has another series expansion. For EPSG::2975, there are numeric differences on centimeter scale. However, the new data fits better with Proj.4 output.
  • Also removed are alternate implementations of LambertConformalConic. (They are all quite similar, though.)
Location:
trunk/src/org/openstreetmap/josm/data/projection
Files:
14 added
11 edited
1 moved

Legend:

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

    r3530 r4285  
    1919    public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
    2020    /**
    21      * Hayford's ellipsoid (ED50 system)
     21     * Hayford's ellipsoid 1909 (ED50 system)
     22     * Proj.4 code: intl
    2223     */
    2324    public static final Ellipsoid hayford =
  • trunk/src/org/openstreetmap/josm/data/projection/Epsg3008.java

    r3692 r4285  
    66import org.openstreetmap.josm.data.Bounds;
    77import org.openstreetmap.josm.data.coor.LatLon;
    8 import org.openstreetmap.josm.data.coor.EastNorth;
     8import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
     9import org.openstreetmap.josm.data.projection.proj.TransverseMercator;
    910
    1011/**
     
    1213 * http://spatialreference.org/ref/epsg/3008/
    1314 *
    14  * @author Hanno Hecker, based on the TransverseMercatorLV.java by Viesturs Zarins
     15 * @author Hanno Hecker
    1516 */
    16 public class Epsg3008 extends TransverseMercator {
     17public class Epsg3008 extends AbstractProjection {
    1718
    18     private final static double UTMScaleFactor = 1.0;
    19     private double UTMCentralMeridianRad;
    20     private double offsetEastMeters = 150000;
    21     private double offsetNorthMeters = 0;
    22 
    23     public Epsg3008()
    24     {
    25         UTMCentralMeridianRad = Math.toRadians(13.5);
     19    public Epsg3008() {
     20        ellps = Ellipsoid.GRS80;
     21        proj = new TransverseMercator(ellps);
     22        datum = GRS80Datum.INSTANCE;
     23        lon_0 = 13.5;
     24        x_0 = 150000;
    2625    }
    27 
    28     @Override public String toString() {
     26   
     27    @Override
     28    public String toString() {
    2929        return tr("SWEREF99 13 30 / EPSG:3008 (Sweden)");
    3030    }
    3131
    32     private int epsgCode() {
     32    @Override
     33    public Integer getEpsgCode() {
    3334        return 3008;
    34     }
    35 
    36     @Override
    37     public String toCode() {
    38         return "EPSG:"+ epsgCode();
    3935    }
    4036
     
    4440    }
    4541
     42    @Override
    4643    public String getCacheDirectoryName() {
    47         return "epsg"+ epsgCode();
     44        return "epsg"+ getEpsgCode();
    4845    }
    4946
     
    5552    }
    5653
    57     @Override
    58     public EastNorth latlon2eastNorth(LatLon p) {
    59         EastNorth a = mapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), UTMCentralMeridianRad);
    60         return new EastNorth(a.east() * UTMScaleFactor + offsetEastMeters, a.north() * UTMScaleFactor + offsetNorthMeters);
    61     }
    62 
    63     @Override
    64     public LatLon eastNorth2latlon(EastNorth p) {
    65         return mapXYToLatLon((p.east() - offsetEastMeters)/UTMScaleFactor, (p.north() - offsetNorthMeters)/UTMScaleFactor, UTMCentralMeridianRad);
    66     }
    67 
    6854}
  • trunk/src/org/openstreetmap/josm/data/projection/Lambert.java

    r3779 r4285  
    66import java.awt.GridBagLayout;
    77import java.awt.event.ActionListener;
    8 import java.io.IOException;
    98import java.io.InputStream;
    109import java.util.Collection;
     
    1716import org.openstreetmap.josm.Main;
    1817import org.openstreetmap.josm.data.Bounds;
    19 import org.openstreetmap.josm.data.coor.EastNorth;
    2018import org.openstreetmap.josm.data.coor.LatLon;
     19import org.openstreetmap.josm.data.projection.proj.LambertConformalConic;
    2120import org.openstreetmap.josm.tools.GBC;
    2221import org.openstreetmap.josm.tools.ImageProvider;
    2322
    2423/**
    25  * This class provides the two methods <code>latlon2eastNorth</code> and <code>eastNorth2latlon</code>
    26  * converting the JOSM LatLon coordinates in WGS84 system (GPS) to and from East North values in
    27  * the projection Lambert conic conform 4 zones using the French geodetic system NTF.
     24 * Lambert conic conform 4 zones using the French geodetic system NTF.
    2825 * This newer version uses the grid translation NTF<->RGF93 provided by IGN for a submillimetric accuracy.
    2926 * (RGF93 is the French geodetic system similar to WGS84 but not mathematically equal)
    3027 * @author Pieren
    3128 */
    32 public class Lambert implements Projection, ProjectionSubPrefs {
     29public class Lambert extends AbstractProjection implements ProjectionSubPrefs {
    3330    /**
    3431     * Lambert I, II, III, and IV projection exponents
    3532     */
    36     public static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 };
     33    private static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 };
    3734
    3835    /**
    3936     * Lambert I, II, III, and IV projection constants
    4037     */
    41     public static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 };
     38    private static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 };
    4239
    4340    /**
    4441     * Lambert I, II, III, and IV false east
    4542     */
    46     public static final double Xs[] = { 600000.0, 600000.0, 600000.0, 234.358 };
     43    private static final double x_0s[] = { 600000.0, 600000.0, 600000.0, 234.358 };
    4744
    4845    /**
    4946     * Lambert I, II, III, and IV false north
    5047     */
    51     public static final double Ys[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 };
    52 
    53     /**
    54      * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian
    55      */
    56     public static final double lg0 = 0.04079234433198; // 2deg20'14.025"
    57 
    58     /**
    59      * precision in iterative schema
    60      */
    61 
    62     public static final double epsilon = 1e-11;
     48    private static final double y_fs[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 };
    6349
    6450    /**
     
    6854    public static final double cMinLatZone1Radian = Math.toRadians(46.1 * 0.9);// lowest latitude of Zone 4 (South Corsica)
    6955
    70     public static final double zoneLimitsDegree[][] = {
     56    public static final double[][] zoneLimitsDegree = {
    7157        {Math.toDegrees(cMaxLatZone1Radian), (53.5 * 0.9)}, // Zone 1  (reference values in grad *0.9)
    7258        {(53.5 * 0.9), (50.5 * 0.9)}, // Zone 2
     
    8672    public static final int DEFAULT_ZONE = 0;
    8773
    88     private int layoutZone = DEFAULT_ZONE;
     74    private int layoutZone;
    8975
    9076    private static NTV2GridShiftFile ntf_rgf93Grid = null;
     
    10086                InputStream is = Main.class.getResourceAsStream("/data/"+gridFileName);
    10187                if (is == null) {
    102                     System.err.println(tr("Warning: failed to open input stream for resource ''/data/{0}''. Cannot load NTF<->RGF93 grid", gridFileName));
    103                     return;
     88                    throw new RuntimeException(tr("Warning: failed to open input stream for resource ''/data/{0}''. Cannot load NTF<->RGF93 grid", gridFileName));
    10489                }
    10590                ntf_rgf93Grid = new NTV2GridShiftFile();
    10691                ntf_rgf93Grid.loadGridShiftFile(is, false);
    107                 //System.out.println("NTF<->RGF93 grid loaded.");
    10892            } catch (Exception e) {
    109                 e.printStackTrace();
     93                throw new RuntimeException(e);
    11094            }
    11195        }
    112     }
    113 
    114     /**
    115      * @param p  WGS84 lat/lon (ellipsoid GRS80) (in degree)
    116      * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
    117      * @throws IOException
    118      */
    119     public EastNorth latlon2eastNorth(LatLon p) {
    120         // translate ellipsoid GRS80 (WGS83) => Clark
    121         LatLon geo = WGS84_to_NTF(p);
    122         double lt = Math.toRadians(geo.lat()); // in radian
    123         double lg = Math.toRadians(geo.lon());
    124 
    125         // check if longitude and latitude are inside the French Lambert zones
    126         if (lt >= cMinLatZone1Radian && lt <= cMaxLatZone1Radian && lg >= cMinLonZonesRadian && lg <= cMaxLonZonesRadian)
    127             return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
    128         return ConicProjection(lt, lg, Xs[0], Ys[0], c[0], n[0]);
    129     }
    130 
    131     public LatLon eastNorth2latlon(EastNorth p) {
    132         LatLon geo;
    133         geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
    134         // translate geodetic system from NTF (ellipsoid Clark) to RGF93/WGS84 (ellipsoid GRS80)
    135         return NTF_to_WGS84(geo);
    136     }
    137 
    138     @Override public String toString() {
     96        updateParameters(DEFAULT_ZONE);
     97    }
     98   
     99    private void updateParameters(int layoutZone) {
     100        this.layoutZone = layoutZone;
     101        ellps = Ellipsoid.clarke;
     102        datum = null; // no datum needed, we have a shift file
     103        nadgrids = ntf_rgf93Grid;
     104        x_0 = x_0s[layoutZone];
     105        lon_0 = 2.0 + 20.0 / 60 + 14.025 / 3600; // 0 grade Paris
     106        if (proj == null) {
     107            proj = new LambertConformalConic();
     108        }
     109        ((LambertConformalConic)proj).updateParametersDirect(
     110                Ellipsoid.clarke, n[layoutZone], c[layoutZone] / ellps.a, y_fs[layoutZone] / ellps.a);
     111    }
     112
     113    @Override
     114    public String toString() {
    139115        return tr("Lambert 4 Zones (France)");
    140116    }
    141117
    142     public String toCode() {
    143         return "EPSG:"+(27561+layoutZone);
     118    @Override
     119    public Integer getEpsgCode() {
     120        return 27561+layoutZone;
    144121    }
    145122
     
    149126    }
    150127
     128    @Override
    151129    public String getCacheDirectoryName() {
    152130        return "lambert";
    153131    }
    154132
    155     /**
    156      * Initializes from geographic coordinates. Note that reference ellipsoid
    157      * used by Lambert is the Clark ellipsoid.
    158      *
    159      * @param lat latitude in grad
    160      * @param lon longitude in grad
    161      * @param Xs  false east (coordinate system origin) in meters
    162      * @param Ys  false north (coordinate system origin) in meters
    163      * @param c   projection constant
    164      * @param n   projection exponent
    165      * @return EastNorth projected coordinates in meter
    166      */
    167     private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) {
    168         double eslt = Ellipsoid.clarke.e * Math.sin(lat);
    169         double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0))
    170                 * Math.pow((1.0 - eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0));
    171         double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0));
    172         double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0));
    173         return new EastNorth(east, north);
    174     }
    175 
    176     /**
    177      * Initializes from projected coordinates (conic projection). Note that
    178      * reference ellipsoid used by Lambert is Clark
    179      *
    180      * @param eastNorth projected coordinates pair in meters
    181      * @param Xs        false east (coordinate system origin) in meters
    182      * @param Ys        false north (coordinate system origin) in meters
    183      * @param c         projection constant
    184      * @param n         projection exponent
    185      * @return LatLon in degree
    186      */
    187     private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) {
    188         double dx = eastNorth.east() - Xs;
    189         double dy = Ys - eastNorth.north();
    190         double R = Math.sqrt(dx * dx + dy * dy);
    191         double gamma = Math.atan(dx / dy);
    192         double l = -1.0 / n * Math.log(Math.abs(R / c));
    193         l = Math.exp(l);
    194         double lon = lg0 + gamma / n;
    195         double lat = 2.0 * Math.atan(l) - Math.PI / 2.0;
    196         double delta = 1.0;
    197         while (delta > epsilon) {
    198             double eslt = Ellipsoid.clarke.e * Math.sin(lat);
    199             double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0 - eslt), Ellipsoid.clarke.e / 2.0) * l) - Math.PI
    200             / 2.0;
    201             delta = Math.abs(nlt - lat);
    202             lat = nlt;
    203         }
    204         return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon)); // in radian
    205     }
    206 
    207     /**
    208      * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert
    209      * geographic, (ellipsoid Clark)
    210      * @throws IOException
    211      */
    212     private LatLon WGS84_to_NTF(LatLon wgs) {
    213         NTV2GridShift gs = new NTV2GridShift(wgs);
    214         if (ntf_rgf93Grid != null) {
    215             ntf_rgf93Grid.gridShiftReverse(gs);
    216             return new LatLon(wgs.lat()+gs.getLatShiftDegrees(), wgs.lon()+gs.getLonShiftPositiveEastDegrees());
    217         }
    218         return new LatLon(0,0);
    219     }
    220 
    221     private LatLon NTF_to_WGS84(LatLon ntf) {
    222         NTV2GridShift gs = new NTV2GridShift(ntf);
    223         if (ntf_rgf93Grid != null) {
    224             ntf_rgf93Grid.gridShiftForward(gs);
    225             return new LatLon(ntf.lat()+gs.getLatShiftDegrees(), ntf.lon()+gs.getLonShiftPositiveEastDegrees());
    226         }
    227         return new LatLon(0,0);
    228     }
    229 
     133    @Override
    230134    public Bounds getWorldBoundsLatLon()
    231135    {
     
    234138                new LatLon(Math.min(zoneLimitsDegree[layoutZone][0] + cMaxOverlappingZonesDegree, Math.toDegrees(cMaxLatZone1Radian)), Math.toDegrees(cMaxLonZonesRadian)));
    235139        return b;
    236     }
    237 
    238     /**
    239      * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
    240      */
    241     public double getDefaultZoomInPPD() {
    242         // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
    243         return 10.0;
    244140    }
    245141
     
    273169    }
    274170
     171    @Override
    275172    public Collection<String> getPreferences(JPanel p) {
    276173        Object prefcb = p.getComponent(2);
     
    281178    }
    282179
     180    @Override
    283181    public void setPreferences(Collection<String> args) {
    284         layoutZone = DEFAULT_ZONE;
     182        int layoutZone = DEFAULT_ZONE;
    285183        if (args != null) {
    286184            try {
     
    295193            } catch(NumberFormatException e) {}
    296194        }
     195        updateParameters(layoutZone);
    297196    }
    298197
     
    306205    }
    307206
     207    @Override
    308208    public Collection<String> getPreferencesFromCode(String code) {
    309209        if (code.startsWith("EPSG:2756") && code.length() == 10) {
  • trunk/src/org/openstreetmap/josm/data/projection/LambertCC9Zones.java

    r4225 r4285  
    1414
    1515import org.openstreetmap.josm.data.Bounds;
    16 import org.openstreetmap.josm.data.coor.EastNorth;
    1716import org.openstreetmap.josm.data.coor.LatLon;
     17import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
     18import org.openstreetmap.josm.data.projection.proj.LambertConformalConic;
    1819import org.openstreetmap.josm.tools.GBC;
    1920import org.openstreetmap.josm.tools.ImageProvider;
    2021
    2122/**
    22  * This class implements the Lambert Conic Conform 9 Zones projection as specified by the IGN
     23 * Lambert Conic Conform 9 Zones projection as specified by the IGN
    2324 * in this document http://professionnels.ign.fr/DISPLAY/000/526/700/5267002/transformation.pdf
    2425 * @author Pieren
    2526 *
    2627 */
    27 public class LambertCC9Zones implements Projection, ProjectionSubPrefs {
    28 
    29     /**
    30      * Lambert 9 zones projection exponents
    31      */
    32     public static final double n[] = { 0.6691500006885269, 0.682018118346418, 0.6946784863203991, 0.7071272481559119,
    33         0.7193606118567315, 0.7313748510399917, 0.7431663060711892, 0.7547313851789208, 0.7660665655489937};
    34 
    35     /**
    36      * Lambert 9 zones projection constants
    37      */
    38     public static final double c[] = { 1.215363305807804E7, 1.2050261119223533E7, 1.195716926884592E7, 1.18737533925172E7,
    39         1.1799460698022118E7, 1.17337838820243E7, 1.16762559948139E7, 1.1626445901183508E7, 1.1583954251630554E7};
    40 
    41     /**
    42      * Lambert 9 zones false east
    43      */
    44     public static final double Xs = 1700000;
    45 
    46     /**
    47      * Lambert 9 zones false north
    48      */
    49     public static final double Ys[] = { 8293467.503439436, 9049604.665107645, 9814691.693461388, 1.0588107871787189E7,
    50         1.1369285637569271E7, 1.2157704903382052E7, 1.2952888086405803E7, 1.3754395745267643E7, 1.4561822739114787E7};
    51 
    52     /**
    53      * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian
    54      */
    55     public static final double lg0 = 0.04079234433198; // 2deg20'14.025"
    56 
    57     /**
    58      * precision in iterative schema
    59      */
    60 
    61     public static final double epsilon = 1e-12;
     28public class LambertCC9Zones extends AbstractProjection implements ProjectionSubPrefs {
    6229
    6330    /**
     
    6734
    6835    public static final double cMinLatZonesDegree = 41.0;
    69     public static final double cMinLatZonesRadian = Math.toRadians(cMinLatZonesDegree);
    70 
    71     public static final double cMinLonZonesRadian = Math.toRadians(-5.0);
    72 
    73     public static final double cMaxLonZonesRadian = Math.toRadians(10.2);
    74 
    75     public static final double lambda0 = Math.toRadians(3);
    76     public static final double e = Ellipsoid.GRS80.e; // but in doc=0.08181919112
    77     public static final double e2 =Ellipsoid.GRS80.e2;
    78     public static final double a = Ellipsoid.GRS80.a;
    7936
    8037    public static final double cMaxOverlappingZones = 1.5;
     
    8239    public static final int DEFAULT_ZONE = 0;
    8340
    84     private static int layoutZone = DEFAULT_ZONE;
     41    private int layoutZone = DEFAULT_ZONE;
    8542
    86     private double L(double phi, double e) {
    87         double sinphi = Math.sin(phi);
    88         return (0.5*Math.log((1+sinphi)/(1-sinphi))) - e/2*Math.log((1+e*sinphi)/(1-e*sinphi));
     43    public LambertCC9Zones() {
     44        this(DEFAULT_ZONE);
    8945    }
    9046
    91     /**
    92      * @param p  WGS84 lat/lon (ellipsoid GRS80) (in degree)
    93      * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
    94      */
    95     public EastNorth latlon2eastNorth(LatLon p) {
    96         double lt = Math.toRadians(p.lat());
    97         double lg = Math.toRadians(p.lon());
    98         if (lt >= cMinLatZonesRadian && lt <= cMaxLatZonesRadian && lg >= cMinLonZonesRadian && lg <= cMaxLonZonesRadian)
    99             return ConicProjection(lt, lg, layoutZone);
    100         return ConicProjection(lt, lg, 0);
     47    public LambertCC9Zones(int layoutZone) {
     48        updateParameters(layoutZone);
    10149    }
    10250
    103     /**
    104      *
    105      * @param lat latitude in grad
    106      * @param lon longitude in grad
    107      * @param nz Lambert CC zone number (from 1 to 9) - 1 !
    108      * @return EastNorth projected coordinates in meter
    109      */
    110     private EastNorth ConicProjection(double lat, double lon, int nz) {
    111         double R = c[nz]*Math.exp(-n[nz]*L(lat,e));
    112         double gamma = n[nz]*(lon-lambda0);
    113         double X = Xs + R*Math.sin(gamma);
    114         double Y = Ys[nz] + -R*Math.cos(gamma);
    115         return new EastNorth(X, Y);
     51    public void updateParameters(int layoutZone) {
     52        ellps = Ellipsoid.GRS80;
     53        datum = GRS80Datum.INSTANCE;
     54        this.layoutZone = layoutZone;
     55        x_0 = 1700000;
     56        y_0 = (layoutZone+1) * 1000000 + 200000;
     57        lon_0 = 3;
     58        double lat_0 = 42 + layoutZone;
     59        double lat_1 = 41.25 + layoutZone;
     60        double lat_2 = 42.75 + layoutZone;
     61        if (proj == null) {
     62            proj = new LambertConformalConic();
     63        }
     64        ((LambertConformalConic)proj).updateParameters2SP(ellps, lat_0, lat_1, lat_2);
    11665    }
    11766
    118     public LatLon eastNorth2latlon(EastNorth p) {
    119         return Geographic(p, layoutZone);
    120     }
    121 
    122     private LatLon Geographic(EastNorth ea, int nz) {
    123         double R = Math.sqrt(Math.pow(ea.getX()-Xs,2)+Math.pow(ea.getY()-Ys[nz], 2));
    124         double gamma = Math.atan((ea.getX()-Xs)/(Ys[nz]-ea.getY()));
    125         double lon = lambda0+gamma/n[nz];
    126         double latIso = (-1/n[nz])*Math.log(Math.abs(R/c[nz]));
    127         double lat = Ellipsoid.GRS80.latitude(latIso, e, epsilon);
    128         return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon));
    129     }
    130 
    131     @Override public String toString() {
     67    @Override
     68    public String toString() {
    13269        return tr("Lambert CC9 Zone (France)");
    13370    }
     
    14077    }
    14178
    142     public String toCode() {
    143         return "EPSG:"+(3942+layoutZone); //CC42 is EPSG:3942 (up to EPSG:3950 for CC50)
     79    @Override
     80    public Integer getEpsgCode() {
     81        return 3942+layoutZone; //CC42 is EPSG:3942 (up to EPSG:3950 for CC50)
    14482    }
    14583
     
    14987    }
    15088
     89    @Override
    15190    public String getCacheDirectoryName() {
    15291        return "lambert";
    15392    }
    15493
    155     /**
    156      * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
    157      */
    158     public double getDefaultZoomInPPD() {
    159         // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
    160         return 10.0;
    161     }
    162 
     94    @Override
    16395    public Bounds getWorldBoundsLatLon()
    16496    {
     
    203135    }
    204136
     137    @Override
    205138    public Collection<String> getPreferences(JPanel p) {
    206139        Object prefcb = p.getComponent(2);
    207140        if(!(prefcb instanceof JComboBox))
    208141            return null;
    209         layoutZone = ((JComboBox)prefcb).getSelectedIndex();
     142        int layoutZone = ((JComboBox)prefcb).getSelectedIndex();
    210143        return Collections.singleton(Integer.toString(layoutZone+1));
    211144    }
    212145
     146    @Override
    213147    public void setPreferences(Collection<String> args)
    214148    {
    215         layoutZone = DEFAULT_ZONE;
     149        int layoutZone = DEFAULT_ZONE;
    216150        if (args != null) {
    217151            try {
     
    226160            } catch(NumberFormatException e) {}
    227161        }
     162        updateParameters(layoutZone);
    228163    }
    229164
     
    237172    }
    238173
     174    @Override
    239175    public Collection<String> getPreferencesFromCode(String code)
    240176    {
     
    246182                if(zoneval >= 0 && zoneval <= 8)
    247183                    return Collections.singleton(String.valueOf(zoneval+1));
    248             } catch(NumberFormatException e) {}
     184            } catch(NumberFormatException ex) {}
    249185        }
    250186        return null;
  • trunk/src/org/openstreetmap/josm/data/projection/LambertEST.java

    r3676 r4285  
    1 //License: GPL. For details, see LICENSE file.
    2 //Thanks to Johan Montagnat and its geoconv java converter application
    3 //(http://www.i3s.unice.fr/~johan/gps/ , published under GPL license)
    4 //from which some code and constants have been reused here.
     1// License: GPL. For details, see LICENSE file.
    52package org.openstreetmap.josm.data.projection;
    63
     
    85
    96import org.openstreetmap.josm.data.Bounds;
    10 import org.openstreetmap.josm.data.coor.EastNorth;
    117import org.openstreetmap.josm.data.coor.LatLon;
     8import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
     9import org.openstreetmap.josm.data.projection.proj.LambertConformalConic;
    1210
    13 public class LambertEST implements Projection {
     11/**
     12 * Estonian Coordinate System of 1997.
     13 *
     14 * Thanks to Johan Montagnat and its geoconv java converter application
     15 * (http://www.i3s.unice.fr/~johan/gps/ , published under GPL license)
     16 * from which some code and constants have been reused here.
     17 */
     18public class LambertEST extends AbstractProjection {
    1419
    15     public static final double ef = 500000; //Easting at false origin            = 5000000 m
    16     public static final double nf = 6375000; //Northing at false origin           = 6375000 m
    17     public static final double lat1 = Math.toRadians(59 + 1.0/3.0); //Latitude of 1st standard parallel  = 59o20`0`` N
    18     public static final double lat2 = Math.toRadians( 58);//Latitude of 2nd standard parallel  = 58o0`0`` N
    19     public static final double latf = Math.toRadians(57.517553930555555555555555555556);//'Latitude of false origin = 57o31`3,19415`` N
    20     public static final double lonf = Math.toRadians( 24.0);
    21     public static final double a = 6378137;
    22     public static final double ee = 0.081819191042815792;
    23     public static final double m1 = Math.cos(lat1) / (Math.sqrt(1 - ee *ee * Math.pow(Math.sin(lat1), 2)));
    24     public static final double m2 = Math.cos(lat2) / (Math.sqrt(1 - ee *ee * Math.pow(Math.sin(lat2), 2)));
    25     public static final double t1 = Math.tan(Math.PI / 4.0 - lat1 / 2.0)
    26     / Math.pow(( (1.0 - ee * Math.sin(lat1)) / (1.0 + ee * Math.sin(lat1))) ,(ee / 2.0));
    27     public static final double t2 = Math.tan(Math.PI / 4.0 - lat2 / 2.0)
    28     / Math.pow(( (1.0 - ee * Math.sin(lat2)) / (1.0 + ee * Math.sin(lat2))) ,(ee / 2.0));
    29     public static final double tf = Math.tan(Math.PI / 4.0 - latf / 2.0)
    30     / Math.pow(( (1.0 - ee * Math.sin(latf)) / (1.0 + ee * Math.sin(latf))) ,(ee / 2.0));
    31     public static final double n  = (Math.log(m1) - Math.log(m2))
    32     / (Math.log(t1) - Math.log(t2));
    33     public static final double f  = m1 / (n * Math.pow(t1, n));
    34     public static final double rf  = a * f * Math.pow(tf, n);
    35 
    36     /**
    37      * precision in iterative schema
    38      */
    39     public static final double epsilon = 1e-11;
    40 
    41     /**
    42      * @param p  WGS84 lat/lon (ellipsoid GRS80) (in degree)
    43      * @return eastnorth projection in Lambert Zone (ellipsoid GRS80)
    44      */
    45     public EastNorth latlon2eastNorth(LatLon p)
    46     {
    47 
    48         double t = Math.tan(Math.PI / 4.0 - Math.toRadians(p.lat()) / 2.0)
    49         / Math.pow(( (1.0 - ee * Math.sin(Math.toRadians(p.lat()))) / (1.0
    50                 + ee * Math.sin(Math.toRadians(p.lat())))) ,(ee / 2.0));
    51         double r = a * f * Math.pow(t, n);
    52         double theta = n * (Math.toRadians(p.lon()) - lonf);
    53 
    54         double x = ef + r * Math.sin(theta);     //587446.7
    55         double y = nf + rf - r * Math.cos(theta); //6485401.6
    56 
    57         return new EastNorth(x,y);
    58     }
    59 
    60     public static double iterateAngle(double e, double t)
    61     {
    62         double a1 = 0.0;
    63         double a2 = 3.1415926535897931;
    64         double a = 1.5707963267948966;
    65         double b = 1.5707963267948966 - (2.0 * Math.atan(t * Math.pow((1.0
    66                 - (e * Math.sin(a))) / (1.0 + (e * Math.sin(a))), e / 2.0)));
    67         while (Math.abs(a-b) > epsilon)
    68         {
    69             a = a1 + ((a2 - a1) / 2.0);
    70             b = 1.5707963267948966 - (2.0 * Math.atan(t * Math.pow((1.0
    71                     - (e * Math.sin(a))) / (1.0 + (e * Math.sin(a))), e / 2.0)));
    72             if (a1 == a2)
    73                 return 0.0;
    74             if (b > a) {
    75                 a1 = a;
    76             } else {
    77                 a2 = a;
    78             }
    79         }
    80         return b;
    81     }
    82 
    83     public LatLon eastNorth2latlon(EastNorth p)
    84     {
    85         double r = Math.sqrt(Math.pow((p.getX() - ef), 2.0) + Math.pow((rf
    86                 - p.getY() + nf), 2.0) ) * Math.signum(n);
    87         double T = Math.pow((r / (a * f)), (1.0/ n)) ;
    88         double theta = Math.atan((p.getX() - ef) / (rf - p.getY() + nf));
    89         double y = (theta / n + lonf) ;
    90         double x = iterateAngle(ee, T);
    91         return new LatLon(Math.toDegrees(x),Math.toDegrees(y));
     20    public LambertEST() {
     21        ellps = Ellipsoid.GRS80;
     22        datum = GRS80Datum.INSTANCE;
     23        lon_0 = 24;
     24        double lat_0 = 57.517553930555555555555555555556;
     25        double lat_1 = 59 + 1.0/3.0;
     26        double lat_2 = 58;
     27        x_0 = 500000;
     28        y_0 = 6375000;
     29        proj = new LambertConformalConic();
     30        ((LambertConformalConic) proj).updateParameters2SP(ellps, lat_0, lat_1, lat_2);
    9231    }
    9332
     
    9736    }
    9837
    99     public String toCode() {
    100         return "EPSG:3301";
     38    @Override
     39    public Integer getEpsgCode() {
     40        return 3301;
    10141    }
    10242
     
    10646    }
    10747
     48    @Override
    10849    public String getCacheDirectoryName() {
    10950        return "lambertest";
    11051    }
    11152
     53    @Override
    11254    public Bounds getWorldBoundsLatLon()
    11355    {
     
    11759    }
    11860
    119     public double getDefaultZoomInPPD() {
    120         return 10;
    121     }
    12261}
  • trunk/src/org/openstreetmap/josm/data/projection/Mercator.java

    r3922 r4285  
    55
    66import org.openstreetmap.josm.data.Bounds;
    7 import org.openstreetmap.josm.data.coor.EastNorth;
    87import org.openstreetmap.josm.data.coor.LatLon;
     8import org.openstreetmap.josm.data.projection.datum.WGS84Datum;
    99
    1010/**
    11  * Implement Mercator Projection code, coded after documentation
    12  * from wikipedia.
     11 * Mercator Projection
    1312 *
    1413 * The center of the mercator projection is always the 0 grad
     
    2019 * @author imi
    2120 */
    22 public class Mercator implements Projection {
     21public class Mercator extends AbstractProjection {
    2322
    24     final double radius = 6378137.0;
    25 
    26     public EastNorth latlon2eastNorth(LatLon p) {
    27         return new EastNorth(
    28                 p.lon()*Math.PI/180*radius,
    29                 Math.log(Math.tan(Math.PI/4+p.lat()*Math.PI/360))*radius);
     23    public Mercator() {
     24        ellps = Ellipsoid.WGS84;
     25        datum = WGS84Datum.INSTANCE;
     26        proj = new org.openstreetmap.josm.data.projection.proj.Mercator();
    3027    }
    3128
    32     public LatLon eastNorth2latlon(EastNorth p) {
    33         return new LatLon(
    34                 Math.atan(Math.sinh(p.north()/radius))*180/Math.PI,
    35                 p.east()/radius*180/Math.PI);
    36     }
    37 
    38     @Override public String toString() {
     29    @Override
     30    public String toString() {
    3931        return tr("Mercator");
    4032    }
    4133
    42     public String toCode() {
    43         return "EPSG:3857"; /* initially they used 3785 but that has been superseded, see http://www.epsg-registry.org/ */
     34    @Override
     35    public Integer getEpsgCode() {
     36        /* initially they used 3785 but that has been superseded,
     37         * see http://www.epsg-registry.org/ */
     38        return 3857;
    4439    }
    4540
     
    4944    }
    5045
     46    @Override
    5147    public String getCacheDirectoryName() {
    5248        return "mercator";
    5349    }
    5450
     51    @Override
    5552    public Bounds getWorldBoundsLatLon()
    5653    {
     
    6057    }
    6158
    62     public double getDefaultZoomInPPD() {
    63         // This will set the scale bar to about 100 km
    64         return 1000.0;/*0.000158*/
    65     }
    6659}
  • trunk/src/org/openstreetmap/josm/data/projection/Puwg.java

    r3779 r4285  
    77import java.awt.GridBagLayout;
    88import java.awt.event.ActionListener;
    9 import java.text.DecimalFormat;
    109import java.util.Collection;
    1110import java.util.Collections;
     
    1615
    1716import org.openstreetmap.josm.data.Bounds;
    18 import org.openstreetmap.josm.data.coor.EastNorth;
    1917import org.openstreetmap.josm.data.coor.LatLon;
     18import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
    2019import org.openstreetmap.josm.tools.GBC;
    2120
     
    2625 * @author steelman
    2726 */
    28 public class Puwg extends UTM implements Projection,ProjectionSubPrefs {
     27public class Puwg extends AbstractProjection implements ProjectionSubPrefs {
     28   
    2929    public static final int DEFAULT_ZONE = 0;
    30     private int zone = DEFAULT_ZONE;
    31 
    32     static PuwgData[] Zones = new PuwgData[]{
     30   
     31    private int zone;
     32
     33    static PuwgData[] Zones = new PuwgData[] {
    3334        new Epsg2180(),
    3435        new Epsg2176(),
     
    3839    };
    3940
    40     private static DecimalFormat decFormatter = new DecimalFormat("###0");
    41 
    42     @Override
    43     public EastNorth latlon2eastNorth(LatLon p) {
     41    public Puwg() {
     42        this(DEFAULT_ZONE);
     43    }
     44
     45    public Puwg(int zone) {
     46        ellps = Ellipsoid.GRS80;
     47        proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
     48        datum = GRS80Datum.INSTANCE;
     49        updateParameters(zone);
     50    }
     51   
     52    public void updateParameters(int zone) {
     53        this.zone = zone;
    4454        PuwgData z = Zones[zone];
    45         double easting = z.getPuwgFalseEasting();
    46         double northing = z.getPuwgFalseNorthing();
    47         double scale = z.getPuwgScaleFactor();
    48         double center = z.getPuwgCentralMeridian(); /* in radians */
    49         EastNorth a = mapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), center);
    50         return new EastNorth(a.east() * scale + easting, a.north() * scale + northing);
    51     }
    52 
    53     @Override
    54     public LatLon eastNorth2latlon(EastNorth p) {
    55         PuwgData z = Zones[zone];
    56         double easting = z.getPuwgFalseEasting();
    57         double northing = z.getPuwgFalseNorthing();
    58         double scale = z.getPuwgScaleFactor();
    59         double center = z.getPuwgCentralMeridian(); /* in radians */
    60         return mapXYToLatLon((p.east() - easting)/scale, (p.north() - northing)/scale, center);
    61     }
    62 
    63     @Override public String toString() {
     55        x_0 = z.getPuwgFalseEasting();
     56        y_0 = z.getPuwgFalseNorthing();
     57        lon_0 = z.getPuwgCentralMeridianDeg();
     58        k_0 = z.getPuwgScaleFactor();
     59    }
     60
     61    @Override
     62    public String toString() {
    6463        return tr("PUWG (Poland)");
    6564    }
    6665
    6766    @Override
    68     public String toCode() {
    69         return Zones[zone].toCode();
     67    public Integer getEpsgCode() {
     68        return Zones[zone].getEpsgCode();
    7069    }
    7170
     
    8382    public Bounds getWorldBoundsLatLon() {
    8483        return Zones[zone].getWorldBoundsLatLon();
    85     }
    86 
    87     @Override
    88     public double getDefaultZoomInPPD() {
    89         // This will set the scale bar to about 100 km
    90         return 0.009;
    91     }
    92 
    93     public String eastToString(EastNorth p) {
    94         return decFormatter.format(p.east());
    95     }
    96 
    97     public String northToString(EastNorth p) {
    98         return decFormatter.format(p.north());
    9984    }
    10085
     
    135120
    136121    @Override
    137     public Collection<String> getPreferencesFromCode(String code)
    138     {
    139         for (PuwgData p : Puwg.Zones)
    140         {
    141             if(code.equals(p.toCode()))
     122    public Collection<String> getPreferencesFromCode(String code) {
     123        for (PuwgData p : Puwg.Zones) {
     124            if (code.equals(p.toCode()))
    142125                return Collections.singleton(code);
    143126        }
     
    146129
    147130    @Override
    148     public void setPreferences(Collection<String> args)
    149     {
    150         zone = DEFAULT_ZONE;
    151         if(args != null)
    152         {
     131    public void setPreferences(Collection<String> args) {
     132        int z = DEFAULT_ZONE;
     133        if (args != null) {
    153134            try {
    154                 for(String s : args)
    155                 {
    156                     for (int i=0; i < Puwg.Zones.length; ++i)
    157                         if(s.equals(Zones[i].toCode())) {
    158                             zone = i;
     135                for (String s : args) {
     136                    for (int i=0; i < Zones.length; ++i)
     137                        if (s.equals(Zones[i].toCode())) {
     138                            z = i;
     139                            break;
    159140                        }
    160141                    break;
     
    162143            } catch (NullPointerException e) {}
    163144        }
     145        updateParameters(z);
    164146    }
    165147}
     
    173155
    174156    // Projection methods
     157    Integer getEpsgCode();
    175158    String toCode();
    176159    String getCacheDirectoryName();
     
    189172    }
    190173
     174    @Override
     175    public Integer getEpsgCode() {
     176        return 2180;
     177    }
     178   
     179    @Override
    191180    public String toCode() {
    192         return "EPSG:2180";
    193     }
    194 
     181        return "EPSG:" + getEpsgCode();
     182    }
     183
     184    @Override
    195185    public String getCacheDirectoryName() {
    196186        return "epsg2180";
    197187    }
    198188
     189    @Override
    199190    public Bounds getWorldBoundsLatLon()
    200191    {
     
    204195    }
    205196
    206     public double getPuwgCentralMeridianDeg() { return Epsg2180CentralMeridian; }
    207     public double getPuwgCentralMeridian() { return Math.toRadians(Epsg2180CentralMeridian); }
    208     public double getPuwgFalseEasting() { return Epsg2180FalseEasting; }
    209     public double getPuwgFalseNorthing() { return Epsg2180FalseNorthing; }
    210     public double getPuwgScaleFactor() { return Epsg2180ScaleFactor; }
     197    @Override public double getPuwgCentralMeridianDeg() { return Epsg2180CentralMeridian; }
     198    @Override public double getPuwgCentralMeridian() { return Math.toRadians(Epsg2180CentralMeridian); }
     199    @Override public double getPuwgFalseEasting() { return Epsg2180FalseEasting; }
     200    @Override public double getPuwgFalseNorthing() { return Epsg2180FalseNorthing; }
     201    @Override public double getPuwgScaleFactor() { return Epsg2180ScaleFactor; }
    211202}
    212203
     
    217208    private static final double PuwgScaleFactor = 0.999923;
    218209    //final private double[] Puwg2000CentralMeridian = {15.0, 18.0, 21.0, 24.0};
    219     final private String[] Puwg2000Code = { "EPSG:2176",  "EPSG:2177", "EPSG:2178", "EPSG:2179"};
    220     final private String[] Puwg2000CDName = { "epsg2176",  "epsg2177", "epsg2178", "epsg2179"};
     210    final private Integer[] Puwg2000Code = { 2176,  2177, 2178, 2179 };
     211    final private String[] Puwg2000CDName = { "epsg2176",  "epsg2177", "epsg2178", "epsg2179" };
    221212
    222213    @Override public String toString() {
     
    224215    }
    225216
     217    @Override
     218    public Integer getEpsgCode() {
     219        return Puwg2000Code[getZoneIndex()];
     220    }
     221
     222    @Override
    226223    public String toCode() {
    227         return Puwg2000Code[getZoneIndex()];
    228     }
    229 
     224        return "EPSG:" + getEpsgCode();
     225    }
     226
     227    @Override
    230228    public String getCacheDirectoryName() {
    231229        return Puwg2000CDName[getZoneIndex()];
    232230    }
    233231
     232    @Override
    234233    public Bounds getWorldBoundsLatLon()
    235234    {
     
    239238    }
    240239
    241     public double getPuwgCentralMeridianDeg() { return getZone() * 3.0; }
    242     public double getPuwgCentralMeridian() { return Math.toRadians(getZone() * 3.0); }
    243     public double getPuwgFalseNorthing() { return PuwgFalseNorthing;}
    244     public double getPuwgFalseEasting() { return 1e6 * getZone() + PuwgFalseEasting; }
    245     public double getPuwgScaleFactor() { return PuwgScaleFactor; }
     240    @Override public double getPuwgCentralMeridianDeg() { return getZone() * 3.0; }
     241    @Override public double getPuwgCentralMeridian() { return Math.toRadians(getZone() * 3.0); }
     242    @Override public double getPuwgFalseNorthing() { return PuwgFalseNorthing;}
     243    @Override public double getPuwgFalseEasting() { return 1e6 * getZone() + PuwgFalseEasting; }
     244    @Override public double getPuwgScaleFactor() { return PuwgScaleFactor; }
    246245    public abstract int getZone();
    247246
  • trunk/src/org/openstreetmap/josm/data/projection/SwissGrid.java

    r3779 r4285  
    1212
    1313import org.openstreetmap.josm.data.Bounds;
    14 import org.openstreetmap.josm.data.coor.EastNorth;
    1514import org.openstreetmap.josm.data.coor.LatLon;
     15import org.openstreetmap.josm.data.projection.datum.ThreeParameterDatum;
     16import org.openstreetmap.josm.data.projection.proj.SwissObliqueMercator;
    1617import org.openstreetmap.josm.gui.widgets.HtmlPanel;
    1718import org.openstreetmap.josm.tools.GBC;
    1819
    1920/**
    20  * Projection for the SwissGrid CH1903 / L03, see http://de.wikipedia.org/wiki/Swiss_Grid.
     21 * SwissGrid CH1903 / L03, see http://de.wikipedia.org/wiki/Swiss_Grid.
     22 *
     23 * Actually, what we have here, is CH1903+ (EPSG:2056), but without
     24 * the additional false easting of 2000km and false northing 1000 km.
    2125 *
    22  * Calculations are based on formula from
    23  * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf
    24  *
    25  * August 2010 update to this formula (rigorous formulas)
    26  * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf
     26 * To get to CH1903, a shift file is required. So currently, there are errors
     27 * up to 1.6m (depending on the location).
     28 * 
     29 * This projection does not have any parameters, it only implements
     30 * ProjectionSubPrefs to show a warning that the grid file correction is not done.
    2731 */
    28 public class SwissGrid implements Projection, ProjectionSubPrefs {
     32public class SwissGrid extends AbstractProjection implements ProjectionSubPrefs {
    2933
    30     private static final double dX = 674.374;
    31     private static final double dY = 15.056;
    32     private static final double dZ = 405.346;
    33 
    34     private static final double phi0 = Math.toRadians(46.0 + 57.0 / 60 + 8.66 / 3600);
    35     private static final double lambda0 = Math.toRadians(7.0 + 26.0 / 60 + 22.50 / 3600);
    36     private static final double R = Ellipsoid.Bessel1841.a * Math.sqrt(1 - Ellipsoid.Bessel1841.e2) / (1 - (Ellipsoid.Bessel1841.e2 * Math.pow(Math.sin(phi0), 2)));
    37     private static final double alpha = Math.sqrt(1 + (Ellipsoid.Bessel1841.eb2 * Math.pow(Math.cos(phi0), 4)));
    38     private static final double b0 = Math.asin(Math.sin(phi0) / alpha);
    39     private static final double K = Math.log(Math.tan(Math.PI / 4 + b0 / 2)) - alpha
    40     * Math.log(Math.tan(Math.PI / 4 + phi0 / 2)) + alpha * Ellipsoid.Bessel1841.e / 2
    41     * Math.log((1 + Ellipsoid.Bessel1841.e * Math.sin(phi0)) / (1 - Ellipsoid.Bessel1841.e * Math.sin(phi0)));
    42 
    43     private static final double xTrans = 200000;
    44     private static final double yTrans = 600000;
    45 
    46     private static final double DELTA_PHI = 1e-11;
    47 
    48     private LatLon correctEllipoideGSR80toBressel1841(LatLon coord) {
    49         double[] XYZ = Ellipsoid.WGS84.latLon2Cart(coord);
    50         XYZ[0] -= dX;
    51         XYZ[1] -= dY;
    52         XYZ[2] -= dZ;
    53         return Ellipsoid.Bessel1841.cart2LatLon(XYZ);
    54     }
    55 
    56     private LatLon correctEllipoideBressel1841toGRS80(LatLon coord) {
    57         double[] XYZ = Ellipsoid.Bessel1841.latLon2Cart(coord);
    58         XYZ[0] += dX;
    59         XYZ[1] += dY;
    60         XYZ[2] += dZ;
    61         return Ellipsoid.WGS84.cart2LatLon(XYZ);
    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     @Override
    69     public EastNorth latlon2eastNorth(LatLon wgs) {
    70         LatLon coord = correctEllipoideGSR80toBressel1841(wgs);
    71         double phi = Math.toRadians(coord.lat());
    72         double lambda = Math.toRadians(coord.lon());
    73 
    74         double S = alpha * Math.log(Math.tan(Math.PI / 4 + phi / 2)) - alpha * Ellipsoid.Bessel1841.e / 2
    75         * Math.log((1 + Ellipsoid.Bessel1841.e * Math.sin(phi)) / (1 - Ellipsoid.Bessel1841.e * Math.sin(phi))) + K;
    76         double b = 2 * (Math.atan(Math.exp(S)) - Math.PI / 4);
    77         double l = alpha * (lambda - lambda0);
    78 
    79         double lb = Math.atan2(Math.sin(l), Math.sin(b0) * Math.tan(b) + Math.cos(b0) * Math.cos(l));
    80         double bb = Math.asin(Math.cos(b0) * Math.sin(b) - Math.sin(b0) * Math.cos(b) * Math.cos(l));
    81 
    82         double y = R * lb;
    83         double x = R / 2 * Math.log((1 + Math.sin(bb)) / (1 - Math.sin(bb)));
    84 
    85         return new EastNorth(y + yTrans, x + xTrans);
    86     }
    87 
    88     /**
    89      * @param xy SwissGrid east/north (in meters)
    90      * @return LatLon WGS84 (in degree)
    91      */
    92     @Override
    93     public LatLon eastNorth2latlon(EastNorth xy) {
    94         double x = xy.north() - xTrans;
    95         double y = xy.east() - yTrans;
    96 
    97         double lb = y / R;
    98         double bb = 2 * (Math.atan(Math.exp(x / R)) - Math.PI / 4);
    99 
    100         double b = Math.asin(Math.cos(b0) * Math.sin(bb) + Math.sin(b0) * Math.cos(bb) * Math.cos(lb));
    101         double l = Math.atan2(Math.sin(lb), Math.cos(b0) * Math.cos(lb) - Math.sin(b0) * Math.tan(bb));
    102 
    103         double lambda = lambda0 + l / alpha;
    104         double phi = b;
    105         double S = 0;
    106 
    107         double prevPhi = -1000;
    108         int iteration = 0;
    109         // iteration to finds S and phi
    110         while (Math.abs(phi - prevPhi) > DELTA_PHI) {
    111             if (++iteration > 30)
    112                 throw new RuntimeException("Two many iterations");
    113             prevPhi = phi;
    114             S = 1 / alpha * (Math.log(Math.tan(Math.PI / 4 + b / 2)) - K) + Ellipsoid.Bessel1841.e
    115             * Math.log(Math.tan(Math.PI / 4 + Math.asin(Ellipsoid.Bessel1841.e * Math.sin(phi)) / 2));
    116             phi = 2 * Math.atan(Math.exp(S)) - Math.PI / 2;
    117         }
    118 
    119         LatLon coord = correctEllipoideBressel1841toGRS80(new LatLon(Math.toDegrees(phi), Math.toDegrees(lambda)));
    120         return coord;
     34    public SwissGrid() {
     35        ellps = Ellipsoid.Bessel1841;
     36        datum = new ThreeParameterDatum("SwissGrid Datum", null, ellps, 674.374, 15.056, 405.346);
     37        x_0 = 600000;
     38        y_0 = 200000;
     39        lon_0 = 7.0 + 26.0/60 + 22.50/3600;
     40        double lat_0 = 46.0 + 57.0/60 + 8.66/3600;
     41        proj = new SwissObliqueMercator(ellps, lat_0);
    12142    }
    12243
     
    12748
    12849    @Override
    129     public String toCode() {
    130         return "EPSG:21781";
     50    public Integer getEpsgCode() {
     51        return 21781;
    13152    }
    13253
     
    14465    public Bounds getWorldBoundsLatLon() {
    14566        return new Bounds(new LatLon(45.7, 5.7), new LatLon(47.9, 10.6));
    146     }
    147 
    148     @Override
    149     public double getDefaultZoomInPPD() {
    150         // This will set the scale bar to about 100 m
    151         return 1.01;
    15267    }
    15368
  • trunk/src/org/openstreetmap/josm/data/projection/TransverseMercatorLV.java

    r3635 r4285  
    66import org.openstreetmap.josm.data.Bounds;
    77import org.openstreetmap.josm.data.coor.LatLon;
     8import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
    89
    910/**
     
    1314 * @author Viesturs Zarins
    1415 */
    15 public class TransverseMercatorLV extends TransverseMercator {
     16public class TransverseMercatorLV extends AbstractProjection {
    1617
    17     public TransverseMercatorLV()
    18     {
    19         setProjectionParameters(24, 500000, -6000000);
     18    public TransverseMercatorLV() {
     19        ellps = Ellipsoid.GRS80;
     20        proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
     21        datum = GRS80Datum.INSTANCE;
     22        lon_0 = 24;
     23        x_0 = 500000;
     24        y_0 = -6000000;
     25        k_0 = 0.9996;
    2026    }
    21 
    22     @Override public String toString() {
     27   
     28    @Override
     29    public String toString() {
    2330        return tr("LKS-92 (Latvia TM)");
    2431    }
    2532
    26     private int epsgCode() {
     33    @Override
     34    public Integer getEpsgCode() {
    2735        return 3059;
    28     }
    29 
    30     @Override
    31     public String toCode() {
    32         return "EPSG:"+ epsgCode();
    3336    }
    3437
     
    3841    }
    3942
     43    @Override
    4044    public String getCacheDirectoryName() {
    41         return "epsg"+ epsgCode();
     45        return "epsg"+ getEpsgCode();
    4246    }
    4347
  • trunk/src/org/openstreetmap/josm/data/projection/UTM.java

    r4275 r4285  
    1919import org.openstreetmap.josm.data.Bounds;
    2020import org.openstreetmap.josm.data.coor.LatLon;
     21import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
    2122import org.openstreetmap.josm.tools.GBC;
    2223
     
    2526 * @author Dirk Stöcker
    2627 * code based on JavaScript from Chuck Taylor
     28 *
    2729 */
    28 public class UTM extends TransverseMercator implements ProjectionSubPrefs {
     30public class UTM extends AbstractProjection implements ProjectionSubPrefs {
    2931
    3032    private static final int DEFAULT_ZONE = 30;
    31     private int zone = DEFAULT_ZONE;
    32 
    33     public enum Hemisphere {North, South}
     33    private int zone;
     34
     35    public enum Hemisphere { North, South }
    3436    private static final Hemisphere DEFAULT_HEMISPHERE = Hemisphere.North;
    35     private Hemisphere hemisphere = DEFAULT_HEMISPHERE;
    36 
    37     private boolean offset = false;
    38 
    39     public UTM()
    40     {
    41         updateParameters();
     37    private Hemisphere hemisphere;
     38
     39    /**
     40     * Applies an additional false easting of 3000000 m if true.
     41     */
     42    private boolean offset;
     43
     44    public UTM() {
     45        this(DEFAULT_ZONE, DEFAULT_HEMISPHERE, false);
     46    }
     47
     48    public UTM(int zone, Hemisphere hemisphere, boolean offset) {
     49        ellps = Ellipsoid.GRS80;
     50        proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
     51        datum = GRS80Datum.INSTANCE;
     52        updateParameters(zone, hemisphere, offset);
     53    }
     54
     55    public void updateParameters(int zone, Hemisphere hemisphere, boolean offset) {
     56        this.zone = zone;
     57        this.hemisphere = hemisphere;
     58        this.offset = offset;
     59        x_0 = 500000 + (offset ? 3000000 : 0);
     60        y_0 = hemisphere == Hemisphere.North ? 0 : 10000000;
     61        lon_0 = getUtmCentralMeridianDeg(zone);
     62        k_0 = 0.9996;
    4263    }
    4364
     
    5677     *
    5778     */
    58     private double UTMCentralMeridianDeg(int zone)
     79    private double getUtmCentralMeridianDeg(int zone)
    5980    {
    6081        return -183.0 + (zone * 6.0);
    6182    }
    6283
    63     @Override public String toString() {
     84    public int getzone() {
     85        return zone;
     86    }
     87
     88    @Override
     89    public String toString() {
    6490        return tr("UTM");
    6591    }
    6692
    67     private void updateParameters() {
    68         setProjectionParameters(this.UTMCentralMeridianDeg(getzone()), getEastOffset(), getNorthOffset());
    69     }
    70 
    71     public int getzone()
    72     {
    73         return zone;
    74     }
    75 
    76     private double getEastOffset() {
    77         return 500000 + (offset?3000000:0);
    78     }
    79 
    80     private double getNorthOffset() {
    81         if (hemisphere == Hemisphere.North)
    82             return 0;
    83         else
    84             return 10000000;
    85     }
    86 
    87     private int epsgCode() {
     93    @Override
     94    public Integer getEpsgCode() {
    8895        return ((offset?325800:32600) + getzone() + (hemisphere == Hemisphere.South?100:0));
    89     }
    90 
    91     public String toCode() {
    92         return "EPSG:"+ epsgCode();
    9396    }
    9497
     
    98101    }
    99102
     103    @Override
    100104    public String getCacheDirectoryName() {
    101         return "epsg"+ epsgCode();
    102     }
    103 
     105        return "epsg"+ getEpsgCode();
     106    }
     107
     108    @Override
    104109    public Bounds getWorldBoundsLatLon()
    105110    {
    106111        if (hemisphere == Hemisphere.North)
    107112            return new Bounds(
    108                     new LatLon(-5.0, UTMCentralMeridianDeg(getzone())-5.0),
    109                     new LatLon(85.0, UTMCentralMeridianDeg(getzone())+5.0));
     113                    new LatLon(-5.0, getUtmCentralMeridianDeg(getzone())-5.0),
     114                    new LatLon(85.0, getUtmCentralMeridianDeg(getzone())+5.0));
    110115        else
    111116            return new Bounds(
    112                     new LatLon(-85.0, UTMCentralMeridianDeg(getzone())-5.0),
    113                     new LatLon(5.0, UTMCentralMeridianDeg(getzone())+5.0));
     117                    new LatLon(-85.0, getUtmCentralMeridianDeg(getzone())-5.0),
     118                    new LatLon(5.0, getUtmCentralMeridianDeg(getzone())+5.0));
    114119    }
    115120
     
    173178    }
    174179
     180    @Override
    175181    public Collection<String> getPreferences(JPanel p) {
    176182        int zone = DEFAULT_ZONE;
     
    199205    }
    200206
    201     public void setPreferences(Collection<String> args)
    202     {
    203         zone = DEFAULT_ZONE;
    204         hemisphere = DEFAULT_HEMISPHERE;
    205         offset = false;
     207    @Override
     208    public void setPreferences(Collection<String> args) {
     209        int zone = DEFAULT_ZONE;
     210        Hemisphere hemisphere = DEFAULT_HEMISPHERE;
     211        boolean offset = false;
    206212
    207213        if(args != null)
     
    223229            }
    224230        }
    225         updateParameters();
    226     }
    227 
     231        updateParameters(zone, hemisphere, offset);
     232    }
     233
     234    @Override
    228235    public String[] allCodes() {
    229236        ArrayList<String> projections = new ArrayList<String>(60*4);
     
    236243        }
    237244        return projections.toArray(new String[0]);
    238 
    239     }
    240 
    241     public Collection<String> getPreferencesFromCode(String code)
    242     {
     245    }
     246
     247    @Override
     248    public Collection<String> getPreferencesFromCode(String code) {
     249
    243250        boolean offset = code.startsWith("EPSG:3258") || code.startsWith("EPSG:3259");
    244251
  • trunk/src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java

    r3779 r4285  
    22package org.openstreetmap.josm.data.projection;
    33
    4 /**
    5  * This class implements all projections for French departements in the Caribbean Sea and
    6  * Indian Ocean using the UTM transvers Mercator projection and specific geodesic settings (7 parameters transformation algorithm).
    7  */
    84import static org.openstreetmap.josm.tools.I18n.tr;
    95
     
    1814
    1915import org.openstreetmap.josm.data.Bounds;
    20 import org.openstreetmap.josm.data.coor.EastNorth;
    2116import org.openstreetmap.josm.data.coor.LatLon;
     17import org.openstreetmap.josm.data.projection.datum.Datum;
     18import org.openstreetmap.josm.data.projection.datum.GRS80Datum;
     19import org.openstreetmap.josm.data.projection.datum.SevenParameterDatum;
     20import org.openstreetmap.josm.data.projection.datum.ThreeParameterDatum;
    2221import org.openstreetmap.josm.tools.GBC;
    2322
    24 public class UTM_France_DOM implements Projection, ProjectionSubPrefs {
     23/**
     24 * This class implements all projections for French departements in the Caribbean Sea and
     25 * Indian Ocean using the UTM transvers Mercator projection and specific geodesic settings (7 parameters transformation algorithm).
     26 *
     27 */
     28public class UTM_France_DOM extends AbstractProjection implements ProjectionSubPrefs {
    2529
    26     private static String FortMarigotName = tr("Guadeloupe Fort-Marigot 1949");
    27     private static String SainteAnneName = tr("Guadeloupe Ste-Anne 1948");
    28     private static String MartiniqueName = tr("Martinique Fort Desaix 1952");
    29     private static String Reunion92Name = tr("Reunion RGR92");
    30     private static String Guyane92Name = tr("Guyane RGFG95");
    31     public static String[] utmGeodesicsNames = { FortMarigotName, SainteAnneName, MartiniqueName, Reunion92Name, Guyane92Name};
     30    private final static String FortMarigotName = tr("Guadeloupe Fort-Marigot 1949");
     31    private final static String SainteAnneName = tr("Guadeloupe Ste-Anne 1948");
     32    private final static String MartiniqueName = tr("Martinique Fort Desaix 1952");
     33    private final static String Reunion92Name = tr("Reunion RGR92");
     34    private final static String Guyane92Name = tr("Guyane RGFG95");
     35    private final static String[] utmGeodesicsNames = { FortMarigotName, SainteAnneName, MartiniqueName, Reunion92Name, Guyane92Name};
    3236
    33     private Bounds FortMarigotBounds = new Bounds( new LatLon(17.6,-63.25), new LatLon(18.5,-62.5));
    34     private Bounds SainteAnneBounds = new Bounds( new LatLon(15.8,-61.9), new LatLon(16.6,-60.9));
    35     private Bounds MartiniqueBounds = new Bounds( new LatLon(14.25,-61.25), new LatLon(15.025,-60.725));
    36     private Bounds ReunionBounds = new Bounds( new LatLon(-25.92,37.58), new LatLon(-10.6, 58.27));
    37     private Bounds GuyaneBounds = new Bounds( new LatLon(2.16 , -54.0), new LatLon(9.06 , -49.62));
    38     private Bounds[] utmBounds = { FortMarigotBounds, SainteAnneBounds, MartiniqueBounds, ReunionBounds, GuyaneBounds};
     37    private final static Bounds FortMarigotBounds = new Bounds( new LatLon(17.6,-63.25), new LatLon(18.5,-62.5));
     38    private final static Bounds SainteAnneBounds = new Bounds( new LatLon(15.8,-61.9), new LatLon(16.6,-60.9));
     39    private final static Bounds MartiniqueBounds = new Bounds( new LatLon(14.25,-61.25), new LatLon(15.025,-60.725));
     40    private final static Bounds ReunionBounds = new Bounds( new LatLon(-25.92,37.58), new LatLon(-10.6, 58.27));
     41    private final static Bounds GuyaneBounds = new Bounds( new LatLon(2.16 , -54.0), new LatLon(9.06 , -49.62));
     42    private final static Bounds[] utmBounds = { FortMarigotBounds, SainteAnneBounds, MartiniqueBounds, ReunionBounds, GuyaneBounds };
    3943
    40     private String FortMarigotEPSG = "EPSG::2969";
    41     private String SainteAnneEPSG = "EPSG::2970";
    42     private String MartiniqueEPSG = "EPSG::2973";
    43     private String ReunionEPSG = "EPSG::2975";
    44     private String GuyaneEPSG = "EPSG::2972";
    45     private String[] utmEPSGs = { FortMarigotEPSG, SainteAnneEPSG, MartiniqueEPSG, ReunionEPSG, GuyaneEPSG};
     44    private final static Integer FortMarigotEPSG = 2969;
     45    private final static Integer SainteAnneEPSG = 2970;
     46    private final static Integer MartiniqueEPSG = 2973;
     47    private final static Integer ReunionEPSG = 2975;
     48    private final static Integer GuyaneEPSG = 2972;
     49    private final static Integer[] utmEPSGs = { FortMarigotEPSG, SainteAnneEPSG, MartiniqueEPSG, ReunionEPSG, GuyaneEPSG };
    4650
    47     /**
    48      * false east in meters (constant)
    49      */
    50     private static final double Xs = 500000.0;
    51     /**
    52      * false north in meters (0 in northern hemisphere, 10000000 in southern
    53      * hemisphere)
    54      */
    55     private static double Ys = 0;
    56     /**
    57      * origin meridian longitude
    58      */
    59     protected double lg0;
     51    private final static Datum FortMarigotDatum = new ThreeParameterDatum("FortMarigot Datum", null, Ellipsoid.hayford, 136.596, 248.148, -429.789);
     52    private final static Datum SainteAnneDatum = new SevenParameterDatum("SainteAnne Datum", null, Ellipsoid.hayford, -472.29, -5.63, -304.12, 0.4362, -0.8374, 0.2563, 1.8984);
     53    private final static Datum MartiniqueDatum = new SevenParameterDatum("Martinique Datum", null, Ellipsoid.hayford, 126.926, 547.939, 130.409, -2.78670, 5.16124, -0.85844, 13.82265);
     54    private final static Datum ReunionDatum = GRS80Datum.INSTANCE;
     55    private final static Datum GuyaneDatum = GRS80Datum.INSTANCE;
     56    private final static Datum[] utmDatums = { FortMarigotDatum, SainteAnneDatum, MartiniqueDatum, ReunionDatum, GuyaneDatum };
     57
     58    private final static int[] utmZones = { 20, 20, 20, 40, 22 };
     59
    6060    /**
    6161     * UTM zone (from 1 to 60)
     
    6969    public static final int DEFAULT_GEODESIC = 0;
    7070
    71     public static int currentGeodesic = DEFAULT_GEODESIC;
     71    public int currentGeodesic;
    7272
    73     /**
    74      * 7 parameters transformation
    75      */
    76     private static double tx = 0.0;
    77     private static double ty = 0.0;
    78     private static double tz = 0.0;
    79     private static double rx = 0;
    80     private static double ry = 0;
    81     private static double rz = 0;
    82     private static double scaleDiff = 0;
    83     /**
    84      * precision in iterative schema
    85      */
    86     public static final double epsilon = 1e-11;
    8773
    88     private void refresh7ParametersTranslation() {
    89         if (currentGeodesic == 0) { // UTM_20N_Guadeloupe_Fort_Marigot
    90             set7ParametersTranslation(new double[]{136.596, 248.148, -429.789},
    91                     new double[]{0, 0, 0},
    92                     0,
    93                     true, 20);
    94         } else if (currentGeodesic == 1) { // UTM_20N_Guadeloupe_Ste_Anne
    95             set7ParametersTranslation(new double[]{-472.29, -5.63, -304.12},
    96                     new double[]{0.4362, -0.8374, 0.2563},
    97                     1.8984E-6,
    98                     true, 20);
    99         } else if (currentGeodesic == 2) { // UTM_20N_Martinique_Fort_Desaix
    100             set7ParametersTranslation(new double[]{126.926, 547.939, 130.409},
    101                     new double[]{-2.78670, 5.16124,  -0.85844},
    102                     13.82265E-6
    103                     , true, 20);
    104         } else if (currentGeodesic == 3) { // UTM_40S_Reunion_RGR92 (translation only required for re-projections from Gauss-Laborde)
    105             set7ParametersTranslation(new double[]{789.524, -626.486, -89.904},
    106                     new double[]{0.6006, 76.7946,  -10.5788},
    107                     -32.3241E-6
    108                     , false, 40);
    109         } else if (currentGeodesic == 4) { // UTM_22N_Guyane_RGFG95 (translation only required for re-projections from CSG67)
    110             set7ParametersTranslation(new double[]{-193.066 , 236.993, 105.447},
    111                     new double[]{0.4814, -0.8074,  0.1276},
    112                     1.5649E-6
    113                     , true, 22);
    114         }
     74    public UTM_France_DOM() {
     75        updateParameters(DEFAULT_GEODESIC);
     76    }
     77   
     78    public void updateParameters(int currentGeodesic) {
     79        this.currentGeodesic = currentGeodesic;
     80        datum = utmDatums[currentGeodesic];
     81        ellps = datum.getEllipsoid();
     82        proj = new org.openstreetmap.josm.data.projection.proj.TransverseMercator(ellps);
     83        isNorth = currentGeodesic != 3;
     84        zone = utmZones[currentGeodesic];
     85        x_0 = 500000;
     86        y_0 = isNorth ? 0.0 : 10000000.0;
     87        lon_0 = 6 * zone - 183;
     88        k_0 = 0.9996;
     89    }
     90   
     91    public int getCurrentGeodesic() {
     92        return currentGeodesic;
    11593    }
    11694
    117     private void set7ParametersTranslation(double[] translation, double[] rotation, double scalediff, boolean north, int utmZone) {
    118         tx = translation[0];
    119         ty = translation[1];
    120         tz = translation[2];
    121         rx = rotation[0]/206264.806247096355; // seconds to radian
    122         ry = rotation[1]/206264.806247096355;
    123         rz = rotation[2]/206264.806247096355;
    124         scaleDiff = scalediff;
    125         isNorth = north;
    126         Ys = isNorth ? 0.0 : 10000000.0;
    127         zone = utmZone;
     95    @Override
     96    public String toString() {
     97        return tr("UTM France (DOM)");
    12898    }
    12999
    130     public EastNorth latlon2eastNorth(LatLon p) {
    131         if (currentGeodesic < 3 ) {
    132             // translate ellipsoid GRS80 (WGS83) => reference ellipsoid geographic
    133             LatLon geo = GRS802Hayford(p);
    134             // reference ellipsoid geographic => UTM projection
    135             return MTProjection(geo, Ellipsoid.hayford.a, Ellipsoid.hayford.e);
    136         } else { // UTM_40S_Reunion_RGR92 or UTM_22N_Guyane_RGFG95
    137             LatLon geo = new LatLon(Math.toRadians(p.lat()), Math.toRadians(p.lon()));
    138             return MTProjection(geo, Ellipsoid.GRS80.a, Ellipsoid.GRS80.e);
    139         }
    140     }
    141 
    142     /**
    143      * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to UTM
    144      * geographic, (ellipsoid Hayford)
    145      */
    146     private LatLon GRS802Hayford(LatLon wgs) {
    147         double lat = Math.toRadians(wgs.lat()); // degree to radian
    148         double lon = Math.toRadians(wgs.lon());
    149         // WGS84 geographic => WGS84 cartesian
    150         double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
    151         double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
    152         double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
    153         double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
    154         // translation
    155         double coord[] = invSevenParametersTransformation(X, Y, Z);
    156         // UTM cartesian => UTM geographic
    157         return Geographic(coord[0], coord[1], coord[2], Ellipsoid.hayford);
    158     }
    159 
    160     /**
    161      * initializes from cartesian coordinates
    162      *
    163      * @param X
    164      *            1st coordinate in meters
    165      * @param Y
    166      *            2nd coordinate in meters
    167      * @param Z
    168      *            3rd coordinate in meters
    169      * @param ell
    170      *            reference ellipsoid
    171      */
    172     private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
    173         double norm = Math.sqrt(X * X + Y * Y);
    174         double lg = 2.0 * Math.atan(Y / (X + norm));
    175         double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
    176         double delta = 1.0;
    177         while (delta > epsilon) {
    178             double s2 = Math.sin(lt);
    179             s2 *= s2;
    180             double l = Math.atan((Z / norm)
    181                     / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
    182             delta = Math.abs(l - lt);
    183             lt = l;
    184         }
    185         double s2 = Math.sin(lt);
    186         s2 *= s2;
    187         // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
    188         return new LatLon(lt, lg);
    189     }
    190 
    191     /**
    192      * initalizes from geographic coordinates
    193      *
    194      * @param coord geographic coordinates triplet
    195      * @param a reference ellipsoid long axis
    196      * @param e reference ellipsoid excentricity
    197      */
    198     private EastNorth MTProjection(LatLon coord, double a, double e) {
    199         double n = 0.9996 * a;
    200         Ys = (coord.lat() >= 0.0) ? 0.0 : 10000000.0;
    201         double r6d = Math.PI / 30.0;
    202         //zone = (int) Math.floor((coord.lon() + Math.PI) / r6d) + 1;
    203         lg0 = r6d * (zone - 0.5) - Math.PI;
    204         double e2 = e * e;
    205         double e4 = e2 * e2;
    206         double e6 = e4 * e2;
    207         double e8 = e4 * e4;
    208         double C[] = {
    209                 1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
    210                 e2/8.0 - e4/96.0 - 9.0*e6/1024.0 - 901.0*e8/184320.0,
    211                 13.0*e4/768.0 + 17.0*e6/5120.0 - 311.0*e8/737280.0,
    212                 61.0*e6/15360.0 + 899.0*e8/430080.0,
    213                 49561.0*e8/41287680.0
    214         };
    215         double s = e * Math.sin(coord.lat());
    216         double l = Math.log(Math.tan(Math.PI/4.0 + coord.lat()/2.0) *
    217                 Math.pow((1.0 - s) / (1.0 + s), e/2.0));
    218         double phi = Math.asin(Math.sin(coord.lon() - lg0) /
    219                 ((Math.exp(l) + Math.exp(-l)) / 2.0));
    220         double ls = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
    221         double lambda = Math.atan(((Math.exp(l) - Math.exp(-l)) / 2.0) /
    222                 Math.cos(coord.lon() - lg0));
    223 
    224         double north = C[0] * lambda;
    225         double east = C[0] * ls;
    226         for(int k = 1; k < 5; k++) {
    227             double r = 2.0 * k * lambda;
    228             double m = 2.0 * k * ls;
    229             double em = Math.exp(m);
    230             double en = Math.exp(-m);
    231             double sr = Math.sin(r)/2.0 * (em + en);
    232             double sm = Math.cos(r)/2.0 * (em - en);
    233             north += C[k] * sr;
    234             east += C[k] * sm;
    235         }
    236         east *= n;
    237         east += Xs;
    238         north *= n;
    239         north += Ys;
    240         return new EastNorth(east, north);
    241     }
    242 
    243     public LatLon eastNorth2latlon(EastNorth p) {
    244         if (currentGeodesic < 3) {
    245             MTProjection(p.east(), p.north(), zone, isNorth);
    246             LatLon geo = Geographic(p, Ellipsoid.hayford.a, Ellipsoid.hayford.e, 0.0 /* z */);
    247 
    248             // reference ellipsoid geographic => reference ellipsoid cartesian
    249             double N = Ellipsoid.hayford.a / (Math.sqrt(1.0 - Ellipsoid.hayford.e2 * Math.sin(geo.lat()) * Math.sin(geo.lat())));
    250             double X = (N /*+ h*/) * Math.cos(geo.lat()) * Math.cos(geo.lon());
    251             double Y = (N /*+ h*/) * Math.cos(geo.lat()) * Math.sin(geo.lon());
    252             double Z = (N * (1.0-Ellipsoid.hayford.e2) /*+ h*/) * Math.sin(geo.lat());
    253             // translation
    254             double coord[] = sevenParametersTransformation(X, Y, Z);
    255             // WGS84 cartesian => WGS84 geographic
    256             LatLon wgs = cart2LatLon(coord[0], coord[1], coord[2], Ellipsoid.GRS80);
    257             return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
    258         } else {
    259             // UTM_40S_Reunion_RGR92 or UTM_22N_Guyane_RGFG95
    260             LatLon geo = Geographic(p, Ellipsoid.GRS80.a, Ellipsoid.GRS80.e, 0.0 /* z */);
    261             double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(geo.lat()) * Math.sin(geo.lat())));
    262             double X = (N /*+ h*/) * Math.cos(geo.lat()) * Math.cos(geo.lon());
    263             double Y = (N /*+ h*/) * Math.cos(geo.lat()) * Math.sin(geo.lon());
    264             double Z = (N * (1.0-Ellipsoid.GRS80.e2) /*+ h*/) * Math.sin(geo.lat());
    265             LatLon wgs = cart2LatLon(X, Y, Z, Ellipsoid.GRS80);
    266             return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
    267         }
    268     }
    269 
    270     /**
    271      * initializes new projection coordinates (in north hemisphere)
    272      *
    273      * @param east east from origin in meters
    274      * @param north north from origin in meters
    275      * @param zone zone number (from 1 to 60)
    276      * @param isNorth true in north hemisphere, false in south hemisphere
    277      */
    278     private void MTProjection(double east, double north, int zone, boolean isNorth) {
    279         Ys = isNorth ? 0.0 : 10000000.0;
    280         double r6d = Math.PI / 30.0;
    281         lg0 = r6d * (zone - 0.5) - Math.PI;
    282     }
    283 
    284     public double scaleFactor() {
    285         return 1/Math.PI/2;
    286     }
    287 
    288     /**
    289      * initalizes from projected coordinates (Mercator transverse projection)
    290      *
    291      * @param coord projected coordinates pair
    292      * @param e reference ellipsoid excentricity
    293      * @param a reference ellipsoid long axis
    294      * @param z altitude in meters
    295      */
    296     private LatLon Geographic(EastNorth coord, double a, double e, double z) {
    297         double n = 0.9996 * a;
    298         double e2 = e * e;
    299         double e4 = e2 * e2;
    300         double e6 = e4 * e2;
    301         double e8 = e4 * e4;
    302         double C[] = {
    303                 1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
    304                 e2/8.0 + e4/48.0 + 7.0*e6/2048.0 + e8/61440.0,
    305                 e4/768.0 + 3.0*e6/1280.0 + 559.0*e8/368640.0,
    306                 17.0*e6/30720.0 + 283.0*e8/430080.0,
    307                 4397.0*e8/41287680.0
    308         };
    309         double l = (coord.north() - Ys) / (n * C[0]);
    310         double ls = (coord.east() - Xs) / (n * C[0]);
    311         double l0 = l;
    312         double ls0 = ls;
    313         for(int k = 1; k < 5; k++) {
    314             double r = 2.0 * k * l0;
    315             double m = 2.0 * k * ls0;
    316             double em = Math.exp(m);
    317             double en = Math.exp(-m);
    318             double sr = Math.sin(r)/2.0 * (em + en);
    319             double sm = Math.cos(r)/2.0 * (em - en);
    320             l -= C[k] * sr;
    321             ls -= C[k] * sm;
    322         }
    323         double lon = lg0 + Math.atan(((Math.exp(ls) - Math.exp(-ls)) / 2.0) /
    324                 Math.cos(l));
    325         double phi = Math.asin(Math.sin(l) /
    326                 ((Math.exp(ls) + Math.exp(-ls)) / 2.0));
    327         l = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
    328         double lat = 2.0 * Math.atan(Math.exp(l)) - Math.PI / 2.0;
    329         double lt0;
    330         do {
    331             lt0 = lat;
    332             double s = e * Math.sin(lat);
    333             lat = 2.0 * Math.atan(Math.pow((1.0 + s) / (1.0 - s), e/2.0) *
    334                     Math.exp(l)) - Math.PI / 2.0;
    335         }
    336         while(Math.abs(lat-lt0) >= epsilon);
    337         //h = z;
    338 
    339         return new LatLon(lat, lon);
    340     }
    341 
    342     /**
    343      * initializes from cartesian coordinates
    344      *
    345      * @param X 1st coordinate in meters
    346      * @param Y 2nd coordinate in meters
    347      * @param Z 3rd coordinate in meters
    348      * @param ell reference ellipsoid
    349      */
    350     private LatLon cart2LatLon(double X, double Y, double Z, Ellipsoid ell) {
    351         double[] XYZ = {X, Y, Z};
    352         LatLon coord = ell.cart2LatLon(XYZ, epsilon);
    353         return new LatLon(Math.toRadians(coord.lat()), Math.toRadians(coord.lon()));
    354     }
    355 
    356     /**
    357      * 7 parameters transformation
    358      * @param coord X, Y, Z in array
    359      * @return transformed X, Y, Z in array
    360      */
    361     private double[] sevenParametersTransformation(double Xa, double Ya, double Za){
    362         double Xb = tx + Xa*(1+scaleDiff) + Za*ry - Ya*rz;
    363         double Yb = ty + Ya*(1+scaleDiff) + Xa*rz - Za*rx;
    364         double Zb = tz + Za*(1+scaleDiff) + Ya*rx - Xa*ry;
    365         return new double[]{Xb, Yb, Zb};
    366     }
    367 
    368     /**
    369      * 7 parameters inverse transformation
    370      * @param coord X, Y, Z in array
    371      * @return transformed X, Y, Z in array
    372      */
    373     private double[] invSevenParametersTransformation(double Xa, double Ya, double Za){
    374         double Xb = (1-scaleDiff)*(-tx + Xa + ((-tz+Za)*(-ry) - (-ty+Ya)*(-rz)));
    375         double Yb = (1-scaleDiff)*(-ty + Ya + ((-tx+Xa)*(-rz) - (-tz+Za)*(-rx)));
    376         double Zb = (1-scaleDiff)*(-tz + Za + ((-ty+Ya)*(-rx) - (-tx+Xa)*(-ry)));
    377         return new double[]{Xb, Yb, Zb};
    378     }
    379 
     100    @Override
    380101    public String getCacheDirectoryName() {
    381102        return this.toString();
    382103    }
    383104
    384     /**
    385      * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
    386      */
    387     public double getDefaultZoomInPPD() {
    388         // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
    389         return 10.0;
    390     }
    391 
     105    @Override
    392106    public Bounds getWorldBoundsLatLon() {
    393107        return utmBounds[currentGeodesic];
    394108    }
    395109
    396     public String toCode() {
     110    @Override
     111    public Integer getEpsgCode() {
    397112        return utmEPSGs[currentGeodesic];
    398113    }
     
    401116    public int hashCode() {
    402117        return getClass().getName().hashCode()+currentGeodesic; // our only real variable
    403     }
    404 
    405     @Override public String toString() {
    406         return (tr("UTM France (DOM)"));
    407     }
    408 
    409     public int getCurrentGeodesic() {
    410         return currentGeodesic;
    411118    }
    412119
     
    426133    }
    427134
     135    @Override
    428136    public Collection<String> getPreferences(JPanel p) {
    429137        Object prefcb = p.getComponent(2);
     
    431139            return null;
    432140        currentGeodesic = ((JComboBox)prefcb).getSelectedIndex();
    433         refresh7ParametersTranslation();
    434141        return Collections.singleton(Integer.toString(currentGeodesic+1));
    435142    }
     
    437144    @Override
    438145    public String[] allCodes() {
    439         return utmEPSGs;
     146        String[] res = new String[utmEPSGs.length];
     147        for (int i=0; i<utmEPSGs.length; ++i) {
     148            res[i] = "EPSG:"+utmEPSGs[i];
     149        }
     150        return res;
    440151    }
    441152
     153    @Override
    442154    public Collection<String> getPreferencesFromCode(String code) {
    443155        for (int i=0; i < utmEPSGs.length; i++ )
    444             if (utmEPSGs[i].endsWith(code))
     156            if (("EPSG:"+utmEPSGs[i]).equals(code))
    445157                return Collections.singleton(Integer.toString(i+1));
    446158        return null;
    447159    }
    448160
     161    @Override
    449162    public void setPreferences(Collection<String> args) {
    450         currentGeodesic = DEFAULT_GEODESIC;
     163        int currentGeodesic = DEFAULT_GEODESIC;
    451164        if (args != null) {
    452165            try {
     
    461174            } catch(NumberFormatException e) {}
    462175        }
    463         refresh7ParametersTranslation();
     176        updateParameters(currentGeodesic);
    464177    }
    465 
    466178}
  • trunk/src/org/openstreetmap/josm/data/projection/proj/TransverseMercator.java

    r4281 r4285  
    11// License: GPL. For details, see LICENSE file.
    2 package org.openstreetmap.josm.data.projection;
    3 
    4 import org.openstreetmap.josm.data.coor.EastNorth;
    5 import org.openstreetmap.josm.data.coor.LatLon;
     2package org.openstreetmap.josm.data.projection.proj;
     3
     4import static java.lang.Math.*;
     5
     6import static org.openstreetmap.josm.tools.I18n.tr;
     7
     8import org.openstreetmap.josm.data.projection.Ellipsoid;
    69
    710/**
    8  * This is a base class to do projections based on Transverse Mercator projection.
     11 * Transverse Mercator projection.
    912 *
    1013 * @author Dirk Stöcker
    1114 * code based on JavaScript from Chuck Taylor
    1215 *
    13  * NOTE: Uses polygon approximation to translate to WGS84.
    1416 */
    15 public abstract class TransverseMercator implements Projection {
    16 
    17     private final static double UTMScaleFactor = 0.9996;
    18 
    19     private double UTMCentralMeridianRad = 0;
    20     private double offsetEastMeters = 500000;
    21     private double offsetNorthMeters = 0;
    22 
    23 
    24     protected void setProjectionParameters(double centralMeridianDegress, double offsetEast, double offsetNorth)
    25     {
    26         UTMCentralMeridianRad = Math.toRadians(centralMeridianDegress);
    27         offsetEastMeters = offsetEast;
    28         offsetNorthMeters = offsetNorth;
    29     }
    30 
    31     /*
    32      * ArcLengthOfMeridian
    33      *
    34      * Computes the ellipsoidal distance from the equator to a point at a
    35      * given latitude.
    36      *
    37      * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
    38      * GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
    39      *
    40      * Inputs:
    41      *     phi - Latitude of the point, in radians.
    42      *
    43      * Globals:
    44      *     Ellipsoid.GRS80.a - Ellipsoid model major axis.
    45      *     Ellipsoid.GRS80.b - Ellipsoid model minor axis.
    46      *
    47      * Returns:
    48      *     The ellipsoidal distance of the point from the equator, in meters.
    49      *
    50      */
    51     private double ArcLengthOfMeridian(double phi)
    52     {
    53         /* Precalculate n */
    54         double n = (Ellipsoid.GRS80.a - Ellipsoid.GRS80.b) / (Ellipsoid.GRS80.a + Ellipsoid.GRS80.b);
    55 
    56         /* Precalculate alpha */
    57         double alpha = ((Ellipsoid.GRS80.a + Ellipsoid.GRS80.b) / 2.0)
    58         * (1.0 + (Math.pow (n, 2.0) / 4.0) + (Math.pow (n, 4.0) / 64.0));
    59 
    60         /* Precalculate beta */
    61         double beta = (-3.0 * n / 2.0) + (9.0 * Math.pow (n, 3.0) / 16.0)
    62         + (-3.0 * Math.pow (n, 5.0) / 32.0);
    63 
    64         /* Precalculate gamma */
    65         double gamma = (15.0 * Math.pow (n, 2.0) / 16.0)
    66         + (-15.0 * Math.pow (n, 4.0) / 32.0);
    67 
    68         /* Precalculate delta */
    69         double delta = (-35.0 * Math.pow (n, 3.0) / 48.0)
    70         + (105.0 * Math.pow (n, 5.0) / 256.0);
    71 
    72         /* Precalculate epsilon */
    73         double epsilon = (315.0 * Math.pow (n, 4.0) / 512.0);
    74 
    75         /* Now calculate the sum of the series and return */
    76         return alpha
    77         * (phi + (beta * Math.sin (2.0 * phi))
    78                 + (gamma * Math.sin (4.0 * phi))
    79                 + (delta * Math.sin (6.0 * phi))
    80                 + (epsilon * Math.sin (8.0 * phi)));
    81     }
    82 
    83     /*
    84      * FootpointLatitude
    85      *
    86      * Computes the footpoint latitude for use in converting transverse
    87      * Mercator coordinates to ellipsoidal coordinates.
    88      *
    89      * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
    90      *   GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
    91      *
    92      * Inputs:
    93      *   y - The UTM northing coordinate, in meters.
    94      *
    95      * Returns:
    96      *   The footpoint latitude, in radians.
    97      *
    98      */
    99     private double FootpointLatitude(double y)
    100     {
    101         /* Precalculate n (Eq. 10.18) */
    102         double n = (Ellipsoid.GRS80.a - Ellipsoid.GRS80.b) / (Ellipsoid.GRS80.a + Ellipsoid.GRS80.b);
    103 
    104         /* Precalculate alpha_ (Eq. 10.22) */
    105         /* (Same as alpha in Eq. 10.17) */
    106         double alpha_ = ((Ellipsoid.GRS80.a + Ellipsoid.GRS80.b) / 2.0)
    107         * (1 + (Math.pow (n, 2.0) / 4) + (Math.pow (n, 4.0) / 64));
    108 
    109         /* Precalculate y_ (Eq. 10.23) */
    110         double y_ = y / alpha_;
    111 
    112         /* Precalculate beta_ (Eq. 10.22) */
    113         double beta_ = (3.0 * n / 2.0) + (-27.0 * Math.pow (n, 3.0) / 32.0)
    114         + (269.0 * Math.pow (n, 5.0) / 512.0);
    115 
    116         /* Precalculate gamma_ (Eq. 10.22) */
    117         double gamma_ = (21.0 * Math.pow (n, 2.0) / 16.0)
    118         + (-55.0 * Math.pow (n, 4.0) / 32.0);
    119 
    120         /* Precalculate delta_ (Eq. 10.22) */
    121         double delta_ = (151.0 * Math.pow (n, 3.0) / 96.0)
    122         + (-417.0 * Math.pow (n, 5.0) / 128.0);
    123 
    124         /* Precalculate epsilon_ (Eq. 10.22) */
    125         double epsilon_ = (1097.0 * Math.pow (n, 4.0) / 512.0);
    126 
    127         /* Now calculate the sum of the series (Eq. 10.21) */
    128         return y_ + (beta_ * Math.sin (2.0 * y_))
    129         + (gamma_ * Math.sin (4.0 * y_))
    130         + (delta_ * Math.sin (6.0 * y_))
    131         + (epsilon_ * Math.sin (8.0 * y_));
    132     }
    133 
    134     /*
    135      * MapLatLonToXY
    136      *
     17public class TransverseMercator implements Proj {
     18
     19    protected double a, b;
     20
     21    public TransverseMercator(Ellipsoid ellps) {
     22        this.a = ellps.a;
     23        this.b = ellps.b;
     24    }
     25   
     26    @Override
     27    public String getName() {
     28        return tr("Transverse Mercator");
     29    }
     30
     31    @Override
     32    public String getProj4Id() {
     33        return "tmerc";
     34    }
     35
     36    /**
    13737     * Converts a latitude/longitude pair to x and y coordinates in the
    13838     * Transverse Mercator projection.  Note that Transverse Mercator is not
     
    14141     * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
    14242     * GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
    143      *
    144      * Inputs:
    145      *    phi - Latitude of the point, in radians.
    146      *    lambda - Longitude of the point, in radians.
    147      *    lambda0 - Longitude of the central meridian to be used, in radians.
    148      *
    149      * Outputs:
    150      *    xy - A 2-element array containing the x and y coordinates
    151      *         of the computed point.
    152      *
    153      * Returns:
    154      *    The function does not return a value.
    155      *
    156      */
    157     public EastNorth mapLatLonToXY(double phi, double lambda, double lambda0)
    158     {
     43     *
     44     * @param phi Latitude of the point, in radians
     45     * @param lambda Longitude of the point, in radians
     46     * @return A 2-element array containing the x and y coordinates
     47     *         of the computed point
     48     */
     49    @Override
     50    public double[] project(double phi, double lambda) {
     51       
    15952        /* Precalculate ep2 */
    160         double ep2 = (Math.pow (Ellipsoid.GRS80.a, 2.0) - Math.pow (Ellipsoid.GRS80.b, 2.0)) / Math.pow (Ellipsoid.GRS80.b, 2.0);
     53        double ep2 = (pow(a, 2.0) - pow(b, 2.0)) / pow(b, 2.0);
    16154
    16255        /* Precalculate nu2 */
    163         double nu2 = ep2 * Math.pow (Math.cos (phi), 2.0);
    164 
    165         /* Precalculate N */
    166         double N = Math.pow (Ellipsoid.GRS80.a, 2.0) / (Ellipsoid.GRS80.b * Math.sqrt (1 + nu2));
     56        double nu2 = ep2 * pow(cos(phi), 2.0);
     57
     58        /* Precalculate N / a */
     59        double N_a = a / (b * sqrt(1 + nu2));
    16760
    16861        /* Precalculate t */
    169         double t = Math.tan (phi);
     62        double t = tan(phi);
    17063        double t2 = t * t;
    17164
    17265        /* Precalculate l */
    173         double l = lambda - lambda0;
     66        double l = lambda;
    17467
    17568        /* Precalculate coefficients for l**n in the equations below
     
    19184        double l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
    19285
    193         return new EastNorth(
     86        return new double[] {
    19487                /* Calculate easting (x) */
    195                 N * Math.cos (phi) * l
    196                 + (N / 6.0 * Math.pow (Math.cos (phi), 3.0) * l3coef * Math.pow (l, 3.0))
    197                 + (N / 120.0 * Math.pow (Math.cos (phi), 5.0) * l5coef * Math.pow (l, 5.0))
    198                 + (N / 5040.0 * Math.pow (Math.cos (phi), 7.0) * l7coef * Math.pow (l, 7.0)),
     88                N_a * cos(phi) * l
     89                + (N_a / 6.0 * pow(cos(phi), 3.0) * l3coef * pow(l, 3.0))
     90                + (N_a / 120.0 * pow(cos(phi), 5.0) * l5coef * pow(l, 5.0))
     91                + (N_a / 5040.0 * pow(cos(phi), 7.0) * l7coef * pow(l, 7.0)),
    19992                /* Calculate northing (y) */
    200                 ArcLengthOfMeridian (phi)
    201                 + (t / 2.0 * N * Math.pow (Math.cos (phi), 2.0) * Math.pow (l, 2.0))
    202                 + (t / 24.0 * N * Math.pow (Math.cos (phi), 4.0) * l4coef * Math.pow (l, 4.0))
    203                 + (t / 720.0 * N * Math.pow (Math.cos (phi), 6.0) * l6coef * Math.pow (l, 6.0))
    204                 + (t / 40320.0 * N * Math.pow (Math.cos (phi), 8.0) * l8coef * Math.pow (l, 8.0)));
    205     }
    206 
    207     /*
    208      * MapXYToLatLon
    209      *
     93                ArcLengthOfMeridian (phi) / a
     94                + (t / 2.0 * N_a * pow(cos(phi), 2.0) * pow(l, 2.0))
     95                + (t / 24.0 * N_a * pow(cos(phi), 4.0) * l4coef * pow(l, 4.0))
     96                + (t / 720.0 * N_a * pow(cos(phi), 6.0) * l6coef * pow(l, 6.0))
     97                + (t / 40320.0 * N_a * pow(cos(phi), 8.0) * l8coef * pow(l, 8.0)) };
     98    }
     99
     100    /**
    210101     * Converts x and y coordinates in the Transverse Mercator projection to
    211102     * a latitude/longitude pair.  Note that Transverse Mercator is not
     
    214105     * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
    215106     *   GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
    216      *
    217      * Inputs:
    218      *   x - The easting of the point, in meters.
    219      *   y - The northing of the point, in meters.
    220      *   lambda0 - Longitude of the central meridian to be used, in radians.
    221      *
    222      * Outputs:
    223      *   philambda - A 2-element containing the latitude and longitude
    224      *               in radians.
    225      *
    226      * Returns:
    227      *   The function does not return a value.
    228      *
     107     *
    229108     * Remarks:
    230109     *   The local variables Nf, nuf2, tf, and tf2 serve the same purpose as
     
    234113     *   x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
    235114     *   to optimize computations.
    236      *
    237      */
    238     public LatLon mapXYToLatLon(double x, double y, double lambda0)
    239     {
     115     *
     116     * @param x The easting of the point, in meters, divided by the semi major axis of the ellipsoid
     117     * @param y The northing of the point, in meters, divided by the semi major axis of the ellipsoid
     118     * @return A 2-element containing the latitude and longitude
     119     *               in radians
     120     */
     121    @Override
     122    public double[] invproject(double x, double y) {
    240123        /* Get the value of phif, the footpoint latitude. */
    241         double phif = FootpointLatitude (y);
     124        double phif = footpointLatitude(y);
    242125
    243126        /* Precalculate ep2 */
    244         double ep2 = (Math.pow (Ellipsoid.GRS80.a, 2.0) - Math.pow (Ellipsoid.GRS80.b, 2.0))
    245         / Math.pow (Ellipsoid.GRS80.b, 2.0);
    246 
    247         /* Precalculate cos (phif) */
    248         double cf = Math.cos (phif);
     127        double ep2 = (a*a - b*b)
     128        / (b*b);
     129
     130        /* Precalculate cos(phif) */
     131        double cf = cos(phif);
    249132
    250133        /* Precalculate nuf2 */
    251         double nuf2 = ep2 * Math.pow (cf, 2.0);
    252 
    253         /* Precalculate Nf and initialize Nfpow */
    254         double Nf = Math.pow (Ellipsoid.GRS80.a, 2.0) / (Ellipsoid.GRS80.b * Math.sqrt (1 + nuf2));
    255         double Nfpow = Nf;
     134        double nuf2 = ep2 * pow(cf, 2.0);
     135
     136        /* Precalculate Nf / a and initialize Nfpow */
     137        double Nf_a = a / (b * sqrt(1 + nuf2));
     138        double Nfpow = Nf_a;
    256139
    257140        /* Precalculate tf */
    258         double tf = Math.tan (phif);
     141        double tf = tan(phif);
    259142        double tf2 = tf * tf;
    260143        double tf4 = tf2 * tf2;
     
    264147        double x1frac = 1.0 / (Nfpow * cf);
    265148
    266         Nfpow *= Nf;   /* now equals Nf**2) */
     149        Nfpow *= Nf_a;   /* now equals Nf**2) */
    267150        double x2frac = tf / (2.0 * Nfpow);
    268151
    269         Nfpow *= Nf;   /* now equals Nf**3) */
     152        Nfpow *= Nf_a;   /* now equals Nf**3) */
    270153        double x3frac = 1.0 / (6.0 * Nfpow * cf);
    271154
    272         Nfpow *= Nf;   /* now equals Nf**4) */
     155        Nfpow *= Nf_a;   /* now equals Nf**4) */
    273156        double x4frac = tf / (24.0 * Nfpow);
    274157
    275         Nfpow *= Nf;   /* now equals Nf**5) */
     158        Nfpow *= Nf_a;   /* now equals Nf**5) */
    276159        double x5frac = 1.0 / (120.0 * Nfpow * cf);
    277160
    278         Nfpow *= Nf;   /* now equals Nf**6) */
     161        Nfpow *= Nf_a;   /* now equals Nf**6) */
    279162        double x6frac = tf / (720.0 * Nfpow);
    280163
    281         Nfpow *= Nf;   /* now equals Nf**7) */
     164        Nfpow *= Nf_a;   /* now equals Nf**7) */
    282165        double x7frac = 1.0 / (5040.0 * Nfpow * cf);
    283166
    284         Nfpow *= Nf;   /* now equals Nf**8) */
     167        Nfpow *= Nf_a;   /* now equals Nf**8) */
    285168        double x8frac = tf / (40320.0 * Nfpow);
    286169
     
    295178        double x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2);
    296179
    297         return new LatLon(
     180        return new double[] {
    298181                /* Calculate latitude */
    299                 Math.toDegrees(
    300182                        phif + x2frac * x2poly * (x * x)
    301                         + x4frac * x4poly * Math.pow (x, 4.0)
    302                         + x6frac * x6poly * Math.pow (x, 6.0)
    303                         + x8frac * x8poly * Math.pow (x, 8.0)),
    304                         Math.toDegrees(
    305                                 /* Calculate longitude */
    306                                 lambda0 + x1frac * x
    307                                 + x3frac * x3poly * Math.pow (x, 3.0)
    308                                 + x5frac * x5poly * Math.pow (x, 5.0)
    309                                 + x7frac * x7poly * Math.pow (x, 7.0)));
    310     }
    311 
    312     @Override
    313     public EastNorth latlon2eastNorth(LatLon p) {
    314         EastNorth a = mapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), UTMCentralMeridianRad);
    315         return new EastNorth(a.east() * UTMScaleFactor + offsetEastMeters, a.north() * UTMScaleFactor + offsetNorthMeters);
    316     }
    317 
    318     @Override
    319     public LatLon eastNorth2latlon(EastNorth p) {
    320         return mapXYToLatLon((p.east() - offsetEastMeters)/UTMScaleFactor, (p.north() - offsetNorthMeters)/UTMScaleFactor, UTMCentralMeridianRad);
    321     }
    322 
    323     @Override
    324     public double getDefaultZoomInPPD() {
    325         // this will set the map scaler to about 1000 m
    326         return 10;
    327     }
     183                        + x4frac * x4poly * pow(x, 4.0)
     184                        + x6frac * x6poly * pow(x, 6.0)
     185                        + x8frac * x8poly * pow(x, 8.0),
     186                        /* Calculate longitude */
     187                        x1frac * x
     188                        + x3frac * x3poly * pow(x, 3.0)
     189                        + x5frac * x5poly * pow(x, 5.0)
     190                        + x7frac * x7poly * pow(x, 7.0) };
     191    }
     192   
     193    /**
     194     * ArcLengthOfMeridian
     195     *
     196     * Computes the ellipsoidal distance from the equator to a point at a
     197     * given latitude.
     198     *
     199     * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
     200     * GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
     201     *
     202     * @param phi Latitude of the point, in radians
     203     * @return The ellipsoidal distance of the point from the equator
     204     *         (in meters, divided by the semi major axis of the ellipsoid)
     205     */
     206    private double ArcLengthOfMeridian(double phi) {
     207        /* Precalculate n */
     208        double n = (a - b) / (a + b);
     209
     210        /* Precalculate alpha */
     211        double alpha = ((a + b) / 2.0)
     212            * (1.0 + (pow(n, 2.0) / 4.0) + (pow(n, 4.0) / 64.0));
     213
     214        /* Precalculate beta */
     215        double beta = (-3.0 * n / 2.0) + (9.0 * pow(n, 3.0) / 16.0)
     216            + (-3.0 * pow(n, 5.0) / 32.0);
     217
     218        /* Precalculate gamma */
     219        double gamma = (15.0 * pow(n, 2.0) / 16.0)
     220            + (-15.0 * pow(n, 4.0) / 32.0);
     221
     222        /* Precalculate delta */
     223        double delta = (-35.0 * pow(n, 3.0) / 48.0)
     224            + (105.0 * pow(n, 5.0) / 256.0);
     225
     226        /* Precalculate epsilon */
     227        double epsilon = (315.0 * pow(n, 4.0) / 512.0);
     228
     229        /* Now calculate the sum of the series and return */
     230        return alpha
     231            * (phi + (beta * sin(2.0 * phi))
     232                    + (gamma * sin(4.0 * phi))
     233                    + (delta * sin(6.0 * phi))
     234                    + (epsilon * sin(8.0 * phi)));
     235    }
     236       
     237    /**
     238     * FootpointLatitude
     239     *
     240     * Computes the footpoint latitude for use in converting transverse
     241     * Mercator coordinates to ellipsoidal coordinates.
     242     *
     243     * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
     244     *   GPS: Theory and Practice, 3rd ed.  New York: Springer-Verlag Wien, 1994.
     245     *
     246     * @param y northing coordinate, in meters, divided by the semi major axis of the ellipsoid
     247     * @return The footpoint latitude, in radians
     248     */
     249    private double footpointLatitude(double y) {
     250        /* Precalculate n (Eq. 10.18) */
     251        double n = (a - b) / (a + b);
     252
     253        /* Precalculate alpha_ (Eq. 10.22) */
     254        /* (Same as alpha in Eq. 10.17) */
     255        double alpha_ = ((a + b) / 2.0)
     256            * (1 + (pow(n, 2.0) / 4) + (pow(n, 4.0) / 64));
     257
     258        /* Precalculate y_ (Eq. 10.23) */
     259        double y_ = y / alpha_ * a;
     260
     261        /* Precalculate beta_ (Eq. 10.22) */
     262        double beta_ = (3.0 * n / 2.0) + (-27.0 * pow(n, 3.0) / 32.0)
     263            + (269.0 * pow(n, 5.0) / 512.0);
     264
     265        /* Precalculate gamma_ (Eq. 10.22) */
     266        double gamma_ = (21.0 * pow(n, 2.0) / 16.0)
     267            + (-55.0 * pow(n, 4.0) / 32.0);
     268
     269        /* Precalculate delta_ (Eq. 10.22) */
     270        double delta_ = (151.0 * pow(n, 3.0) / 96.0)
     271            + (-417.0 * pow(n, 5.0) / 128.0);
     272
     273        /* Precalculate epsilon_ (Eq. 10.22) */
     274        double epsilon_ = (1097.0 * pow(n, 4.0) / 512.0);
     275
     276        /* Now calculate the sum of the series (Eq. 10.21) */
     277        return y_ + (beta_ * sin(2.0 * y_))
     278            + (gamma_ * sin(4.0 * y_))
     279            + (delta_ * sin(6.0 * y_))
     280            + (epsilon_ * sin(8.0 * y_));
     281    }
     282   
    328283}
Note: See TracChangeset for help on using the changeset viewer.