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; 008 009/** 010 * Lambert Azimuthal Equal Area (EPSG code 9820). 011 * <p> 012 * This class has been derived from the implementation of the Geotools project; 013 * git 8cbf52d, org.geotools.referencing.operation.projection.LambertAzimuthalEqualArea 014 * at the time of migration. 015 * <p> 016 * <b>References:</b> 017 * <ul> 018 * <li> A. Annoni, C. Luzet, E.Gubler and J. Ihde - Map Projections for Europe</li> 019 * <li> John P. Snyder (Map Projections - A Working Manual, 020 * U.S. Geological Survey Professional Paper 1395)</li> 021 * </ul> 022 * 023 * @author Gerald Evenden (for original code in Proj4) 024 * @author Beate Stollberg 025 * @author Martin Desruisseaux 026 * 027 * @see <A HREF="http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html">Lambert Azimuthal Equal-Area Projection</A> 028 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_azimuthal_equal_area.html">"Lambert_Azimuthal_Equal_Area"</A> 029 */ 030public class LambertAzimuthalEqualArea extends AbstractProj { 031 032 /** Maximum difference allowed when comparing real numbers. */ 033 private static final double EPSILON = 1E-7; 034 035 /** Epsilon for the comparison of small quantities. */ 036 private static final double FINE_EPSILON = 1E-10; 037 038 /** Epsilon for the comparison of latitudes. */ 039 private static final double EPSILON_LATITUDE = 1E-10; 040 041 /** Constants for authalic latitude. */ 042 private static final double P00 = 0.33333333333333333333; 043 private static final double P01 = 0.17222222222222222222; 044 private static final double P02 = 0.10257936507936507936; 045 private static final double P10 = 0.06388888888888888888; 046 private static final double P11 = 0.06640211640211640211; 047 private static final double P20 = 0.01641501294219154443; 048 049 /** The projection mode. */ 050 private enum Mode { OBLIQUE, EQUATORIAL, NORTH_POLE, SOUTH_POLE } 051 052 /** The projection mode for this particular instance. */ 053 private Mode mode; 054 055 /** Constant parameters. */ 056 private double sinb1, cosb1, xmf, ymf, mmf, qp, dd, rq; 057 058 /** Coefficients for authalic latitude. */ 059 private double aPA0, aPA1, aPA2; 060 061 private double latitudeOfOrigin; 062 063 @Override 064 public String getName() { 065 return tr("Lambert Azimuthal Equal Area"); 066 } 067 068 @Override 069 public String getProj4Id() { 070 return "laea"; 071 } 072 073 @Override 074 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 075 super.initialize(params); 076 077 if (params.lat0 == null) 078 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 079 080 latitudeOfOrigin = Math.toRadians(params.lat0); 081 /* 082 * Detects the mode (oblique, etc.). 083 */ 084 final double t = Math.abs(latitudeOfOrigin); 085 if (Math.abs(t - Math.PI/2) < EPSILON_LATITUDE) { 086 mode = latitudeOfOrigin < 0.0 ? Mode.SOUTH_POLE : Mode.NORTH_POLE; 087 } else if (Math.abs(t) < EPSILON_LATITUDE) { 088 mode = Mode.EQUATORIAL; 089 } else { 090 mode = Mode.OBLIQUE; 091 } 092 /* 093 * Computes the constants for authalic latitude. 094 */ 095 final double es2 = e2 * e2; 096 final double es3 = e2 * es2; 097 aPA0 = P02 * es3 + P01 * es2 + P00 * e2; 098 aPA1 = P11 * es3 + P10 * es2; 099 aPA2 = P20 * es3; 100 101 final double sinphi; 102 qp = qsfn(1); 103 rq = Math.sqrt(0.5 * qp); 104 mmf = 0.5 / (1 - e2); 105 sinphi = Math.sin(latitudeOfOrigin); 106 sinb1 = qsfn(sinphi) / qp; 107 cosb1 = Math.sqrt(1.0 - sinb1 * sinb1); 108 switch (mode) { 109 case NORTH_POLE: // Fall through 110 case SOUTH_POLE: { 111 dd = 1.0; 112 xmf = ymf = rq; 113 break; 114 } 115 case EQUATORIAL: { 116 dd = 1.0 / rq; 117 xmf = 1.0; 118 ymf = 0.5 * qp; 119 break; 120 } 121 case OBLIQUE: { 122 dd = Math.cos(latitudeOfOrigin) / 123 (Math.sqrt(1.0 - e2 * sinphi * sinphi) * rq * cosb1); 124 xmf = rq * dd; 125 ymf = rq / dd; 126 break; 127 } 128 default: { 129 throw new AssertionError(mode); 130 } 131 } 132 } 133 134 @Override 135 public double[] project(final double phi, final double lambda) { 136 final double coslam = Math.cos(lambda); 137 final double sinlam = Math.sin(lambda); 138 final double sinphi = Math.sin(phi); 139 double q = qsfn(sinphi); 140 final double sinb, cosb, b, c, x, y; 141 switch (mode) { 142 case OBLIQUE: { 143 sinb = q / qp; 144 cosb = Math.sqrt(1.0 - sinb * sinb); 145 c = 1.0 + sinb1 * sinb + cosb1 * cosb * coslam; 146 b = Math.sqrt(2.0 / c); 147 y = ymf * b * (cosb1 * sinb - sinb1 * cosb * coslam); 148 x = xmf * b * cosb * sinlam; 149 break; 150 } 151 case EQUATORIAL: { 152 sinb = q / qp; 153 cosb = Math.sqrt(1.0 - sinb * sinb); 154 c = 1.0 + cosb * coslam; 155 b = Math.sqrt(2.0 / c); 156 y = ymf * b * sinb; 157 x = xmf * b * cosb * sinlam; 158 break; 159 } 160 case NORTH_POLE: { 161 c = (Math.PI / 2) + phi; 162 q = qp - q; 163 if (q >= 0.0) { 164 b = Math.sqrt(q); 165 x = b * sinlam; 166 y = coslam * -b; 167 } else { 168 x = y = 0.; 169 } 170 break; 171 } 172 case SOUTH_POLE: { 173 c = phi - (Math.PI / 2); 174 q = qp + q; 175 if (q >= 0.0) { 176 b = Math.sqrt(q); 177 x = b * sinlam; 178 y = coslam * +b; 179 } else { 180 x = y = 0.; 181 } 182 break; 183 } 184 default: { 185 throw new AssertionError(mode); 186 } 187 } 188 if (Math.abs(c) < EPSILON_LATITUDE) { 189 return new double[] {0, 0}; // this is an error, we should handle it somehow 190 } 191 return new double[] {x, y}; 192 } 193 194 @Override 195 public double[] invproject(double x, double y) { 196 final double lambda, phi; 197 switch (mode) { 198 case EQUATORIAL: // Fall through 199 case OBLIQUE: { 200 x /= dd; 201 y *= dd; 202 final double rho = Math.hypot(x, y); 203 if (rho < FINE_EPSILON) { 204 lambda = 0.0; 205 phi = latitudeOfOrigin; 206 } else { 207 double sCe, cCe, ab; 208 sCe = 2.0 * Math.asin(0.5 * rho / rq); 209 cCe = Math.cos(sCe); 210 sCe = Math.sin(sCe); 211 x *= sCe; 212 if (mode == Mode.OBLIQUE) { 213 ab = cCe * sinb1 + y * sCe * cosb1 / rho; 214 y = rho * cosb1 * cCe - y * sinb1 * sCe; 215 } else { 216 ab = y * sCe / rho; 217 y = rho * cCe; 218 } 219 lambda = Math.atan2(x, y); 220 phi = authlat(Math.asin(ab)); 221 } 222 break; 223 } 224 case NORTH_POLE: { 225 y = -y; 226 // Fall through 227 } 228 case SOUTH_POLE: { 229 final double q = x*x + y*y; 230 if (q == 0) { 231 lambda = 0.; 232 phi = latitudeOfOrigin; 233 } else { 234 double ab = 1.0 - q / qp; 235 if (mode == Mode.SOUTH_POLE) { 236 ab = -ab; 237 } 238 lambda = Math.atan2(x, y); 239 phi = authlat(Math.asin(ab)); 240 } 241 break; 242 } 243 default: { 244 throw new AssertionError(mode); 245 } 246 } 247 return new double[] {phi, lambda}; 248 } 249 250 251 /** 252 * Calculates <var>q</var>, Snyder equation (3-12) 253 * 254 * @param sinphi sin of the latitude <var>q</var> is calculated for. 255 * @return <var>q</var> from Snyder equation (3-12). 256 */ 257 private double qsfn(final double sinphi) { 258 if (e >= EPSILON) { 259 final double con = e * sinphi; 260 return (1.0 - e2) * (sinphi / (1.0 - con*con) - 261 (0.5 / e) * Math.log((1.0 - con) / (1.0 + con))); 262 } else { 263 return sinphi + sinphi; 264 } 265 } 266 267 /** 268 * Determines latitude from authalic latitude. 269 * @param beta authalic latitude 270 * @return corresponding latitude 271 */ 272 private double authlat(final double beta) { 273 final double t = beta + beta; 274 return beta + aPA0 * Math.sin(t) + aPA1 * Math.sin(t+t) + aPA2 * Math.sin(t+t+t); 275 } 276 277 @Override 278 public Bounds getAlgorithmBounds() { 279 return new Bounds(-89, -174, 89, 174, false); 280 } 281}