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 * Albers Equal Area Projection (EPSG code 9822). This is a conic projection with parallels being 011 * unequally spaced arcs of concentric circles, more closely spaced at north and south edges of the 012 * map. Merideans are equally spaced radii of the same circles and intersect parallels at right 013 * angles. As the name implies, this projection minimizes distortion in areas. 014 * <p> 015 * The {@code "standard_parallel_2"} parameter is optional and will be given the same value as 016 * {@code "standard_parallel_1"} if not set (creating a 1 standard parallel projection). 017 * <p> 018 * <b>NOTE:</b> 019 * formulae used below are from a port, to Java, of the {@code proj4} 020 * package of the USGS survey. USGS work is acknowledged here. 021 * <p> 022 * This class has been derived from the implementation of the Geotools project; 023 * git 8cbf52d, org.geotools.referencing.operation.projection.AlbersEqualArea 024 * at the time of migration. 025 * <p> 026 * <b>References:</b> 027 * <ul> 028 * <li> Proj-4.4.7 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br> 029 * Relevent files are: PJ_aea.c, pj_fwd.c and pj_inv.c </li> 030 * <li> John P. Snyder (Map Projections - A Working Manual, 031 * U.S. Geological Survey Professional Paper 1395, 1987)</li> 032 * <li> "Coordinate Conversions and Transformations including Formulas", 033 * EPSG Guidence Note Number 7, Version 19.</li> 034 * </ul> 035 * 036 * @author Gerald I. Evenden (for original code in Proj4) 037 * @author Rueben Schulz 038 * 039 * @see <A HREF="http://mathworld.wolfram.com/AlbersEqual-AreaConicProjection.html">Albers Equal-Area Conic Projection on MathWorld</A> 040 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/albers_equal_area_conic.html">"Albers_Conic_Equal_Area" on RemoteSensing.org</A> 041 * @see <A HREF="http://srmwww.gov.bc.ca/gis/bceprojection.html">British Columbia Albers Standard Projection</A> 042 * 043 * @since 9419 044 */ 045public class AlbersEqualArea extends AbstractProj { 046 047 /** 048 * Maximum number of iterations for iterative computations. 049 */ 050 private static final int MAXIMUM_ITERATIONS = 15; 051 052 /** 053 * Difference allowed in iterative computations. 054 */ 055 private static final double ITERATION_TOLERANCE = 1E-10; 056 057 /** 058 * Maximum difference allowed when comparing real numbers. 059 */ 060 private static final double EPSILON = 1E-6; 061 062 /** 063 * Constants used by the spherical and elliptical Albers projection. 064 */ 065 private double n, c, rho0; 066 067 /** 068 * An error condition indicating iteration will not converge for the 069 * inverse ellipse. See Snyder (14-20) 070 */ 071 private double ec; 072 073 /** 074 * Standards parallel 1 in radians. 075 */ 076 private double phi1; 077 078 /** 079 * Standards parallel 2 in radians. 080 */ 081 private double phi2; 082 083 @Override 084 public String getName() { 085 return tr("Albers Equal Area"); 086 } 087 088 @Override 089 public String getProj4Id() { 090 return "aea"; 091 } 092 093 @Override 094 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 095 super.initialize(params); 096 if (params.lat0 == null) 097 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 098 if (params.lat1 == null) 099 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_1")); 100 101 double lat0 = Math.toRadians(params.lat0); 102 phi1 = Math.toRadians(params.lat1); 103 phi2 = params.lat2 == null ? phi1 : Math.toRadians(params.lat2); 104 105 // Compute Constants 106 if (Math.abs(phi1 + phi2) < EPSILON) { 107 throw new ProjectionConfigurationException(tr("standard parallels are opposite")); 108 } 109 double sinphi = Math.sin(phi1); 110 double cosphi = Math.cos(phi1); 111 double n = sinphi; 112 boolean secant = Math.abs(phi1 - phi2) >= EPSILON; 113 double m1 = msfn(sinphi, cosphi); 114 double q1 = qsfn(sinphi); 115 if (secant) { // secant cone 116 sinphi = Math.sin(phi2); 117 cosphi = Math.cos(phi2); 118 double m2 = msfn(sinphi, cosphi); 119 double q2 = qsfn(sinphi); 120 n = (m1 * m1 - m2 * m2) / (q2 - q1); 121 } 122 c = m1 * m1 + n * q1; 123 rho0 = Math.sqrt(c - n * qsfn(Math.sin(lat0))) /n; 124 ec = 1.0 - .5 * (1.0-e2) * 125 Math.log((1.0 - e) / (1.0 + e)) / e; 126 this.n = n; 127 } 128 129 @Override 130 public double[] project(double y, double x) { 131 x *= n; 132 double rho = c - n * qsfn(Math.sin(y)); 133 if (rho < 0.0) { 134 if (rho > -EPSILON) { 135 rho = 0.0; 136 } else { 137 throw new RuntimeException(); 138 } 139 } 140 rho = Math.sqrt(rho) / n; 141 y = rho0 - rho * Math.cos(x); 142 x = rho * Math.sin(x); 143 return new double[] {x, y}; 144 } 145 146 @Override 147 public double[] invproject(double x, double y) { 148 y = rho0 - y; 149 double rho = Math.hypot(x, y); 150 if (rho > EPSILON) { 151 if (n < 0.0) { 152 rho = -rho; 153 x = -x; 154 y = -y; 155 } 156 x = Math.atan2(x, y) / n; 157 y = rho * n; 158 y = (c - y*y) / n; 159 if (Math.abs(y) <= ec) { 160 y = phi1(y); 161 } else { 162 y = (y < 0.0) ? -Math.PI/2.0 : Math.PI/2.0; 163 } 164 } else { 165 x = 0.0; 166 y = n > 0.0 ? Math.PI/2.0 : -Math.PI/2.0; 167 } 168 return new double[] {y, x}; 169 } 170 171 /** 172 * Iteratively solves equation (3-16) from Snyder. 173 * 174 * @param qs arcsin(q/2), used in the first step of iteration 175 * @return the latitude 176 */ 177 public double phi1(final double qs) { 178 final double tone_es = 1 - e2; 179 double phi = Math.asin(0.5 * qs); 180 if (e < EPSILON) { 181 return phi; 182 } 183 for (int i = 0; i < MAXIMUM_ITERATIONS; i++) { 184 final double sinpi = Math.sin(phi); 185 final double cospi = Math.cos(phi); 186 final double con = e * sinpi; 187 final double com = 1.0 - con*con; 188 final double dphi = 0.5 * com*com / cospi * 189 (qs/tone_es - sinpi / com + 0.5/e * Math.log((1. - con) / (1. + con))); 190 phi += dphi; 191 if (Math.abs(dphi) <= ITERATION_TOLERANCE) { 192 return phi; 193 } 194 } 195 throw new RuntimeException("no convergence for q="+qs); 196 } 197 198 /** 199 * Calculates q, Snyder equation (3-12) 200 * 201 * @param sinphi sin of the latitude q is calculated for 202 * @return q from Snyder equation (3-12) 203 */ 204 private double qsfn(final double sinphi) { 205 final double one_es = 1 - e2; 206 if (e >= EPSILON) { 207 final double con = e * sinphi; 208 return one_es * (sinphi / (1. - con*con) - 209 (0.5/e) * Math.log((1.-con) / (1.+con))); 210 } else { 211 return sinphi + sinphi; 212 } 213 } 214 215 @Override 216 public Bounds getAlgorithmBounds() { 217 return new Bounds(-89, -180, 89, 180, false); 218 } 219}