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 * The polar case of the stereographic projection. 012 * <p> 013 * In the proj.4 library, the code "stere" covers several variants of the 014 * Stereographic projection, depending on the latitude of natural origin 015 * (parameter lat_0). 016 * <p> 017 * 018 * In this file, only the polar case is implemented. This corresponds to 019 * EPSG:9810 (Polar Stereographic Variant A) and EPSG:9829 (Polar Stereographic 020 * Variant B). 021 * <p> 022 * 023 * It is required, that the latitude of natural origin has the value +/-90 degrees. 024 * <p> 025 * 026 * This class has been derived from the implementation of the Geotools project; 027 * git 8cbf52d, org.geotools.referencing.operation.projection.PolarStereographic 028 * at the time of migration. 029 * <p> 030 * 031 * <b>References:</b> 032 * <ul> 033 * <li>John P. Snyder (Map Projections - A Working Manual,<br> 034 * U.S. Geological Survey Professional Paper 1395, 1987)</li> 035 * <li>"Coordinate Conversions and Transformations including Formulas",<br> 036 * EPSG Guidence Note Number 7, Version 19.</li> 037 * <li>Gerald Evenden. <A HREF="http://members.bellatlantic.net/~vze2hc4d/proj4/sterea.pdf"> 038 * "Supplementary PROJ.4 Notes - Oblique Stereographic Alternative"</A></li> 039 * <li>Krakiwsky, E.J., D.B. Thomson, and R.R. Steeves. 1977. A Manual 040 * For Geodetic Coordinate Transformations in the Maritimes. 041 * Geodesy and Geomatics Engineering, UNB. Technical Report No. 48.</li> 042 * <li>Thomson, D.B., M.P. Mepham and R.R. Steeves. 1977. 043 * The Stereographic Double Projection. 044 * Geodesy and Geomatics Engineereng, UNB. Technical Report No. 46.</li> 045 * </ul> 046 * 047 * @author André Gosselin 048 * @author Martin Desruisseaux (PMO, IRD) 049 * @author Rueben Schulz 050 * 051 * @see <A HREF="http://mathworld.wolfram.com/StereographicProjection.html">Stereographic projection on MathWorld</A> 052 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/polar_stereographic.html">Polar_Stereographic</A> 053 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/oblique_stereographic.html">Oblique_Stereographic</A> 054 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/stereographic.html">Stereographic</A> 055 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/random_issues.html#stereographic">Some Random Stereographic Issues</A> 056 * 057 * @see DoubleStereographic 058 * @since 9419 059 */ 060public class PolarStereographic extends AbstractProj { 061 /** 062 * Maximum number of iterations for iterative computations. 063 */ 064 private static final int MAXIMUM_ITERATIONS = 15; 065 066 /** 067 * Difference allowed in iterative computations. 068 */ 069 private static final double ITERATION_TOLERANCE = 1E-10; 070 071 /** 072 * Maximum difference allowed when comparing real numbers. 073 */ 074 private static final double EPSILON = 1E-8; 075 076 /** 077 * A constant used in the transformations. 078 */ 079 private double k0; 080 081 /** 082 * {@code true} if this projection is for the south pole, or {@code false} 083 * if it is for the north pole. 084 */ 085 boolean southPole; 086 087 @Override 088 public String getName() { 089 return tr("Polar Stereographic"); 090 } 091 092 @Override 093 public String getProj4Id() { 094 return "stere"; 095 } 096 097 @Override 098 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 099 super.initialize(params); 100 if (params.lat0 == null) 101 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 102 if (params.lat0 != 90.0 && params.lat0 != -90.0) 103 throw new ProjectionConfigurationException( 104 tr("Polar Stereographic: Parameter ''{0}'' must be 90 or -90.", "lat_0")); 105 // Latitude of true scale, in radians 106 double latitudeTrueScale; 107 if (params.lat_ts == null) { 108 latitudeTrueScale = (params.lat0 < 0) ? -Math.PI/2 : Math.PI/2; 109 } else { 110 latitudeTrueScale = Utils.toRadians(params.lat_ts); 111 } 112 southPole = latitudeTrueScale < 0; 113 114 // Computes coefficients 115 double latitudeTrueScaleAbs = Math.abs(latitudeTrueScale); 116 if (Math.abs(latitudeTrueScaleAbs - Math.PI/2) >= EPSILON) { 117 final double t = Math.sin(latitudeTrueScaleAbs); 118 k0 = msfn(t, Math.cos(latitudeTrueScaleAbs)) / 119 tsfn(latitudeTrueScaleAbs, t); // Derives from (21-32 and 21-33) 120 } else { 121 // True scale at pole (part of (21-33)) 122 k0 = 2.0 / Math.sqrt(Math.pow(1+e, 1+e)* 123 Math.pow(1-e, 1-e)); 124 } 125 } 126 127 @Override 128 public double[] project(double y, double x) { 129 final double sinlat = Math.sin(y); 130 final double coslon = Math.cos(x); 131 final double sinlon = Math.sin(x); 132 if (southPole) { 133 final double rho = k0 * tsfn(-y, -sinlat); 134 x = rho * sinlon; 135 y = rho * coslon; 136 } else { 137 final double rho = k0 * tsfn(y, sinlat); 138 x = rho * sinlon; 139 y = -rho * coslon; 140 } 141 return new double[] {x, y}; 142 } 143 144 @Override 145 public double[] invproject(double x, double y) { 146 final double rho = Math.hypot(x, y); 147 if (southPole) { 148 y = -y; 149 } 150 /* 151 * Compute latitude using iterative technique. 152 */ 153 final double t = rho/k0; 154 final double halfe = e/2.0; 155 double phi0 = 0; 156 for (int i = MAXIMUM_ITERATIONS;;) { 157 final double esinphi = e * Math.sin(phi0); 158 final double phi = (Math.PI/2) - 2.0*Math.atan(t*Math.pow((1-esinphi)/(1+esinphi), halfe)); 159 if (Math.abs(phi-phi0) < ITERATION_TOLERANCE) { 160 x = (Math.abs(rho) < EPSILON) ? 0.0 : Math.atan2(x, -y); 161 y = southPole ? -phi : phi; 162 break; 163 } 164 phi0 = phi; 165 if (--i < 0) { 166 throw new IllegalStateException("no convergence for x="+x+", y="+y); 167 } 168 } 169 return new double[] {y, x}; 170 } 171 172 @Override 173 public Bounds getAlgorithmBounds() { 174 final double cut = 60; 175 if (southPole) { 176 return new Bounds(-90, -180, cut, 180, false); 177 } else { 178 return new Bounds(-cut, -180, 90, 180, false); 179 } 180 } 181}