Changeset 1169 in josm for trunk/src/org/openstreetmap/josm/data/projection
- Timestamp:
- 2008-12-23T15:07:05+01:00 (17 years ago)
- Location:
- trunk/src/org/openstreetmap/josm/data/projection
- Files:
-
- 4 edited
-
Epsg4326.java (modified) (2 diffs)
-
Lambert.java (modified) (1 diff)
-
Mercator.java (modified) (2 diffs)
-
Projection.java (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/org/openstreetmap/josm/data/projection/Epsg4326.java
r1032 r1169 12 12 public class Epsg4326 implements Projection { 13 13 14 15 16 14 public EastNorth latlon2eastNorth(LatLon p) { 15 return new EastNorth(p.lon(), p.lat()); 16 } 17 17 18 19 20 18 public LatLon eastNorth2latlon(EastNorth p) { 19 return new LatLon(p.north(), p.east()); 20 } 21 21 22 23 24 22 @Override public String toString() { 23 return tr("EPSG:4326"); 24 } 25 25 26 26 public String getCacheDirectoryName() { … … 28 28 } 29 29 30 31 30 public double scaleFactor() { 31 return 1.0/360; 32 32 } 33 33 34 35 36 34 @Override public boolean equals(Object o) { 35 return o instanceof Epsg4326; 36 } 37 37 38 39 40 38 @Override public int hashCode() { 39 return Epsg4326.class.hashCode(); 40 } 41 41 } -
trunk/src/org/openstreetmap/josm/data/projection/Lambert.java
r865 r1169 14 14 15 15 public class Lambert implements Projection { 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 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 = 1e-11; 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 private static int currentZone = 0; 69 70 private static boolean dontDisplayErrors = false; 71 72 /** 73 * @param p WGS84 lat/lon (ellipsoid GRS80) (in degree) 74 * @return eastnorth projection in Lambert Zone (ellipsoid Clark) 75 */ 76 public EastNorth latlon2eastNorth(LatLon p) { 77 // translate ellipsoid GRS80 (WGS83) => Clark 78 LatLon geo = GRS802Clark(p); 79 double lt = geo.lat(); // in radian 80 double lg = geo.lon(); 81 82 // check if longitude and latitude are inside the french Lambert zones 83 currentZone = 0; 84 boolean outOfLambertZones = false; 85 if (lt >= zoneLimits[3] && lt <= cMaxLatZone1 && lg >= cMinLonZones && lg <= cMaxLonZones) { 86 // zone I 87 if (lt > zoneLimits[0]) 88 currentZone = 0; 89 // zone II 90 else if (lt > zoneLimits[1]) 91 currentZone = 1; 92 // zone III 93 else if (lt > zoneLimits[2]) 94 currentZone = 2; 95 // zone III or IV 96 else if (lt > zoneLimits[3]) 97 // Note: zone IV is dedicated to Corsica island and extends from 47.8 to 98 // 45.9 degrees of latitude. There is an overlap with zone III that can be 99 // solved only with longitude (covers Corsica if lon > 7.2 degree) 100 if (lg < Math.toRadians(8 * 0.9)) 101 currentZone = 2; 102 else 103 currentZone = 3; 104 } else { 105 outOfLambertZones = true; // possible when MAX_LAT is used 106 if (p.lat() != 0 && Math.abs(p.lat()) != Projection.MAX_LAT 107 && p.lon() != 0 && Math.abs(p.lon()) != Projection.MAX_LON 108 && dontDisplayErrors == false) { 109 JOptionPane.showMessageDialog(Main.parent, 110 tr("The projection \"{0}\" is designed for\n" 111 + "latitudes between 46.1° and 57° only.\n" 112 + "Use another projection system if you are not using\n" 113 + "a french WMS server.\n" 114 + "Do not upload any data after this message.", this.toString())); 115 dontDisplayErrors = true; 116 } 117 } 118 if (!outOfLambertZones) { 119 if (layoutZone == -1) { 120 layoutZone = currentZone; 121 dontDisplayErrors = false; 122 } else if (layoutZone != currentZone) { 123 if ((currentZone < layoutZone && Math.abs(zoneLimits[currentZone] - lt) > cMaxOverlappingZones) 124 || (currentZone > layoutZone && Math.abs(zoneLimits[layoutZone] - lt) > cMaxOverlappingZones)) { 125 JOptionPane.showMessageDialog(Main.parent, 126 tr("IMPORTANT : data positionned far away from\n" 127 + "the current Lambert zone limits.\n" 128 + "Do not upload any data after this message.\n" 129 + "Undo your last action, Save your work \n" 130 + "and Start a new layer on the new zone.")); 131 layoutZone = -1; 132 dontDisplayErrors = true; 133 } else { 134 System.out.println("temporarily extend Lambert zone " + layoutZone + " projection at lat,lon:" 135 + lt + "," + lg); 136 } 137 } 138 } 139 if (layoutZone == -1) { 140 return ConicProjection(lt, lg, Xs[currentZone], Ys[currentZone], c[currentZone], n[currentZone]); 141 } // else 142 return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]); 143 } 144 145 public LatLon eastNorth2latlon(EastNorth p) { 146 LatLon geo; 147 if (layoutZone == -1) 148 // possible until the Lambert zone is determined by latlon2eastNorth() with a valid LatLon 149 geo = Geographic(p, Xs[currentZone], Ys[currentZone], c[currentZone], n[currentZone]); 150 else 151 geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]); 152 // translate ellipsoid Clark => GRS80 (WGS83) 153 LatLon wgs = Clark2GRS80(geo); 154 return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon())); 155 } 156 157 @Override public String toString() { 158 return tr("Lambert Zone (France)"); 159 } 160 161 public String getCacheDirectoryName() { 162 return "lambert"; 163 } 164 165 public double scaleFactor() { 166 return 1.0 / 360; 167 } 168 169 @Override 170 public boolean equals(Object o) { 171 return o instanceof Lambert; 172 } 173 174 @Override 175 public int hashCode() { 176 return Lambert.class.hashCode(); 177 } 178 179 /** 180 * Initializes from geographic coordinates. Note that reference ellipsoid 181 * used by Lambert is the Clark ellipsoid. 182 * 183 * @param lat latitude in grad 184 * @param lon longitude in grad 185 * @param Xs false east (coordinate system origin) in meters 186 * @param Ys false north (coordinate system origin) in meters 187 * @param c projection constant 188 * @param n projection exponent 189 * @return EastNorth projected coordinates in meter 190 */ 191 private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) { 192 double eslt = Ellipsoid.clarke.e * Math.sin(lat); 193 double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0)) 194 * Math.pow((1.0 - eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0)); 195 double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0)); 196 double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0)); 197 return new EastNorth(east, north); 198 } 199 200 /** 201 * Initializes from projected coordinates (conic projection). Note that 202 * reference ellipsoid used by Lambert is Clark 203 * 204 * @param coord projected coordinates pair in meters 205 * @param Xs false east (coordinate system origin) in meters 206 * @param Ys false north (coordinate system origin) in meters 207 * @param c projection constant 208 * @param n projection exponent 209 * @return LatLon in radian 210 */ 211 private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) { 212 double dx = eastNorth.east() - Xs; 213 double dy = Ys - eastNorth.north(); 214 double R = Math.sqrt(dx * dx + dy * dy); 215 double gamma = Math.atan(dx / dy); 216 double l = -1.0 / n * Math.log(Math.abs(R / c)); 217 l = Math.exp(l); 218 double lon = lg0 + gamma / n; 219 double lat = 2.0 * Math.atan(l) - Math.PI / 2.0; 220 double delta = 1.0; 221 while (delta > epsilon) { 222 double eslt = Ellipsoid.clarke.e * Math.sin(lat); 223 double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0 - eslt), Ellipsoid.clarke.e / 2.0) * l) - Math.PI 224 / 2.0; 225 delta = Math.abs(nlt - lat); 226 lat = nlt; 227 } 228 return new LatLon(lat, lon); // in radian 229 } 230 231 /** 232 * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert 233 * geographic, (ellipsoid Clark) 234 * 235 * @param wgs 236 * @return 237 */ 238 private LatLon GRS802Clark(LatLon wgs) { 239 double lat = Math.toRadians(wgs.lat()); // degree to radian 240 double lon = Math.toRadians(wgs.lon()); 241 // WGS84 geographic => WGS84 cartesian 242 double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat))); 243 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon); 244 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon); 245 double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat); 246 // WGS84 => Lambert ellipsoide similarity 247 X += 168.0; 248 Y += 60.0; 249 Z += -320.0; 250 // Lambert cartesian => Lambert geographic 251 return Geographic(X, Y, Z, Ellipsoid.clarke); 252 } 253 254 /** 255 * @param lambert 256 * @return 257 */ 258 private LatLon Clark2GRS80(LatLon lambert) { 259 double lat = lambert.lat(); // in radian 260 double lon = lambert.lon(); 261 // Lambert geographic => Lambert cartesian 262 double N = Ellipsoid.clarke.a / (Math.sqrt(1.0 - Ellipsoid.clarke.e2 * Math.sin(lat) * Math.sin(lat))); 263 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon); 264 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon); 265 double Z = (N * (1.0 - Ellipsoid.clarke.e2)/* + height */) * Math.sin(lat); 266 // Lambert => WGS84 ellipsoide similarity 267 X += -168.0; 268 Y += -60.0; 269 Z += 320.0; 270 // WGS84 cartesian => WGS84 geographic 271 return Geographic(X, Y, Z, Ellipsoid.GRS80); 272 } 273 274 /** 275 * initializes from cartesian coordinates 276 * 277 * @param X 278 * 1st coordinate in meters 279 * @param Y 280 * 2nd coordinate in meters 281 * @param Z 282 * 3rd coordinate in meters 283 * @param ell 284 * reference ellipsoid 285 */ 286 private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) { 287 double norm = Math.sqrt(X * X + Y * Y); 288 double lg = 2.0 * Math.atan(Y / (X + norm)); 289 double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z))))); 290 double delta = 1.0; 291 while (delta > epsilon) { 292 double s2 = Math.sin(lt); 293 s2 *= s2; 294 double l = Math.atan((Z / norm) 295 / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2))))); 296 delta = Math.abs(l - lt); 297 lt = l; 298 } 299 double s2 = Math.sin(lt); 300 s2 *= s2; 301 // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2); 302 return new LatLon(lt, lg); 303 } 304 304 305 305 } -
trunk/src/org/openstreetmap/josm/data/projection/Mercator.java
r1032 r1169 10 10 * Implement Mercator Projection code, coded after documentation 11 11 * from wikipedia. 12 * 13 * The center of the mercator projection is always the 0 grad 12 * 13 * The center of the mercator projection is always the 0 grad 14 14 * coordinate. 15 * 15 * 16 16 * @author imi 17 17 */ 18 18 public class Mercator implements Projection { 19 19 20 21 22 23 24 20 public EastNorth latlon2eastNorth(LatLon p) { 21 return new EastNorth( 22 p.lon()*Math.PI/180, 23 Math.log(Math.tan(Math.PI/4+p.lat()*Math.PI/360))); 24 } 25 25 26 27 28 29 30 26 public LatLon eastNorth2latlon(EastNorth p) { 27 return new LatLon( 28 Math.atan(Math.sinh(p.north()))*180/Math.PI, 29 p.east()*180/Math.PI); 30 } 31 31 32 33 34 32 @Override public String toString() { 33 return tr("Mercator"); 34 } 35 35 36 36 public String getCacheDirectoryName() { … … 38 38 } 39 39 40 41 40 public double scaleFactor() { 41 return 1/Math.PI/2; 42 42 } 43 43 44 45 46 44 @Override public boolean equals(Object o) { 45 return o instanceof Mercator; 46 } 47 47 48 49 50 48 @Override public int hashCode() { 49 return Mercator.class.hashCode(); 50 } 51 51 } -
trunk/src/org/openstreetmap/josm/data/projection/Projection.java
r656 r1169 6 6 7 7 /** 8 * Classes implementing this are able to convert lat/lon values to 8 * Classes implementing this are able to convert lat/lon values to 9 9 * planar screen coordinates. 10 * 10 * 11 11 * @author imi 12 12 */ 13 13 public interface Projection { 14 14 15 /** 16 * Maximum latitude representable. 17 */ 18 public static final double MAX_LAT = 85.05112877980659; // Mercator squares the world 19 20 /** 21 * Maximum longditude representable. 22 */ 23 public static final double MAX_LON = 180; 15 /** 16 * Maximum latitude representable. 17 */ 18 public static final double MAX_LAT = 85.05112877980659; // Mercator squares the world 24 19 25 26 * Minimum difference in location to not be represented as the same position.27 28 SERVER_PRECISION = 1e12;20 /** 21 * Maximum longditude representable. 22 */ 23 public static final double MAX_LON = 180; 29 24 30 /** 31 * List of all available projections. 32 */ 33 public static Projection[] allProjections = new Projection[]{ 34 new Epsg4326(), 35 new Mercator(), 36 new Lambert() 37 }; 38 39 /** 40 * Convert from lat/lon to northing/easting. 41 * 42 * @param p The geo point to convert. x/y members of the point are filled. 43 */ 44 EastNorth latlon2eastNorth(LatLon p); 45 46 /** 47 * Convert from norting/easting to lat/lon. 48 * 49 * @param p The geo point to convert. lat/lon members of the point are filled. 50 */ 51 LatLon eastNorth2latlon(EastNorth p); 25 /** 26 * Minimum difference in location to not be represented as the same position. 27 */ 28 public static final double MAX_SERVER_PRECISION = 1e12; 52 29 53 /** 54 * Describe the projection converter in one or two words. 55 */ 56 String toString(); 57 30 /** 31 * List of all available projections. 32 */ 33 public static Projection[] allProjections = new Projection[]{ 34 new Epsg4326(), 35 new Mercator(), 36 new Lambert() 37 }; 38 39 /** 40 * Convert from lat/lon to northing/easting. 41 * 42 * @param p The geo point to convert. x/y members of the point are filled. 43 */ 44 EastNorth latlon2eastNorth(LatLon p); 45 46 /** 47 * Convert from norting/easting to lat/lon. 48 * 49 * @param p The geo point to convert. lat/lon members of the point are filled. 50 */ 51 LatLon eastNorth2latlon(EastNorth p); 52 53 /** 54 * Describe the projection converter in one or two words. 55 */ 56 String toString(); 57 58 58 /** 59 59 * Get a filename compatible string (for the cache directory) 60 60 */ 61 61 String getCacheDirectoryName(); 62 62 63 63 /** 64 * The factor to multiply with an easting coordinate to get from "easting 64 * The factor to multiply with an easting coordinate to get from "easting 65 65 * units per pixel" to "meters per pixel" 66 66 */
Note:
See TracChangeset
for help on using the changeset viewer.
