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