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, 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 sinphi = Math.sin(latitudeOfOrigin); 105 sinb1 = qsfn(sinphi) / qp; 106 cosb1 = Math.sqrt(1.0 - sinb1 * sinb1); 107 switch (mode) { 108 case NORTH_POLE: // Fall through 109 case SOUTH_POLE: { 110 dd = 1.0; 111 xmf = ymf = rq; 112 break; 113 } 114 case EQUATORIAL: { 115 dd = 1.0 / rq; 116 xmf = 1.0; 117 ymf = 0.5 * qp; 118 break; 119 } 120 case OBLIQUE: { 121 dd = Math.cos(latitudeOfOrigin) / 122 (Math.sqrt(1.0 - e2 * sinphi * sinphi) * rq * cosb1); 123 xmf = rq * dd; 124 ymf = rq / dd; 125 break; 126 } 127 default: { 128 throw new AssertionError(mode); 129 } 130 } 131 } 132 133 @Override 134 public double[] project(final double phi, final double lambda) { 135 final double coslam = Math.cos(lambda); 136 final double sinlam = Math.sin(lambda); 137 final double sinphi = Math.sin(phi); 138 double q = qsfn(sinphi); 139 final double sinb, cosb, b, c, x, y; 140 switch (mode) { 141 case OBLIQUE: { 142 sinb = q / qp; 143 cosb = Math.sqrt(1.0 - sinb * sinb); 144 c = 1.0 + sinb1 * sinb + cosb1 * cosb * coslam; 145 b = Math.sqrt(2.0 / c); 146 y = ymf * b * (cosb1 * sinb - sinb1 * cosb * coslam); 147 x = xmf * b * cosb * sinlam; 148 break; 149 } 150 case EQUATORIAL: { 151 sinb = q / qp; 152 cosb = Math.sqrt(1.0 - sinb * sinb); 153 c = 1.0 + cosb * coslam; 154 b = Math.sqrt(2.0 / c); 155 y = ymf * b * sinb; 156 x = xmf * b * cosb * sinlam; 157 break; 158 } 159 case NORTH_POLE: { 160 c = (Math.PI / 2) + phi; 161 q = qp - q; 162 if (q >= 0.0) { 163 b = Math.sqrt(q); 164 x = b * sinlam; 165 y = coslam * -b; 166 } else { 167 x = y = 0.; 168 } 169 break; 170 } 171 case SOUTH_POLE: { 172 c = phi - (Math.PI / 2); 173 q = qp + q; 174 if (q >= 0.0) { 175 b = Math.sqrt(q); 176 x = b * sinlam; 177 y = coslam * +b; 178 } else { 179 x = y = 0.; 180 } 181 break; 182 } 183 default: { 184 throw new AssertionError(mode); 185 } 186 } 187 if (Math.abs(c) < EPSILON_LATITUDE) { 188 return new double[] {0, 0}; // this is an error, we should handle it somehow 189 } 190 return new double[] {x, y}; 191 } 192 193 @Override 194 public double[] invproject(double x, double y) { 195 final double lambda, phi; 196 switch (mode) { 197 case EQUATORIAL: // Fall through 198 case OBLIQUE: { 199 x /= dd; 200 y *= dd; 201 final double rho = Math.hypot(x, y); 202 if (rho < FINE_EPSILON) { 203 lambda = 0.0; 204 phi = latitudeOfOrigin; 205 } else { 206 double sCe, cCe, ab; 207 sCe = 2.0 * Math.asin(0.5 * rho / rq); 208 cCe = Math.cos(sCe); 209 sCe = Math.sin(sCe); 210 x *= sCe; 211 if (mode == Mode.OBLIQUE) { 212 ab = cCe * sinb1 + y * sCe * cosb1 / rho; 213 y = rho * cosb1 * cCe - y * sinb1 * sCe; 214 } else { 215 ab = y * sCe / rho; 216 y = rho * cCe; 217 } 218 lambda = Math.atan2(x, y); 219 phi = authlat(Math.asin(ab)); 220 } 221 break; 222 } 223 case NORTH_POLE: { 224 y = -y; 225 // Fall through 226 } 227 case SOUTH_POLE: { 228 final double q = x*x + y*y; 229 if (q == 0) { 230 lambda = 0.; 231 phi = latitudeOfOrigin; 232 } else { 233 double ab = 1.0 - q / qp; 234 if (mode == Mode.SOUTH_POLE) { 235 ab = -ab; 236 } 237 lambda = Math.atan2(x, y); 238 phi = authlat(Math.asin(ab)); 239 } 240 break; 241 } 242 default: { 243 throw new AssertionError(mode); 244 } 245 } 246 return new double[] {phi, lambda}; 247 } 248 249 250 /** 251 * Calculates <var>q</var>, Snyder equation (3-12) 252 * 253 * @param sinphi sin of the latitude <var>q</var> is calculated for. 254 * @return <var>q</var> from Snyder equation (3-12). 255 */ 256 private double qsfn(final double sinphi) { 257 if (e >= EPSILON) { 258 final double con = e * sinphi; 259 return (1.0 - e2) * (sinphi / (1.0 - con*con) - 260 (0.5 / e) * Math.log((1.0 - con) / (1.0 + con))); 261 } else { 262 return sinphi + sinphi; 263 } 264 } 265 266 /** 267 * Determines latitude from authalic latitude. 268 * @param beta authalic latitude 269 * @return corresponding latitude 270 */ 271 private double authlat(final double beta) { 272 final double t = beta + beta; 273 return beta + aPA0 * Math.sin(t) + aPA1 * Math.sin(t+t) + aPA2 * Math.sin(t+t+t); 274 } 275 276 @Override 277 public Bounds getAlgorithmBounds() { 278 return new Bounds(-89, -174, 89, 174, false); 279 } 280}