Changeset 810 in josm for trunk/src/org/openstreetmap/josm/data/projection
 Timestamp:
 20080820T17:14:17+02:00 (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/org/openstreetmap/josm/data/projection/Lambert.java
r687 r810 1 1 //License: GPL. For details, see LICENSE file. 2 2 //Thanks to Johan Montagnat and its geoconv java converter application 3 //(http://www.i3s.unice.fr/~johan/gps/ , published under GPL license) 3 //(http://www.i3s.unice.fr/~johan/gps/ , published under GPL license) 4 4 //from which some code and constants have been reused here. 5 5 package org.openstreetmap.josm.data.projection; … … 14 14 15 15 public class Lambert implements Projection { 16 /** 17 * Lambert I, II, III, and IV projection exponents 18 */ 19 public static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 }; 20 21 /** 22 * Lambert I, II, III, and IV projection constants 23 */ 24 public static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 }; 25 26 /** 27 * Lambert I, II, III, and IV false east 28 */ 29 public static final double Xs[] = { 600000.0, 600000.0, 600000.0, 234.358 }; 30 31 /** 32 * Lambert I, II, III, and IV false north 33 */ 34 public static final double Ys[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 }; 35 36 /** 37 * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian 38 */ 39 public static final double lg0 = 0.04079234433198; // 2deg20'14.025" 40 41 /** 42 * precision in iterative schema 43 */ 44 45 public static final double epsilon = 1e11; 46 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(1.5 * 0.9); 65 66 public static int layoutZone = 1; 67 68 /** 69 * @param p WGS84 lat/lon (ellipsoid GRS80) (in degree) 70 * @return eastnorth projection in Lambert Zone (ellipsoid Clark) 71 */ 72 public EastNorth latlon2eastNorth(LatLon p) { 73 // translate ellipsoid GRS80 (WGS83) => Clark 74 LatLon geo = GRS802Clark(p); 75 double lt = geo.lat(); // in radian 76 double lg = geo.lon(); 77 78 // check if longitude and latitude are inside the french Lambert zones 79 int currentZone = 0; 80 boolean outOfLambertZones = false; 81 if (lt >= zoneLimits[3] && lt <= cMaxLatZone1 && lg >= cMinLonZones && lg <= cMaxLonZones) { 82 // zone I 83 if (lt > zoneLimits[0]) 84 currentZone = 0; 85 // zone II 86 else if (lt > zoneLimits[1]) 87 currentZone = 1; 88 // zone III 89 else if (lt > zoneLimits[2]) 90 currentZone = 2; 91 // zone III or IV 92 else if (lt > zoneLimits[3]) 93 // Note: zone IV is dedicated to Corsica island and extends from 47.8 to 94 // 45.9 degrees of latitude. There is an overlap with zone III that can be 95 // solved only with longitude (covers Corsica if lon > 7.2 degree) 96 if (lg < Math.toRadians(8 * 0.9)) 97 currentZone = 2; 98 else 99 currentZone = 3; 100 } else { 101 outOfLambertZones = true; // possible when MAX_LAT is used 102 } 103 if (layoutZone == 1) 104 layoutZone = currentZone; 105 else if (layoutZone != currentZone && !outOfLambertZones) { 106 if ((currentZone < layoutZone && Math.abs(zoneLimits[currentZone]  lt) > cMaxOverlappingZones) 107  (currentZone > layoutZone && Math.abs(zoneLimits[layoutZone]  lt) > cMaxOverlappingZones)) { 108 JOptionPane.showMessageDialog(Main.parent, 109 tr("IMPORTANT : data positionned far away from\n" 110 +"the current Lambert zone limits.\n" 111 +"Undo your last action, Save your work \n" 112 +"and Start a new layer on the new zone.")); 113 layoutZone = 1; 114 } else { 115 System.out.println("temporarily extends Lambert zone " + layoutZone 116 + " projection at lat,lon:" + lt + "," + lg); 117 } 118 } 119 if (layoutZone == 1) { 120 return ConicProjection(lt, lg, Xs[currentZone], Ys[currentZone], c[currentZone], n[currentZone]); 121 } // else 122 return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]); 123 } 124 125 public LatLon eastNorth2latlon(EastNorth p) { 126 LatLon geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]); 127 // translate ellipsoid Clark => GRS80 (WGS83) 128 LatLon wgs = Clark2GRS80(geo); 129 return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon())); 130 } 131 132 @Override 133 public String toString() { 134 return "Lambert"; 135 } 136 137 public String getCacheDirectoryName() { 138 return "lambert"; 139 } 140 141 public double scaleFactor() { 142 return 1.0 / 360; 143 } 144 145 @Override 146 public boolean equals(Object o) { 147 return o instanceof Lambert; 148 } 149 150 @Override 151 public int hashCode() { 152 return Lambert.class.hashCode(); 153 } 154 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 coord 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 radian 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(lat, lon); // in radian 205 } 206 207 /** 208 * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert 209 * geographic, (ellipsoid Clark) 210 * 211 * @param wgs 212 * @return 213 */ 214 private LatLon GRS802Clark(LatLon wgs) { 215 double lat = Math.toRadians(wgs.lat()); // degree to radian 216 double lon = Math.toRadians(wgs.lon()); 217 // WGS84 geographic => WGS84 cartesian 218 double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0  Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat))); 219 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon); 220 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon); 221 double Z = (N * (1.0  Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat); 222 // WGS84 => Lambert ellipsoide similarity 223 X += 168.0; 224 Y += 60.0; 225 Z += 320.0; 226 // Lambert cartesian => Lambert geographic 227 return Geographic(X, Y, Z, Ellipsoid.clarke); 228 } 229 230 /** 231 * @param lambert 232 * @return 233 */ 234 private LatLon Clark2GRS80(LatLon lambert) { 235 double lat = lambert.lat(); // in radian 236 double lon = lambert.lon(); 237 // Lambert geographic => Lambert cartesian 238 double N = Ellipsoid.clarke.a / (Math.sqrt(1.0  Ellipsoid.clarke.e2 * Math.sin(lat) * Math.sin(lat))); 239 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon); 240 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon); 241 double Z = (N * (1.0  Ellipsoid.clarke.e2)/* + height */) * Math.sin(lat); 242 // Lambert => WGS84 ellipsoide similarity 243 X += 168.0; 244 Y += 60.0; 245 Z += 320.0; 246 // WGS84 cartesian => WGS84 geographic 247 return Geographic(X, Y, Z, Ellipsoid.GRS80); 248 } 249 250 /** 251 * initializes from cartesian coordinates 252 * 253 * @param X 254 * 1st coordinate in meters 255 * @param Y 256 * 2nd coordinate in meters 257 * @param Z 258 * 3rd coordinate in meters 259 * @param ell 260 * reference ellipsoid 261 */ 262 private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) { 263 double norm = Math.sqrt(X * X + Y * Y); 264 double lg = 2.0 * Math.atan(Y / (X + norm)); 265 double lt = Math.atan(Z / (norm * (1.0  (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z))))); 266 double delta = 1.0; 267 while (delta > epsilon) { 268 double s2 = Math.sin(lt); 269 s2 *= s2; 270 double l = Math.atan((Z / norm) 271 / (1.0  (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0  ell.e2 * s2))))); 272 delta = Math.abs(l  lt); 273 lt = l; 274 } 275 double s2 = Math.sin(lt); 276 s2 *= s2; 277 // h = norm / Math.cos(lt)  ell.a / Math.sqrt(1.0  ell.e2 * s2); 278 return new LatLon(lt, lg); 279 } 16 /** 17 * Lambert I, II, III, and IV projection exponents 18 */ 19 public static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 }; 20 21 /** 22 * Lambert I, II, III, and IV projection constants 23 */ 24 public static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 }; 25 26 /** 27 * Lambert I, II, III, and IV false east 28 */ 29 public static final double Xs[] = { 600000.0, 600000.0, 600000.0, 234.358 }; 30 31 /** 32 * Lambert I, II, III, and IV false north 33 */ 34 public static final double Ys[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 }; 35 36 /** 37 * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian 38 */ 39 public static final double lg0 = 0.04079234433198; // 2deg20'14.025" 40 41 /** 42 * precision in iterative schema 43 */ 44 45 public static final double epsilon = 1e11; 46 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(1.5 * 0.9); 65 66 public static int layoutZone = 1; 67 68 /** 69 * @param p WGS84 lat/lon (ellipsoid GRS80) (in degree) 70 * @return eastnorth projection in Lambert Zone (ellipsoid Clark) 71 */ 72 public EastNorth latlon2eastNorth(LatLon p) { 73 // translate ellipsoid GRS80 (WGS83) => Clark 74 LatLon geo = GRS802Clark(p); 75 double lt = geo.lat(); // in radian 76 double lg = geo.lon(); 77 78 // check if longitude and latitude are inside the french Lambert zones 79 int currentZone = 0; 80 boolean outOfLambertZones = false; 81 if (lt >= zoneLimits[3] && lt <= cMaxLatZone1 && lg >= cMinLonZones && lg <= cMaxLonZones) { 82 // zone I 83 if (lt > zoneLimits[0]) 84 currentZone = 0; 85 // zone II 86 else if (lt > zoneLimits[1]) 87 currentZone = 1; 88 // zone III 89 else if (lt > zoneLimits[2]) 90 currentZone = 2; 91 // zone III or IV 92 else if (lt > zoneLimits[3]) 93 // Note: zone IV is dedicated to Corsica island and extends from 47.8 to 94 // 45.9 degrees of latitude. There is an overlap with zone III that can be 95 // solved only with longitude (covers Corsica if lon > 7.2 degree) 96 if (lg < Math.toRadians(8 * 0.9)) 97 currentZone = 2; 98 else 99 currentZone = 3; 100 } else { 101 outOfLambertZones = true; // possible when MAX_LAT is used 102 } 103 if (!outOfLambertZones) { 104 if (layoutZone == 1) 105 layoutZone = currentZone; 106 else if (layoutZone != currentZone) { 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("IMPORTANT : data positionned far away from\n" 111 + "the current Lambert zone limits.\n" 112 + "Undo your last action, Save your work \n" 113 + "and Start a new layer on the new zone.")); 114 layoutZone = 1; 115 } else { 116 System.out.println("temporarily extends Lambert zone " 117 + layoutZone + " projection at lat,lon:" + lt + "," 118 + lg); 119 } 120 } 121 } 122 if (layoutZone == 1) { 123 return ConicProjection(lt, lg, Xs[currentZone], Ys[currentZone], c[currentZone], n[currentZone]); 124 } // else 125 return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]); 126 } 127 128 public LatLon eastNorth2latlon(EastNorth p) { 129 LatLon geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]); 130 // translate ellipsoid Clark => GRS80 (WGS83) 131 LatLon wgs = Clark2GRS80(geo); 132 return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon())); 133 } 134 135 @Override 136 public String toString() { 137 return "Lambert"; 138 } 139 140 public String getCacheDirectoryName() { 141 return "lambert"; 142 } 143 144 public double scaleFactor() { 145 return 1.0 / 360; 146 } 147 148 @Override 149 public boolean equals(Object o) { 150 return o instanceof Lambert; 151 } 152 153 @Override 154 public int hashCode() { 155 return Lambert.class.hashCode(); 156 } 157 158 /** 159 * Initializes from geographic coordinates. Note that reference ellipsoid 160 * used by Lambert is the Clark ellipsoid. 161 * 162 * @param lat latitude in grad 163 * @param lon longitude in grad 164 * @param Xs false east (coordinate system origin) in meters 165 * @param Ys false north (coordinate system origin) in meters 166 * @param c projection constant 167 * @param n projection exponent 168 * @return EastNorth projected coordinates in meter 169 */ 170 private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) { 171 double eslt = Ellipsoid.clarke.e * Math.sin(lat); 172 double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0)) 173 * Math.pow((1.0  eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0)); 174 double east = Xs + c * Math.exp(n * l) * Math.sin(n * (lon  lg0)); 175 double north = Ys  c * Math.exp(n * l) * Math.cos(n * (lon  lg0)); 176 return new EastNorth(east, north); 177 } 178 179 /** 180 * Initializes from projected coordinates (conic projection). Note that 181 * reference ellipsoid used by Lambert is Clark 182 * 183 * @param coord projected coordinates pair in meters 184 * @param Xs false east (coordinate system origin) in meters 185 * @param Ys false north (coordinate system origin) in meters 186 * @param c projection constant 187 * @param n projection exponent 188 * @return LatLon in radian 189 */ 190 private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) { 191 double dx = eastNorth.east()  Xs; 192 double dy = Ys  eastNorth.north(); 193 double R = Math.sqrt(dx * dx + dy * dy); 194 double gamma = Math.atan(dx / dy); 195 double l = 1.0 / n * Math.log(Math.abs(R / c)); 196 l = Math.exp(l); 197 double lon = lg0 + gamma / n; 198 double lat = 2.0 * Math.atan(l)  Math.PI / 2.0; 199 double delta = 1.0; 200 while (delta > epsilon) { 201 double eslt = Ellipsoid.clarke.e * Math.sin(lat); 202 double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0  eslt), Ellipsoid.clarke.e / 2.0) * l)  Math.PI 203 / 2.0; 204 delta = Math.abs(nlt  lat); 205 lat = nlt; 206 } 207 return new LatLon(lat, lon); // in radian 208 } 209 210 /** 211 * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert 212 * geographic, (ellipsoid Clark) 213 * 214 * @param wgs 215 * @return 216 */ 217 private LatLon GRS802Clark(LatLon wgs) { 218 double lat = Math.toRadians(wgs.lat()); // degree to radian 219 double lon = Math.toRadians(wgs.lon()); 220 // WGS84 geographic => WGS84 cartesian 221 double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0  Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat))); 222 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon); 223 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon); 224 double Z = (N * (1.0  Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat); 225 // WGS84 => Lambert ellipsoide similarity 226 X += 168.0; 227 Y += 60.0; 228 Z += 320.0; 229 // Lambert cartesian => Lambert geographic 230 return Geographic(X, Y, Z, Ellipsoid.clarke); 231 } 232 233 /** 234 * @param lambert 235 * @return 236 */ 237 private LatLon Clark2GRS80(LatLon lambert) { 238 double lat = lambert.lat(); // in radian 239 double lon = lambert.lon(); 240 // Lambert geographic => Lambert cartesian 241 double N = Ellipsoid.clarke.a / (Math.sqrt(1.0  Ellipsoid.clarke.e2 * Math.sin(lat) * Math.sin(lat))); 242 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon); 243 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon); 244 double Z = (N * (1.0  Ellipsoid.clarke.e2)/* + height */) * Math.sin(lat); 245 // Lambert => WGS84 ellipsoide similarity 246 X += 168.0; 247 Y += 60.0; 248 Z += 320.0; 249 // WGS84 cartesian => WGS84 geographic 250 return Geographic(X, Y, Z, Ellipsoid.GRS80); 251 } 252 253 /** 254 * initializes from cartesian coordinates 255 * 256 * @param X 257 * 1st coordinate in meters 258 * @param Y 259 * 2nd coordinate in meters 260 * @param Z 261 * 3rd coordinate in meters 262 * @param ell 263 * reference ellipsoid 264 */ 265 private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) { 266 double norm = Math.sqrt(X * X + Y * Y); 267 double lg = 2.0 * Math.atan(Y / (X + norm)); 268 double lt = Math.atan(Z / (norm * (1.0  (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z))))); 269 double delta = 1.0; 270 while (delta > epsilon) { 271 double s2 = Math.sin(lt); 272 s2 *= s2; 273 double l = Math.atan((Z / norm) 274 / (1.0  (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0  ell.e2 * s2))))); 275 delta = Math.abs(l  lt); 276 lt = l; 277 } 278 double s2 = Math.sin(lt); 279 s2 *= s2; 280 // h = norm / Math.cos(lt)  ell.a / Math.sqrt(1.0  ell.e2 * s2); 281 return new LatLon(lt, lg); 282 } 280 283 281 284 }
Note: See TracChangeset
for help on using the changeset viewer.