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