Ignore:
Timestamp:
2011-08-07T12:17:39+02:00 (13 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/proj
Files:
1 added
1 moved

Legend:

Unmodified
Added
Removed

  • 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.