001// License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.data.projection.proj; 003 004import static org.openstreetmap.josm.tools.I18n.tr; 005 006import org.openstreetmap.josm.data.Bounds; 007import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; 008import org.openstreetmap.josm.tools.JosmRuntimeException; 009 010/** 011 * Azimuthal Equidistant projection. 012 * <p> 013 * This implementation does not include the Guam or Micronesia variants. 014 * 015 * @author Gerald Evenden (original PROJ.4 implementation in C) 016 * @author Ben Caradoc-Davies (Transient Software Limited) 017 * @see <a href="https://pubs.er.usgs.gov/publication/pp1395"><em>Map Projections: A Working Manual</em>, Snyder (1987), pages 191-202</a> 018 * @see <a href="http://geotiff.maptools.org/proj_list/azimuthal_equidistant.html">PROJ.4 notes on parameters</a> 019 * @see <a href="https://github.com/OSGeo/proj.4/blob/master/src/PJ_aeqd.c">PROJ.4 implemention in C</a> 020 * @see <a href="https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection">Wikipedia</a> 021 * @see <a href="http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html">Wolfram Alpha</a> 022 * @since 13598 023 */ 024public class AzimuthalEquidistant extends AbstractProj { 025 026 /** 027 * Less strict tolerance. 028 */ 029 public static final double EPS10 = 1.e-10; 030 031 /** 032 * Stricter tolerance. 033 */ 034 public static final double TOL = 1.e-14; 035 036 /** 037 * Half of π. 038 */ 039 public static final double HALF_PI = Math.PI / 2; 040 041 /** 042 * The four possible modes or aspects of the projection. 043 */ 044 public enum Mode { 045 /** North pole */ 046 NORTH_POLAR, 047 /** South pole */ 048 SOUTH_POLAR, 049 /** Equator */ 050 EQUATORIAL, 051 /** Oblique */ 052 OBLIQUE; 053 } 054 055 /** 056 * Length of semi-major axis, in metres. This is named '<var>a</var>' or '<var>R</var>' 057 * (Radius in spherical cases) in Snyder. 058 */ 059 protected double semiMajor; 060 061 /** 062 * Length of semi-minor axis, in metres. This is named '<var>b</var>' in Snyder. 063 */ 064 protected double semiMinor; 065 066 /** 067 * Central longitude in <u>radians</u>. Default value is 0, the Greenwich meridian. 068 * This is called '<var>lambda0</var>' in Snyder. 069 */ 070 protected double centralMeridian; 071 072 /** 073 * Latitude of origin in <u>radians</u>. Default value is 0, the equator. 074 * This is called '<var>phi0</var>' in Snyder. 075 */ 076 protected double latitudeOfOrigin; 077 078 /** 079 * The mode or aspect of the projection. 080 */ 081 protected Mode mode; 082 083 /** 084 * Geodesic calculator used for this projection. Not used and set to null for polar projections. 085 */ 086 //protected Geodesic geodesic; // See https://josm.openstreetmap.de/ticket/16129#comment:21 087 088 /** 089 * The sine of the central latitude of the projection. 090 */ 091 protected double sinph0; 092 093 /** 094 * The cosine of the central latitude of the projection. 095 */ 096 protected double cosph0; 097 098 /** 099 * Meridian distance from the equator to the pole. Not used and set to NaN for non-polar projections. 100 */ 101 protected double mp; 102 103 @Override 104 public String getName() { 105 return tr("Azimuthal Equidistant"); 106 } 107 108 @Override 109 public String getProj4Id() { 110 return "aeqd"; 111 } 112 113 @Override 114 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 115 super.initialize(params); 116 if (params.lon0 == null) 117 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lon_0")); 118 if (params.lat0 == null) 119 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 120 centralMeridian = Math.toRadians(params.lon0); 121 latitudeOfOrigin = Math.toRadians(params.lat0); 122 semiMajor = params.ellps.a; 123 semiMinor = params.ellps.b; 124 if (Math.abs(latitudeOfOrigin - HALF_PI) < EPS10) { 125 mode = Mode.NORTH_POLAR; 126 mp = mlfn(HALF_PI, 1, 0); 127 sinph0 = 1; 128 cosph0 = 0; 129 } else if (Math.abs(latitudeOfOrigin + HALF_PI) < EPS10) { 130 mode = Mode.SOUTH_POLAR; 131 mp = mlfn(-HALF_PI, -1, 0); 132 sinph0 = -1; 133 cosph0 = 0; 134 } else if (Math.abs(latitudeOfOrigin) < EPS10) { 135 mode = Mode.EQUATORIAL; 136 mp = Double.NaN; 137 sinph0 = 0; 138 cosph0 = 1; 139 //geodesic = new Geodesic(semiMajor, (semiMajor - semiMinor) / semiMajor); 140 throw new ProjectionConfigurationException("Equatorial AzimuthalEquidistant not yet supported"); 141 } else { 142 mode = Mode.OBLIQUE; 143 mp = Double.NaN; 144 sinph0 = Math.sin(latitudeOfOrigin); 145 cosph0 = Math.cos(latitudeOfOrigin); 146 //geodesic = new Geodesic(semiMajor, (semiMajor - semiMinor) / semiMajor); 147 throw new ProjectionConfigurationException("Oblique AzimuthalEquidistant not yet supported"); 148 } 149 } 150 151 @Override 152 public Bounds getAlgorithmBounds() { 153 return new Bounds(-89, -180, 89, 180, false); 154 } 155 156 @Override 157 public double[] project(double latRad, double lonRad) { 158 return spherical ? projectSpherical(latRad, lonRad) : projectEllipsoidal(latRad, lonRad); 159 } 160 161 @Override 162 public double[] invproject(double east, double north) { 163 return spherical ? invprojectSpherical(east, north) : invprojectEllipsoidal(east, north); 164 } 165 166 double[] projectSpherical(double latRad, double lonRad) { 167 double x = 0; 168 double y = 0; 169 double sinphi = Math.sin(latRad); 170 double cosphi = Math.cos(latRad); 171 double coslam = Math.cos(lonRad); 172 switch (mode) { 173 case EQUATORIAL: 174 case OBLIQUE: 175 if (mode == Mode.EQUATORIAL) { 176 y = cosphi * coslam; 177 } else { // Oblique 178 y = sinph0 * sinphi + cosph0 * cosphi * coslam; 179 } 180 if (Math.abs(Math.abs(y) - 1) < TOL) { 181 if (y < 0) { 182 throw new JosmRuntimeException("TOLERANCE_ERROR"); 183 } else { 184 x = 0; 185 y = 0; 186 } 187 } else { 188 y = Math.acos(y); 189 y /= Math.sin(y); 190 x = y * cosphi * Math.sin(lonRad); 191 y *= (mode == Mode.EQUATORIAL) ? sinphi 192 : (cosph0 * sinphi - sinph0 * cosphi * coslam); 193 } 194 break; 195 case NORTH_POLAR: 196 latRad = -latRad; 197 coslam = -coslam; 198 // fall through 199 case SOUTH_POLAR: 200 if (Math.abs(latRad - HALF_PI) < EPS10) { 201 throw new JosmRuntimeException("TOLERANCE_ERROR"); 202 } 203 y = HALF_PI + latRad; 204 x = y * Math.sin(lonRad); 205 y *= coslam; 206 break; 207 } 208 return new double[] {x, y}; 209 } 210 211 double[] invprojectSpherical(double east, double north) { 212 double x = east; 213 double y = north; 214 double lambda = 0; 215 double phi = 0; 216 double c_rh = Math.hypot(x, y); 217 if (c_rh > Math.PI) { 218 if (c_rh - EPS10 > Math.PI) { 219 throw new JosmRuntimeException("TOLERANCE_ERROR"); 220 } 221 } else if (c_rh < EPS10) { 222 phi = latitudeOfOrigin; 223 lambda = 0.; 224 } else { 225 if (mode == Mode.OBLIQUE || mode == Mode.EQUATORIAL) { 226 double sinc = Math.sin(c_rh); 227 double cosc = Math.cos(c_rh); 228 if (mode == Mode.EQUATORIAL) { 229 phi = aasin(y * sinc / c_rh); 230 x *= sinc; 231 y = cosc * c_rh; 232 } else { // Oblique 233 phi = aasin(cosc * sinph0 + y * sinc * cosph0 / c_rh); 234 y = (cosc - sinph0 * Math.sin(phi)) * c_rh; 235 x *= sinc * cosph0; 236 } 237 lambda = (y == 0) ? 0 : Math.atan2(x, y); 238 } else if (mode == Mode.NORTH_POLAR) { 239 phi = HALF_PI - c_rh; 240 lambda = Math.atan2(x, -y); 241 } else { // South Polar 242 phi = c_rh - HALF_PI; 243 lambda = Math.atan2(x, y); 244 } 245 } 246 return new double[] {phi, lambda}; 247 } 248 249 double[] projectEllipsoidal(double latRad, double lonRad) { 250 double x = 0; 251 double y = 0; 252 double coslam = Math.cos(lonRad); 253 double cosphi = Math.cos(latRad); 254 double sinphi = Math.sin(latRad); 255 switch (mode) { 256 case NORTH_POLAR: 257 coslam = -coslam; 258 // fall through 259 case SOUTH_POLAR: 260 double rho = Math.abs(mp - mlfn(latRad, sinphi, cosphi)); 261 x = rho * Math.sin(lonRad); 262 y = rho * coslam; 263 break; 264 case EQUATORIAL: 265 case OBLIQUE: 266 if (Math.abs(lonRad) < EPS10 && Math.abs(latRad - latitudeOfOrigin) < EPS10) { 267 x = 0; 268 y = 0; 269 break; 270 } 271 /*GeodesicData g = geodesic.Inverse(Math.toDegrees(latitudeOfOrigin), 272 Math.toDegrees(centralMeridian), Math.toDegrees(latRad), 273 Math.toDegrees(lonRad + centralMeridian)); 274 double azi1 = Math.toRadians(g.azi1); 275 x = g.s12 * Math.sin(azi1) / semiMajor; 276 y = g.s12 * Math.cos(azi1) / semiMajor;*/ 277 break; 278 } 279 return new double[] {x, y}; 280 } 281 282 double[] invprojectEllipsoidal(double east, double north) { 283 double x = east; 284 double y = north; 285 double lambda = 0; 286 double phi = 0; 287 double c = Math.hypot(x, y); 288 if (c < EPS10) { 289 phi = latitudeOfOrigin; 290 lambda = 0; 291 } else { 292 if (mode == Mode.OBLIQUE || mode == Mode.EQUATORIAL) { 293 /*double x2 = x * semiMajor; 294 double y2 = y * semiMajor; 295 double azi1 = Math.atan2(x2, y2); 296 double s12 = Math.sqrt(x2 * x2 + y2 * y2); 297 GeodesicData g = geodesic.Direct(Math.toDegrees(latitudeOfOrigin), 298 Math.toDegrees(centralMeridian), Math.toDegrees(azi1), s12); 299 phi = Math.toRadians(g.lat2); 300 lambda = Math.toRadians(g.lon2);*/ 301 lambda -= centralMeridian; 302 } else { // Polar 303 phi = invMlfn((mode == Mode.NORTH_POLAR) ? (mp - c) : (mp + c)); 304 lambda = Math.atan2(x, (mode == Mode.NORTH_POLAR) ? -y : y); 305 } 306 } 307 return new double[] {phi, lambda}; 308 } 309}