001// License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.data.projection.proj; 003 004import static java.lang.Math.abs; 005import static java.lang.Math.atan; 006import static java.lang.Math.cos; 007import static java.lang.Math.exp; 008import static java.lang.Math.log; 009import static java.lang.Math.pow; 010import static java.lang.Math.sin; 011import static java.lang.Math.sqrt; 012import static org.openstreetmap.josm.tools.I18n.tr; 013import static org.openstreetmap.josm.tools.Utils.toRadians; 014 015import org.openstreetmap.josm.data.Bounds; 016import org.openstreetmap.josm.data.projection.CustomProjection.Param; 017import org.openstreetmap.josm.data.projection.Ellipsoid; 018import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; 019 020/** 021 * Lambert Conical Conformal Projection. Areas and shapes are deformed as one moves away from standard parallels. 022 * The angles are true in a limited area. This projection is used for the charts of North America, France and Belgium. 023 * <p> 024 * This implementation provides transforms for two cases of the lambert conic conformal projection: 025 * </p> 026 * <ul> 027 * <li>{@code Lambert_Conformal_Conic_1SP} (EPSG code 9801)</li> 028 * <li>{@code Lambert_Conformal_Conic_2SP} (EPSG code 9802)</li> 029 * </ul> 030 * <p> 031 * For the 1SP case the latitude of origin is used as the standard parallel (SP). 032 * To use 1SP with a latitude of origin different from the SP, use the 2SP and set the SP1 to the single SP. 033 * The {@code "standard_parallel_2"} parameter is optional and will be given the same value 034 * as {@code "standard_parallel_1"} if not set (creating a 1 standard parallel projection). 035 * </p> 036 * <b>References:</b> 037 * <ul> 038 * <li>John P. Snyder (Map Projections - A Working Manual,<br>U.S. Geological Survey Professional Paper 1395, 1987)</li> 039 * <li>"Coordinate Conversions and Transformations including Formulas",<br>EPSG Guidence Note Number 7, Version 19.</li> 040 * </ul> 041 * 042 * @author Pieren 043 * @author André Gosselin 044 * @author Martin Desruisseaux (PMO, IRD) 045 * @author Rueben Schulz 046 * 047 * @see <A HREF="http://mathworld.wolfram.com/LambertConformalConicProjection.html">Lambert conformal conic projection on MathWorld</A> 048 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html">lambert_conic_conformal_1sp</A> 049 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html">lambert_conic_conformal_2sp</A> 050 * 051 * @since 13639 (align implementation with proj.4 / GeoTools) 052 * @since 4285 (reworked from Lambert / LambertCC9Zones) 053 * @since 2304 (initial implementation by Pieren) 054 */ 055public class LambertConformalConic extends AbstractProj { 056 057 /** ellipsoid */ 058 protected Ellipsoid ellps; 059 060 /** 061 * Base class of Lambert Conformal Conic parameters. 062 */ 063 public static class Parameters { 064 /** latitude of origin */ 065 public final double latitudeOrigin; 066 067 /** 068 * Constructs a new {@code Parameters}. 069 * @param latitudeOrigin latitude of origin 070 */ 071 protected Parameters(double latitudeOrigin) { 072 this.latitudeOrigin = latitudeOrigin; 073 } 074 } 075 076 /** 077 * Parameters with a single standard parallel. 078 */ 079 public static class Parameters1SP extends Parameters { 080 /** 081 * Constructs a new {@code Parameters1SP}. 082 * @param latitudeOrigin latitude of origin 083 */ 084 public Parameters1SP(double latitudeOrigin) { 085 super(latitudeOrigin); 086 } 087 } 088 089 /** 090 * Parameters with two standard parallels. 091 */ 092 public static class Parameters2SP extends Parameters { 093 /** first standard parallel */ 094 public final double standardParallel1; 095 /** second standard parallel */ 096 public final double standardParallel2; 097 098 /** 099 * Constructs a new {@code Parameters2SP}. 100 * @param latitudeOrigin latitude of origin 101 * @param standardParallel1 first standard parallel 102 * @param standardParallel2 second standard parallel 103 */ 104 public Parameters2SP(double latitudeOrigin, double standardParallel1, double standardParallel2) { 105 super(latitudeOrigin); 106 this.standardParallel1 = standardParallel1; 107 this.standardParallel2 = standardParallel2; 108 } 109 } 110 111 private Parameters params; 112 113 /** 114 * projection exponent 115 */ 116 protected double n; 117 118 /** 119 * projection factor 120 */ 121 protected double f; 122 123 /** 124 * radius of the parallel of latitude of the false origin (2SP) or at natural origin (1SP) 125 */ 126 protected double r0; 127 128 /** 129 * precision in iterative schema 130 */ 131 protected static final double epsilon = 1e-12; 132 133 @Override 134 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 135 super.initialize(params); 136 ellps = params.ellps; 137 if (params.lat0 == null) 138 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", Param.lat_0.key)); 139 if (params.lat1 != null && params.lat2 != null) { 140 this.params = new Parameters2SP(params.lat0, params.lat1, params.lat2); 141 initialize2SP(toRadians(params.lat0), toRadians(params.lat1), toRadians(params.lat2)); 142 } else { 143 this.params = new Parameters1SP(params.lat0); 144 initialize1SP(toRadians(params.lat0)); 145 } 146 } 147 148 /** 149 * Initialize for LCC with 2 standard parallels. 150 * 151 * @param lat0 latitude of false origin (in radians) 152 * @param lat1 latitude of first standard parallel (in radians) 153 * @param lat2 latitude of second standard parallel (in radians) 154 */ 155 private void initialize2SP(double lat0, double lat1, double lat2) { 156 157 final double cosphi1 = cos(lat1); 158 final double sinphi1 = sin(lat1); 159 160 final double cosphi2 = cos(lat2); 161 final double sinphi2 = sin(lat2); 162 163 final double m1 = msfn(sinphi1, cosphi1); 164 final double m2 = msfn(sinphi2, cosphi2); 165 166 final double t0 = tsfn(lat0, sin(lat0)); 167 final double t1 = tsfn(lat1, sinphi1); 168 final double t2 = tsfn(lat2, sinphi2); 169 170 n = log(m1/m2) / log(t1/t2); 171 f = m1 * pow(t1, -n) / n; 172 r0 = f * pow(t0, n); 173 } 174 175 /** 176 * Initialize for LCC with 1 standard parallel. 177 * 178 * @param lat0 latitude of natural origin (in radians) 179 */ 180 private void initialize1SP(double lat0) { 181 final double m0 = msfn(sin(lat0), cos(lat0)); 182 final double t0 = tsfn(lat0, sin(lat0)); 183 184 n = sin(lat0); 185 f = m0 * pow(t0, -n) / n; 186 r0 = f * pow(t0, n); 187 } 188 189 @Override 190 public String getName() { 191 return tr("Lambert Conformal Conic"); 192 } 193 194 @Override 195 public String getProj4Id() { 196 return "lcc"; 197 } 198 199 @Override 200 public double[] project(double phi, double lambda) { 201 double sinphi = sin(phi); 202 double l = (0.5*log((1+sinphi)/(1-sinphi))) - e/2*log((1+e*sinphi)/(1-e*sinphi)); 203 double r = f*exp(-n*l); 204 double gamma = n*lambda; 205 double x = r*sin(gamma); 206 double y = r0 - r*cos(gamma); 207 return new double[] {x, y}; 208 } 209 210 @Override 211 public double[] invproject(double east, double north) { 212 double r = sqrt(pow(east, 2) + pow(north-r0, 2)); 213 double gamma = atan(east / (r0-north)); 214 double lambda = gamma/n; 215 double latIso = (-1/n) * log(abs(r/f)); 216 double phi = ellps.latitude(latIso, e, epsilon); 217 return new double[] {phi, lambda}; 218 } 219 220 /** 221 * Returns projection parameters. 222 * @return projection parameters 223 */ 224 public final Parameters getParameters() { 225 return params; 226 } 227 228 @Override 229 public Bounds getAlgorithmBounds() { 230 double lat; 231 if (params instanceof Parameters2SP) { 232 Parameters2SP p2p = (Parameters2SP) params; 233 lat = (p2p.standardParallel1 + p2p.standardParallel2) / 2; 234 } else { 235 lat = params.latitudeOrigin; 236 } 237 double minlat = Math.max(lat - 60, -89); 238 double maxlat = Math.min(lat + 60, 89); 239 return new Bounds(minlat, -85, maxlat, 85, false); 240 } 241}