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 @Override 074 public String getName() { 075 return tr("Albers Equal Area"); 076 } 077 078 @Override 079 public String getProj4Id() { 080 return "aea"; 081 } 082 083 @Override 084 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 085 super.initialize(params); 086 if (params.lat0 == null) 087 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 088 if (params.lat1 == null) 089 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_1")); 090 091 double lat0 = Math.toRadians(params.lat0); 092 // Standards parallels in radians. 093 double phi1 = Math.toRadians(params.lat1); 094 double phi2 = params.lat2 == null ? phi1 : Math.toRadians(params.lat2); 095 096 // Compute Constants 097 if (Math.abs(phi1 + phi2) < EPSILON) { 098 throw new ProjectionConfigurationException(tr("standard parallels are opposite")); 099 } 100 double sinphi = Math.sin(phi1); 101 double cosphi = Math.cos(phi1); 102 double n = sinphi; 103 boolean secant = Math.abs(phi1 - phi2) >= EPSILON; 104 double m1 = msfn(sinphi, cosphi); 105 double q1 = qsfn(sinphi); 106 if (secant) { // secant cone 107 sinphi = Math.sin(phi2); 108 cosphi = Math.cos(phi2); 109 double m2 = msfn(sinphi, cosphi); 110 double q2 = qsfn(sinphi); 111 n = (m1 * m1 - m2 * m2) / (q2 - q1); 112 } 113 c = m1 * m1 + n * q1; 114 rho0 = Math.sqrt(c - n * qsfn(Math.sin(lat0))) /n; 115 ec = 1.0 - .5 * (1.0-e2) * 116 Math.log((1.0 - e) / (1.0 + e)) / e; 117 this.n = n; 118 } 119 120 @Override 121 public double[] project(double y, double x) { 122 x *= n; 123 double rho = c - n * qsfn(Math.sin(y)); 124 if (rho < 0.0) { 125 if (rho > -EPSILON) { 126 rho = 0.0; 127 } else { 128 throw new AssertionError(); 129 } 130 } 131 rho = Math.sqrt(rho) / n; 132 y = rho0 - rho * Math.cos(x); 133 x = rho * Math.sin(x); 134 return new double[] {x, y}; 135 } 136 137 @Override 138 public double[] invproject(double x, double y) { 139 y = rho0 - y; 140 double rho = Math.hypot(x, y); 141 if (rho > EPSILON) { 142 if (n < 0.0) { 143 rho = -rho; 144 x = -x; 145 y = -y; 146 } 147 x = Math.atan2(x, y) / n; 148 y = rho * n; 149 y = (c - y*y) / n; 150 if (Math.abs(y) <= ec) { 151 y = phi1(y); 152 } else { 153 y = (y < 0.0) ? -Math.PI/2.0 : Math.PI/2.0; 154 } 155 } else { 156 x = 0.0; 157 y = n > 0.0 ? Math.PI/2.0 : -Math.PI/2.0; 158 } 159 return new double[] {y, x}; 160 } 161 162 /** 163 * Iteratively solves equation (3-16) from Snyder. 164 * 165 * @param qs arcsin(q/2), used in the first step of iteration 166 * @return the latitude 167 */ 168 public double phi1(final double qs) { 169 final double toneEs = 1 - e2; 170 double phi = Math.asin(0.5 * qs); 171 if (e < EPSILON) { 172 return phi; 173 } 174 for (int i = 0; i < MAXIMUM_ITERATIONS; i++) { 175 final double sinpi = Math.sin(phi); 176 final double cospi = Math.cos(phi); 177 final double con = e * sinpi; 178 final double com = 1.0 - con*con; 179 final double dphi = 0.5 * com*com / cospi * 180 (qs/toneEs - sinpi / com + 0.5/e * Math.log((1. - con) / (1. + con))); 181 phi += dphi; 182 if (Math.abs(dphi) <= ITERATION_TOLERANCE) { 183 return phi; 184 } 185 } 186 throw new IllegalStateException("no convergence for qs="+qs); 187 } 188 189 /** 190 * Calculates q, Snyder equation (3-12) 191 * 192 * @param sinphi sin of the latitude q is calculated for 193 * @return q from Snyder equation (3-12) 194 */ 195 private double qsfn(final double sinphi) { 196 final double oneEs = 1 - e2; 197 if (e >= EPSILON) { 198 final double con = e * sinphi; 199 return oneEs * (sinphi / (1. - con*con) - 200 (0.5/e) * Math.log((1.-con) / (1.+con))); 201 } else { 202 return sinphi + sinphi; 203 } 204 } 205 206 @Override 207 public Bounds getAlgorithmBounds() { 208 return new Bounds(-89, -180, 89, 180, false); 209 } 210}