1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
3 |
|
---|
4 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
5 |
|
---|
6 | import org.openstreetmap.josm.data.Bounds;
|
---|
7 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
8 | import org.openstreetmap.josm.tools.Utils;
|
---|
9 |
|
---|
10 | /**
|
---|
11 | * Albers Equal Area Projection (EPSG code 9822). This is a conic projection with parallels being
|
---|
12 | * unequally spaced arcs of concentric circles, more closely spaced at north and south edges of the
|
---|
13 | * map. Merideans are equally spaced radii of the same circles and intersect parallels at right
|
---|
14 | * angles. As the name implies, this projection minimizes distortion in areas.
|
---|
15 | * <p>
|
---|
16 | * The {@code "standard_parallel_2"} parameter is optional and will be given the same value as
|
---|
17 | * {@code "standard_parallel_1"} if not set (creating a 1 standard parallel projection).
|
---|
18 | * <p>
|
---|
19 | * <b>NOTE:</b>
|
---|
20 | * formulae used below are from a port, to Java, of the {@code proj4}
|
---|
21 | * package of the USGS survey. USGS work is acknowledged here.
|
---|
22 | * <p>
|
---|
23 | * This class has been derived from the implementation of the Geotools project;
|
---|
24 | * git 8cbf52d, org.geotools.referencing.operation.projection.AlbersEqualArea
|
---|
25 | * at the time of migration.
|
---|
26 | * <p>
|
---|
27 | * <b>References:</b>
|
---|
28 | * <ul>
|
---|
29 | * <li> Proj-4.4.7 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
|
---|
30 | * Relevent files are: PJ_aea.c, pj_fwd.c and pj_inv.c </li>
|
---|
31 | * <li> John P. Snyder (Map Projections - A Working Manual,
|
---|
32 | * U.S. Geological Survey Professional Paper 1395, 1987)</li>
|
---|
33 | * <li> "Coordinate Conversions and Transformations including Formulas",
|
---|
34 | * EPSG Guidence Note Number 7, Version 19.</li>
|
---|
35 | * </ul>
|
---|
36 | *
|
---|
37 | * @author Gerald I. Evenden (for original code in Proj4)
|
---|
38 | * @author Rueben Schulz
|
---|
39 | *
|
---|
40 | * @see <A HREF="http://mathworld.wolfram.com/AlbersEqual-AreaConicProjection.html">Albers Equal-Area Conic Projection on MathWorld</A>
|
---|
41 | * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/albers_equal_area_conic.html">"Albers_Conic_Equal_Area" on RemoteSensing.org</A>
|
---|
42 | * @see <A HREF="http://srmwww.gov.bc.ca/gis/bceprojection.html">British Columbia Albers Standard Projection</A>
|
---|
43 | *
|
---|
44 | * @since 9419
|
---|
45 | */
|
---|
46 | public class AlbersEqualArea extends AbstractProj {
|
---|
47 |
|
---|
48 | /**
|
---|
49 | * Maximum number of iterations for iterative computations.
|
---|
50 | */
|
---|
51 | private static final int MAXIMUM_ITERATIONS = 15;
|
---|
52 |
|
---|
53 | /**
|
---|
54 | * Difference allowed in iterative computations.
|
---|
55 | */
|
---|
56 | private static final double ITERATION_TOLERANCE = 1E-10;
|
---|
57 |
|
---|
58 | /**
|
---|
59 | * Maximum difference allowed when comparing real numbers.
|
---|
60 | */
|
---|
61 | private static final double EPSILON = 1E-6;
|
---|
62 |
|
---|
63 | /**
|
---|
64 | * Constants used by the spherical and elliptical Albers projection.
|
---|
65 | */
|
---|
66 | private double n, c, rho0;
|
---|
67 |
|
---|
68 | /**
|
---|
69 | * An error condition indicating iteration will not converge for the
|
---|
70 | * inverse ellipse. See Snyder (14-20)
|
---|
71 | */
|
---|
72 | private double ec;
|
---|
73 |
|
---|
74 | @Override
|
---|
75 | public String getName() {
|
---|
76 | return tr("Albers Equal Area");
|
---|
77 | }
|
---|
78 |
|
---|
79 | @Override
|
---|
80 | public String getProj4Id() {
|
---|
81 | return "aea";
|
---|
82 | }
|
---|
83 |
|
---|
84 | @Override
|
---|
85 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
86 | super.initialize(params);
|
---|
87 | if (params.lat0 == null)
|
---|
88 | throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
|
---|
89 | if (params.lat1 == null)
|
---|
90 | throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_1"));
|
---|
91 |
|
---|
92 | double lat0 = Utils.toRadians(params.lat0);
|
---|
93 | // Standards parallels in radians.
|
---|
94 | double phi1 = Utils.toRadians(params.lat1);
|
---|
95 | double phi2 = params.lat2 == null ? phi1 : Utils.toRadians(params.lat2);
|
---|
96 |
|
---|
97 | // Compute Constants
|
---|
98 | if (Math.abs(phi1 + phi2) < EPSILON) {
|
---|
99 | 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 | this.n = n;
|
---|
119 | }
|
---|
120 |
|
---|
121 | @Override
|
---|
122 | public double[] project(double y, double x) {
|
---|
123 | x *= n;
|
---|
124 | double rho = c - n * qsfn(Math.sin(y));
|
---|
125 | if (rho < 0.0) {
|
---|
126 | if (rho > -EPSILON) {
|
---|
127 | rho = 0.0;
|
---|
128 | } else {
|
---|
129 | throw new AssertionError();
|
---|
130 | }
|
---|
131 | }
|
---|
132 | rho = Math.sqrt(rho) / n;
|
---|
133 | // CHECKSTYLE.OFF: SingleSpaceSeparator
|
---|
134 | y = rho0 - rho * Math.cos(x);
|
---|
135 | x = rho * Math.sin(x);
|
---|
136 | // CHECKSTYLE.ON: SingleSpaceSeparator
|
---|
137 | return new double[] {x, y};
|
---|
138 | }
|
---|
139 |
|
---|
140 | @Override
|
---|
141 | public double[] invproject(double x, double y) {
|
---|
142 | y = rho0 - y;
|
---|
143 | double rho = Math.hypot(x, y);
|
---|
144 | if (rho > EPSILON) {
|
---|
145 | if (n < 0.0) {
|
---|
146 | rho = -rho;
|
---|
147 | x = -x;
|
---|
148 | y = -y;
|
---|
149 | }
|
---|
150 | x = Math.atan2(x, y) / n;
|
---|
151 | y = rho * n;
|
---|
152 | y = (c - y*y) / n;
|
---|
153 | if (Math.abs(y) <= ec) {
|
---|
154 | y = phi1(y);
|
---|
155 | } else {
|
---|
156 | y = (y < 0.0) ? -Math.PI/2.0 : Math.PI/2.0;
|
---|
157 | }
|
---|
158 | } else {
|
---|
159 | x = 0.0;
|
---|
160 | y = n > 0.0 ? Math.PI/2.0 : -Math.PI/2.0;
|
---|
161 | }
|
---|
162 | return new double[] {y, x};
|
---|
163 | }
|
---|
164 |
|
---|
165 | /**
|
---|
166 | * Iteratively solves equation (3-16) from Snyder.
|
---|
167 | *
|
---|
168 | * @param qs arcsin(q/2), used in the first step of iteration
|
---|
169 | * @return the latitude
|
---|
170 | */
|
---|
171 | public double phi1(final double qs) {
|
---|
172 | final double toneEs = 1 - e2;
|
---|
173 | double phi = Math.asin(0.5 * qs);
|
---|
174 | if (e < EPSILON) {
|
---|
175 | return phi;
|
---|
176 | }
|
---|
177 | for (int i = 0; i < MAXIMUM_ITERATIONS; i++) {
|
---|
178 | final double sinpi = Math.sin(phi);
|
---|
179 | final double cospi = Math.cos(phi);
|
---|
180 | final double con = e * sinpi;
|
---|
181 | final double com = 1.0 - con*con;
|
---|
182 | final double dphi = 0.5 * com*com / cospi *
|
---|
183 | (qs/toneEs - sinpi / com + 0.5/e * Math.log((1. - con) / (1. + con)));
|
---|
184 | phi += dphi;
|
---|
185 | if (Math.abs(dphi) <= ITERATION_TOLERANCE) {
|
---|
186 | return phi;
|
---|
187 | }
|
---|
188 | }
|
---|
189 | throw new IllegalStateException("no convergence for qs="+qs);
|
---|
190 | }
|
---|
191 |
|
---|
192 | /**
|
---|
193 | * Calculates q, Snyder equation (3-12)
|
---|
194 | *
|
---|
195 | * @param sinphi sin of the latitude q is calculated for
|
---|
196 | * @return q from Snyder equation (3-12)
|
---|
197 | */
|
---|
198 | private double qsfn(final double sinphi) {
|
---|
199 | final double oneEs = 1 - e2;
|
---|
200 | if (e >= EPSILON) {
|
---|
201 | final double con = e * sinphi;
|
---|
202 | return oneEs * (sinphi / (1. - con*con) -
|
---|
203 | (0.5/e) * Math.log((1.-con) / (1.+con)));
|
---|
204 | } else {
|
---|
205 | return sinphi + sinphi;
|
---|
206 | }
|
---|
207 | }
|
---|
208 |
|
---|
209 | @Override
|
---|
210 | public Bounds getAlgorithmBounds() {
|
---|
211 | return new Bounds(-89, -180, 89, 180, false);
|
---|
212 | }
|
---|
213 | }
|
---|