Ignore:
Timestamp:
2008-06-25T01:49:12+02:00 (16 years ago)
Author:
framm
Message:
  • fix in Lambert projection by Pieren <pieren3@…>
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r656 r657  
    55package org.openstreetmap.josm.data.projection;
    66
     7import static org.openstreetmap.josm.tools.I18n.tr;
     8
     9import javax.swing.JOptionPane;
     10
     11import org.openstreetmap.josm.Main;
    712import org.openstreetmap.josm.data.coor.EastNorth;
    813import org.openstreetmap.josm.data.coor.LatLon;
    … …  
    1217     * Lambert I, II, III, and IV projection exponents
    1318     */
    14     public static final double n[] = {
    15         0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322
    16     };
     19    public static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 };
     20
    1721    /**
    1822     * Lambert I, II, III, and IV projection constants
    1923     */
    20     public static final double c[] = {
    21         11603796.98, 11745793.39, 11947992.52, 12136281.99
    22     };
     24    public static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 };
     25
    2326    /**
    2427     * Lambert I, II, III, and IV false east
    2528     */
    26     public static final double Xs[] = {
    27         600000.0, 600000.0, 600000.0, 234.358
    28     };
     29    public static final double Xs[] = { 600000.0, 600000.0, 600000.0, 234.358 };
     30
    2931    /**
    3032     * Lambert I, II, III, and IV false north
    3133     */
    32     public static final double Ys[] = {
    33         5657616.674, 6199695.768, 6791905.085, 7239161.542
    34     };
     34    public static final double Ys[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 };
     35
    3536    /**
    3637     * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian
    3738     */
    38     public static final double lg0 = 0.04079234433198; // 2 deg 20 min 14.025 sec
     39    public static final double lg0 = 0.04079234433198; // 2deg20'14.025"
     40
    3941    /**
    4042     * precision in iterative schema
    4143     */
     44
    4245    public static final double epsilon = 1e-11;
    4346
    44     public static int zone = 0;
    45 
     47    /**
     48     * France is divided in 4 Lambert projection zones (1,2,3 + 4th for Corsica)
     49     */
     50    public static final double cMaxLatZone1 = Math.toRadians(57 * 0.9);
     51
     52    public static final double zoneLimits[] = { Math.toRadians(53.5 * 0.9), // between Zone 1 and Zone 2 (in grad *0.9)
     53            Math.toRadians(50.5 * 0.9), // between Zone 2 and Zone 3
     54            Math.toRadians(47.51963 * 0.9), // between Zone 3 and Zone 4
     55            Math.toRadians(46.17821 * 0.9) };// lowest latitude of Zone 4
     56
     57    public static final double cMinLonZones = Math.toRadians(-4.9074074074074059 * 0.9);
     58
     59    public static final double cMaxLonZones = Math.toRadians(10.2 * 0.9);
     60
     61    /**
     62     *  Because josm cannot work correctly if two zones are displayed, we allow some overlapping
     63     */
     64    public static final double cMaxOverlappingZones = Math.toRadians(0.5 * 0.9);
     65
     66    public static int layoutZone = -1;
     67
     68    /**
     69     * @param p
     70     *            a WGS84 lat/lon (ellipsoid GRS80) (in degree)
     71     * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
     72     */
    4673    public EastNorth latlon2eastNorth(LatLon p) {
     74        // translate ellipsoid GRS80 (WGS83) => Clark
     75        LatLon geo = GRS802Clark(p);
     76        double lt = geo.lat(); // in radian
     77        double lg = geo.lon();
     78
    4779        // check if longitude and latitude are inside the french Lambert zones
    48         double lt = p.lat();
    49         double lg = p.lon();
    50         if(lt > 57*9/10 || lt < 46.17821*9/10 ||
    51                 lg > 10.2*9/10 || lg < -4.9074074074074059*9/10) {
    52             if (lt > 57*9/10) {
    53                 // out of Lambert zones, possible when MAX_LAT is used- keep current zone
     80        int currentZone = 0;
     81        boolean outOfLambertZones = false;
     82        if (lt >= zoneLimits[3] && lt <= cMaxLatZone1 && lg >= cMinLonZones && lg <= cMaxLonZones) {
     83            // zone I
     84            if (lt > zoneLimits[0])
     85                currentZone = 0;
     86            // zone II
     87            else if (lt > zoneLimits[1])
     88                currentZone = 1;
     89            // zone III
     90            else if (lt > zoneLimits[2])
     91                currentZone = 2;
     92            // zone III or IV
     93            else if (lt > zoneLimits[3])
     94                // Note: zone IV is dedicated to Corsica island and extends from 47.8 to
     95                // 45.9 degrees of latitude. There is an overlap with zone III that can be
     96                // solved only with longitude (covers Corsica if lon > 7.2 degree)
     97                if (lg < Math.toRadians(8 * 0.9))
     98                    currentZone = 2;
     99                else
     100                    currentZone = 3;
     101        } else {
     102            outOfLambertZones = true; // possible when MAX_LAT is used
     103        }
     104        if (layoutZone == -1)
     105            layoutZone = currentZone;
     106        else if (layoutZone != currentZone && !outOfLambertZones) {
     107            if ((currentZone < layoutZone && Math.abs(zoneLimits[currentZone] - lt) > cMaxOverlappingZones)
     108                    || (currentZone > layoutZone && Math.abs(zoneLimits[layoutZone] - lt) > cMaxOverlappingZones)) {
     109                JOptionPane.showMessageDialog(Main.parent,
     110                        tr("Some data are positionned far away from current Lambert zone limits.\n"
     111                                +"Split long ways to avoid distortions."));
     112                layoutZone = -1;
     113            } else {
     114                System.out.println("temporarily extends Lambert zone " + layoutZone
     115                                + " projection at lat,lon:" + lt + "," + lg);
    54116            }
    55             // zone I
    56             else if (lt > 53.5*9/10)
    57                 zone = 0;
    58             // zone II
    59             else if(lt > 50.5*9/10)
    60                 zone = 1;
    61             // zone III
    62             else if(lt > 47.51963*9/10)
    63                 zone = 2;
    64             else if (lt > 46.17821*9/10)
    65                 // zone III or IV
    66                 // Note: zone IV is dedicated to Corsica island and extends
    67                 // from 47.8 to 45.9 degrees of latitude. There is an overlap with
    68                 // zone III that can  be solved only with longitude
    69                 // (covers Corsica if lon > 7.2 degree)
    70                 if (lg < 8*9/10)
    71                     zone = 2;
    72                 else
    73                     zone = 3;
    74         } else {
    75             // out of Lambert zones- keep current one
    76         }
    77         EastNorth eastNorth = ConicProjection(lt, lg, Xs[zone], Ys[zone], c[zone], n[zone]);
    78         return eastNorth;
     117        }
     118        if (layoutZone == -1) {
     119            return ConicProjection(lt, lg, Xs[currentZone], Ys[currentZone], c[currentZone], n[currentZone]);
     120        } // else
     121        return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
    79122    }
    80123
    81124    public LatLon eastNorth2latlon(EastNorth p) {
    82         return Geographic(p, Xs[zone], Ys[zone], c[zone], n[zone]);
    83     }
    84 
    85     @Override public String toString() {
     125        LatLon geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
     126        // translate ellipsoid Clark => GRS80 (WGS83)
     127        LatLon wgs = Clark2GRS80(geo);
     128        return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
     129    }
     130
     131    @Override
     132    public String toString() {
    86133        return "Lambert";
    87134    }
    … …  
    92139
    93140    public double scaleFactor() {
    94         return 1.0/360;
    95     }
    96 
    97     @Override public boolean equals(Object o) {
     141        return 1.0 / 360;
     142    }
     143
     144    @Override
     145    public boolean equals(Object o) {
    98146        return o instanceof Lambert;
    99147    }
    100148
    101     @Override public int hashCode() {
     149    @Override
     150    public int hashCode() {
    102151        return Lambert.class.hashCode();
    103152    }
    104153
    105154    /**
    106      * Initializes from geographic coordinates.
    107      * Note that reference ellipsoid used by Lambert is the Clark ellipsoid.
    108      *
    109      * @param lat latitude in degree
    110      * @param lon longitude in degree
    111      * @param Xs false east (coordinate system origin) in meters
    112      * @param Ys false north (coordinate system origin) in meters
    113      * @param c projection constant
    114      * @param n projection exponent
    115      * @return EastNorth projected coordinates in meter
    116      */
    117     public EastNorth ConicProjection(double lat, double lon, double Xs, double Ys,
    118             double c, double n) {
    119         lat = Math.toRadians(lat);
    120         lon = Math.toRadians(lon);
     155     * Initializes from geographic coordinates. Note that reference ellipsoid
     156     * used by Lambert is the Clark ellipsoid.
     157     *
     158     * @param lat latitude in grad
     159     * @param lon longitude in grad
     160     * @param Xs  false east (coordinate system origin) in meters
     161     * @param Ys  false north (coordinate system origin) in meters
     162     * @param c   projection constant
     163     * @param n   projection exponent
     164     * @return EastNorth projected coordinates in meter
     165     */
     166    private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) {
    121167        double eslt = Ellipsoid.clarke.e * Math.sin(lat);
    122         double l = Math.log(Math.tan(Math.PI/4.0 + (lat/2.0)) *
    123                 Math.pow((1.0 - eslt)/(1.0 + eslt), Ellipsoid.clarke.e/2.0));
     168        double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0))
     169                * Math.pow((1.0 - eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0));
    124170        double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0));
    125171        double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0));
    126172        return new EastNorth(east, north);
    127173    }
    128     /**
    129      * Initializes from projected coordinates (conic projection).
    130      * Note that reference ellipsoid used by Lambert is Clark
    131      *
     174
     175    /**
     176     * Initializes from projected coordinates (conic projection). Note that
     177     * reference ellipsoid used by Lambert is Clark
     178     *
    132179     * @param coord projected coordinates pair in meters
    133      * @param Xs false east (coordinate system origin) in meters
    134      * @param Ys false north (coordinate system origin) in meters
    135      * @param c projection constant
    136      * @param n projection exponent
    137      * @return LatLon in degrees
    138      */
    139     public LatLon Geographic(EastNorth eastNorth, double Xs, double Ys,
    140             double c, double n) {
     180     * @param Xs    false east (coordinate system origin) in meters
     181     * @param Ys    false north (coordinate system origin) in meters
     182     * @param c     projection constant
     183     * @param n     projection exponent
     184     * @return LatLon in radian
     185     */
     186    private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) {
    141187        double dx = eastNorth.east() - Xs;
    142188        double dy = Ys - eastNorth.north();
    143         double R = Math.sqrt(dx*dx + dy*dy);
     189        double R = Math.sqrt(dx * dx + dy * dy);
    144190        double gamma = Math.atan(dx / dy);
    145         double l = -1.0/n * Math.log(Math.abs(R / c));
     191        double l = -1.0 / n * Math.log(Math.abs(R / c));
    146192        l = Math.exp(l);
    147193        double lon = lg0 + gamma / n;
    148         double lat = 2.0 * Math.atan(l) - Math.PI/2.0;
     194        double lat = 2.0 * Math.atan(l) - Math.PI / 2.0;
    149195        double delta = 1.0;
    150         while(delta > epsilon) {
     196        while (delta > epsilon) {
    151197            double eslt = Ellipsoid.clarke.e * Math.sin(lat);
    152             double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt)/(1.0 - eslt), Ellipsoid.clarke.e/2.0)
    153                     * l) - Math.PI/2.0;
     198            double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0 - eslt), Ellipsoid.clarke.e / 2.0) * l) - Math.PI
     199                    / 2.0;
    154200            delta = Math.abs(nlt - lat);
    155201            lat = nlt;
    156202        }
    157         return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon));
    158     }
     203        return new LatLon(lat, lon); // in radian
     204    }
     205
     206    /**
     207     * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert
     208     * geographic, (ellipsoid Clark)
     209     *
     210     * @param wgs
     211     * @return
     212     */
     213    private LatLon GRS802Clark(LatLon wgs) {
     214        double lat = Math.toRadians(wgs.lat()); // degree to radian
     215        double lon = Math.toRadians(wgs.lon());
     216        // WGS84 geographic => WGS84 cartesian
     217        double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
     218        double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
     219        double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
     220        double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
     221        // WGS84 => Lambert ellipsoide similarity
     222        X += 168.0;
     223        Y += 60.0;
     224        Z += -320.0;
     225        // Lambert cartesian => Lambert geographic
     226        return Geographic(X, Y, Z, Ellipsoid.clarke);
     227    }
     228
     229    /**
     230     * @param lambert
     231     * @return
     232     */
     233    private LatLon Clark2GRS80(LatLon lambert) {
     234        double lat = lambert.lat(); // in radian
     235        double lon = lambert.lon();
     236        // Lambert geographic => Lambert cartesian
     237        double N = Ellipsoid.clarke.a / (Math.sqrt(1.0 - Ellipsoid.clarke.e2 * Math.sin(lat) * Math.sin(lat)));
     238        double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
     239        double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
     240        double Z = (N * (1.0 - Ellipsoid.clarke.e2)/* + height */) * Math.sin(lat);
     241        // Lambert => WGS84 ellipsoide similarity
     242        X += -168.0;
     243        Y += -60.0;
     244        Z += 320.0;
     245        // WGS84 cartesian => WGS84 geographic
     246        return Geographic(X, Y, Z, Ellipsoid.GRS80);
     247    }
     248
     249    /**
     250     * initializes from cartesian coordinates
     251     *
     252     * @param X
     253     *            1st coordinate in meters
     254     * @param Y
     255     *            2nd coordinate in meters
     256     * @param Z
     257     *            3rd coordinate in meters
     258     * @param ell
     259     *            reference ellipsoid
     260     */
     261    private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
     262        double norm = Math.sqrt(X * X + Y * Y);
     263        double lg = 2.0 * Math.atan(Y / (X + norm));
     264        double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
     265        double delta = 1.0;
     266        while (delta > epsilon) {
     267            double s2 = Math.sin(lt);
     268            s2 *= s2;
     269            double l = Math.atan((Z / norm)
     270                    / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
     271            delta = Math.abs(l - lt);
     272            lt = l;
     273        }
     274        double s2 = Math.sin(lt);
     275        s2 *= s2;
     276        // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
     277        return new LatLon(lt, lg);
     278    }
     279
    159280}
Note: See TracChangeset for help on using the changeset viewer.