Index: trunk/src/org/openstreetmap/josm/data/projection/Lambert.java =================================================================== --- trunk/src/org/openstreetmap/josm/data/projection/Lambert.java (revision 656) +++ trunk/src/org/openstreetmap/josm/data/projection/Lambert.java (revision 657) @@ -5,4 +5,9 @@ package org.openstreetmap.josm.data.projection; +import static org.openstreetmap.josm.tools.I18n.tr; + +import javax.swing.JOptionPane; + +import org.openstreetmap.josm.Main; import org.openstreetmap.josm.data.coor.EastNorth; import org.openstreetmap.josm.data.coor.LatLon; @@ -12,76 +17,118 @@ * Lambert I, II, III, and IV projection exponents */ - public static final double n[] = { - 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 - }; + public static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 }; + /** * Lambert I, II, III, and IV projection constants */ - public static final double c[] = { - 11603796.98, 11745793.39, 11947992.52, 12136281.99 - }; + public static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 }; + /** * Lambert I, II, III, and IV false east */ - public static final double Xs[] = { - 600000.0, 600000.0, 600000.0, 234.358 - }; + public static final double Xs[] = { 600000.0, 600000.0, 600000.0, 234.358 }; + /** * Lambert I, II, III, and IV false north */ - public static final double Ys[] = { - 5657616.674, 6199695.768, 6791905.085, 7239161.542 - }; + public static final double Ys[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 }; + /** * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian */ - public static final double lg0 = 0.04079234433198; // 2 deg 20 min 14.025 sec + public static final double lg0 = 0.04079234433198; // 2deg20'14.025" + /** * precision in iterative schema */ + public static final double epsilon = 1e-11; - public static int zone = 0; - + /** + * France is divided in 4 Lambert projection zones (1,2,3 + 4th for Corsica) + */ + public static final double cMaxLatZone1 = Math.toRadians(57 * 0.9); + + public static final double zoneLimits[] = { Math.toRadians(53.5 * 0.9), // between Zone 1 and Zone 2 (in grad *0.9) + Math.toRadians(50.5 * 0.9), // between Zone 2 and Zone 3 + Math.toRadians(47.51963 * 0.9), // between Zone 3 and Zone 4 + Math.toRadians(46.17821 * 0.9) };// lowest latitude of Zone 4 + + public static final double cMinLonZones = Math.toRadians(-4.9074074074074059 * 0.9); + + public static final double cMaxLonZones = Math.toRadians(10.2 * 0.9); + + /** + * Because josm cannot work correctly if two zones are displayed, we allow some overlapping + */ + public static final double cMaxOverlappingZones = Math.toRadians(0.5 * 0.9); + + public static int layoutZone = -1; + + /** + * @param p + * a WGS84 lat/lon (ellipsoid GRS80) (in degree) + * @return eastnorth projection in Lambert Zone (ellipsoid Clark) + */ public EastNorth latlon2eastNorth(LatLon p) { + // translate ellipsoid GRS80 (WGS83) => Clark + LatLon geo = GRS802Clark(p); + double lt = geo.lat(); // in radian + double lg = geo.lon(); + // check if longitude and latitude are inside the french Lambert zones - double lt = p.lat(); - double lg = p.lon(); - if(lt > 57*9/10 || lt < 46.17821*9/10 || - lg > 10.2*9/10 || lg < -4.9074074074074059*9/10) { - if (lt > 57*9/10) { - // out of Lambert zones, possible when MAX_LAT is used- keep current zone + int currentZone = 0; + boolean outOfLambertZones = false; + if (lt >= zoneLimits[3] && lt <= cMaxLatZone1 && lg >= cMinLonZones && lg <= cMaxLonZones) { + // zone I + if (lt > zoneLimits[0]) + currentZone = 0; + // zone II + else if (lt > zoneLimits[1]) + currentZone = 1; + // zone III + else if (lt > zoneLimits[2]) + currentZone = 2; + // zone III or IV + else if (lt > zoneLimits[3]) + // Note: zone IV is dedicated to Corsica island and extends from 47.8 to + // 45.9 degrees of latitude. There is an overlap with zone III that can be + // solved only with longitude (covers Corsica if lon > 7.2 degree) + if (lg < Math.toRadians(8 * 0.9)) + currentZone = 2; + else + currentZone = 3; + } else { + outOfLambertZones = true; // possible when MAX_LAT is used + } + if (layoutZone == -1) + layoutZone = currentZone; + else if (layoutZone != currentZone && !outOfLambertZones) { + if ((currentZone < layoutZone && Math.abs(zoneLimits[currentZone] - lt) > cMaxOverlappingZones) + || (currentZone > layoutZone && Math.abs(zoneLimits[layoutZone] - lt) > cMaxOverlappingZones)) { + JOptionPane.showMessageDialog(Main.parent, + tr("Some data are positionned far away from current Lambert zone limits.\n" + +"Split long ways to avoid distortions.")); + layoutZone = -1; + } else { + System.out.println("temporarily extends Lambert zone " + layoutZone + + " projection at lat,lon:" + lt + "," + lg); } - // zone I - else if (lt > 53.5*9/10) - zone = 0; - // zone II - else if(lt > 50.5*9/10) - zone = 1; - // zone III - else if(lt > 47.51963*9/10) - zone = 2; - else if (lt > 46.17821*9/10) - // zone III or IV - // Note: zone IV is dedicated to Corsica island and extends - // from 47.8 to 45.9 degrees of latitude. There is an overlap with - // zone III that can be solved only with longitude - // (covers Corsica if lon > 7.2 degree) - if (lg < 8*9/10) - zone = 2; - else - zone = 3; - } else { - // out of Lambert zones- keep current one - } - EastNorth eastNorth = ConicProjection(lt, lg, Xs[zone], Ys[zone], c[zone], n[zone]); - return eastNorth; + } + if (layoutZone == -1) { + return ConicProjection(lt, lg, Xs[currentZone], Ys[currentZone], c[currentZone], n[currentZone]); + } // else + return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]); } public LatLon eastNorth2latlon(EastNorth p) { - return Geographic(p, Xs[zone], Ys[zone], c[zone], n[zone]); - } - - @Override public String toString() { + LatLon geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]); + // translate ellipsoid Clark => GRS80 (WGS83) + LatLon wgs = Clark2GRS80(geo); + return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon())); + } + + @Override + public String toString() { return "Lambert"; } @@ -92,68 +139,142 @@ public double scaleFactor() { - return 1.0/360; - } - - @Override public boolean equals(Object o) { + return 1.0 / 360; + } + + @Override + public boolean equals(Object o) { return o instanceof Lambert; } - @Override public int hashCode() { + @Override + public int hashCode() { return Lambert.class.hashCode(); } /** - * Initializes from geographic coordinates. - * Note that reference ellipsoid used by Lambert is the Clark ellipsoid. - * - * @param lat latitude in degree - * @param lon longitude in degree - * @param Xs false east (coordinate system origin) in meters - * @param Ys false north (coordinate system origin) in meters - * @param c projection constant - * @param n projection exponent - * @return EastNorth projected coordinates in meter - */ - public EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, - double c, double n) { - lat = Math.toRadians(lat); - lon = Math.toRadians(lon); + * Initializes from geographic coordinates. Note that reference ellipsoid + * used by Lambert is the Clark ellipsoid. + * + * @param lat latitude in grad + * @param lon longitude in grad + * @param Xs false east (coordinate system origin) in meters + * @param Ys false north (coordinate system origin) in meters + * @param c projection constant + * @param n projection exponent + * @return EastNorth projected coordinates in meter + */ + private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) { double eslt = Ellipsoid.clarke.e * Math.sin(lat); - double l = Math.log(Math.tan(Math.PI/4.0 + (lat/2.0)) * - Math.pow((1.0 - eslt)/(1.0 + eslt), Ellipsoid.clarke.e/2.0)); + double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0)) + * Math.pow((1.0 - eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0)); double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0)); double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0)); return new EastNorth(east, north); } - /** - * Initializes from projected coordinates (conic projection). - * Note that reference ellipsoid used by Lambert is Clark - * + + /** + * Initializes from projected coordinates (conic projection). Note that + * reference ellipsoid used by Lambert is Clark + * * @param coord projected coordinates pair in meters - * @param Xs false east (coordinate system origin) in meters - * @param Ys false north (coordinate system origin) in meters - * @param c projection constant - * @param n projection exponent - * @return LatLon in degrees - */ - public LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, - double c, double n) { + * @param Xs false east (coordinate system origin) in meters + * @param Ys false north (coordinate system origin) in meters + * @param c projection constant + * @param n projection exponent + * @return LatLon in radian + */ + private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) { double dx = eastNorth.east() - Xs; double dy = Ys - eastNorth.north(); - double R = Math.sqrt(dx*dx + dy*dy); + double R = Math.sqrt(dx * dx + dy * dy); double gamma = Math.atan(dx / dy); - double l = -1.0/n * Math.log(Math.abs(R / c)); + double l = -1.0 / n * Math.log(Math.abs(R / c)); l = Math.exp(l); double lon = lg0 + gamma / n; - double lat = 2.0 * Math.atan(l) - Math.PI/2.0; + double lat = 2.0 * Math.atan(l) - Math.PI / 2.0; double delta = 1.0; - while(delta > epsilon) { + while (delta > epsilon) { double eslt = Ellipsoid.clarke.e * Math.sin(lat); - double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt)/(1.0 - eslt), Ellipsoid.clarke.e/2.0) - * l) - Math.PI/2.0; + double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0 - eslt), Ellipsoid.clarke.e / 2.0) * l) - Math.PI + / 2.0; delta = Math.abs(nlt - lat); lat = nlt; } - return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon)); - } + return new LatLon(lat, lon); // in radian + } + + /** + * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert + * geographic, (ellipsoid Clark) + * + * @param wgs + * @return + */ + private LatLon GRS802Clark(LatLon wgs) { + double lat = Math.toRadians(wgs.lat()); // degree to radian + double lon = Math.toRadians(wgs.lon()); + // WGS84 geographic => WGS84 cartesian + double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat))); + double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon); + double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon); + double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat); + // WGS84 => Lambert ellipsoide similarity + X += 168.0; + Y += 60.0; + Z += -320.0; + // Lambert cartesian => Lambert geographic + return Geographic(X, Y, Z, Ellipsoid.clarke); + } + + /** + * @param lambert + * @return + */ + private LatLon Clark2GRS80(LatLon lambert) { + double lat = lambert.lat(); // in radian + double lon = lambert.lon(); + // Lambert geographic => Lambert cartesian + double N = Ellipsoid.clarke.a / (Math.sqrt(1.0 - Ellipsoid.clarke.e2 * Math.sin(lat) * Math.sin(lat))); + double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon); + double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon); + double Z = (N * (1.0 - Ellipsoid.clarke.e2)/* + height */) * Math.sin(lat); + // Lambert => WGS84 ellipsoide similarity + X += -168.0; + Y += -60.0; + Z += 320.0; + // WGS84 cartesian => WGS84 geographic + return Geographic(X, Y, Z, Ellipsoid.GRS80); + } + + /** + * initializes from cartesian coordinates + * + * @param X + * 1st coordinate in meters + * @param Y + * 2nd coordinate in meters + * @param Z + * 3rd coordinate in meters + * @param ell + * reference ellipsoid + */ + private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) { + double norm = Math.sqrt(X * X + Y * Y); + double lg = 2.0 * Math.atan(Y / (X + norm)); + double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z))))); + double delta = 1.0; + while (delta > epsilon) { + double s2 = Math.sin(lt); + s2 *= s2; + double l = Math.atan((Z / norm) + / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2))))); + delta = Math.abs(l - lt); + lt = l; + } + double s2 = Math.sin(lt); + s2 *= s2; + // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2); + return new LatLon(lt, lg); + } + }