- Timestamp:
- 2008-06-25T01:49:12+02:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/org/openstreetmap/josm/data/projection/Lambert.java
r656 r657 5 5 package org.openstreetmap.josm.data.projection; 6 6 7 import static org.openstreetmap.josm.tools.I18n.tr; 8 9 import javax.swing.JOptionPane; 10 11 import org.openstreetmap.josm.Main; 7 12 import org.openstreetmap.josm.data.coor.EastNorth; 8 13 import org.openstreetmap.josm.data.coor.LatLon; … … 12 17 * Lambert I, II, III, and IV projection exponents 13 18 */ 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 17 21 /** 18 22 * Lambert I, II, III, and IV projection constants 19 23 */ 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 23 26 /** 24 27 * Lambert I, II, III, and IV false east 25 28 */ 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 29 31 /** 30 32 * Lambert I, II, III, and IV false north 31 33 */ 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 35 36 /** 36 37 * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian 37 38 */ 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 39 41 /** 40 42 * precision in iterative schema 41 43 */ 44 42 45 public static final double epsilon = 1e-11; 43 46 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 */ 46 73 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 47 79 // 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); 54 116 } 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]); 79 122 } 80 123 81 124 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() { 86 133 return "Lambert"; 87 134 } … … 92 139 93 140 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) { 98 146 return o instanceof Lambert; 99 147 } 100 148 101 @Override public int hashCode() { 149 @Override 150 public int hashCode() { 102 151 return Lambert.class.hashCode(); 103 152 } 104 153 105 154 /** 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) { 121 167 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)); 124 170 double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0)); 125 171 double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0)); 126 172 return new EastNorth(east, north); 127 173 } 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 * 132 179 * @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) { 141 187 double dx = eastNorth.east() - Xs; 142 188 double dy = Ys - eastNorth.north(); 143 double R = Math.sqrt(dx *dx + dy*dy);189 double R = Math.sqrt(dx * dx + dy * dy); 144 190 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)); 146 192 l = Math.exp(l); 147 193 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; 149 195 double delta = 1.0; 150 while (delta > epsilon) {196 while (delta > epsilon) { 151 197 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; 154 200 delta = Math.abs(nlt - lat); 155 201 lat = nlt; 156 202 } 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 159 280 }
Note:
See TracChangeset
for help on using the changeset viewer.