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