001// License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.data.projection.proj; 003 004import static java.lang.Math.PI; 005import static java.lang.Math.abs; 006import static java.lang.Math.atan; 007import static java.lang.Math.cos; 008import static java.lang.Math.exp; 009import static java.lang.Math.log; 010import static java.lang.Math.pow; 011import static java.lang.Math.sin; 012import static java.lang.Math.sqrt; 013import static java.lang.Math.tan; 014import static org.openstreetmap.josm.tools.I18n.tr; 015 016import org.openstreetmap.josm.data.Bounds; 017import org.openstreetmap.josm.data.projection.CustomProjection.Param; 018import org.openstreetmap.josm.data.projection.Ellipsoid; 019import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; 020import org.openstreetmap.josm.tools.Utils; 021 022/** 023 * Implementation of the Lambert Conformal Conic projection. 024 * 025 * @author Pieren 026 */ 027public class LambertConformalConic extends AbstractProj { 028 029 /** ellipsoid */ 030 protected Ellipsoid ellps; 031 032 /** 033 * Base class of Lambert Conformal Conic parameters. 034 */ 035 public static class Parameters { 036 /** latitude of origin */ 037 public final double latitudeOrigin; 038 039 /** 040 * Constructs a new {@code Parameters}. 041 * @param latitudeOrigin latitude of origin 042 */ 043 protected Parameters(double latitudeOrigin) { 044 this.latitudeOrigin = latitudeOrigin; 045 } 046 } 047 048 /** 049 * Parameters with a single standard parallel. 050 */ 051 public static class Parameters1SP extends Parameters { 052 /** 053 * Constructs a new {@code Parameters1SP}. 054 * @param latitudeOrigin latitude of origin 055 */ 056 public Parameters1SP(double latitudeOrigin) { 057 super(latitudeOrigin); 058 } 059 } 060 061 /** 062 * Parameters with two standard parallels. 063 */ 064 public static class Parameters2SP extends Parameters { 065 /** first standard parallel */ 066 public final double standardParallel1; 067 /** second standard parallel */ 068 public final double standardParallel2; 069 070 /** 071 * Constructs a new {@code Parameters2SP}. 072 * @param latitudeOrigin latitude of origin 073 * @param standardParallel1 first standard parallel 074 * @param standardParallel2 second standard parallel 075 */ 076 public Parameters2SP(double latitudeOrigin, double standardParallel1, double standardParallel2) { 077 super(latitudeOrigin); 078 this.standardParallel1 = standardParallel1; 079 this.standardParallel2 = standardParallel2; 080 } 081 } 082 083 private Parameters params; 084 085 /** 086 * projection exponent 087 */ 088 protected double n; 089 /** 090 * projection factor 091 */ 092 protected double f; 093 /** 094 * radius of the parallel of latitude of the false origin (2SP) or at 095 * natural origin (1SP) 096 */ 097 protected double r0; 098 099 /** 100 * precision in iterative schema 101 */ 102 protected static final double epsilon = 1e-12; 103 104 @Override 105 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 106 super.initialize(params); 107 ellps = params.ellps; 108 if (params.lat0 == null) 109 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", Param.lat_0.key)); 110 if (params.lat1 != null && params.lat2 != null) { 111 initialize2SP(params.lat0, params.lat1, params.lat2); 112 } else { 113 initialize1SP(params.lat0); 114 } 115 } 116 117 /** 118 * Initialize for LCC with 2 standard parallels. 119 * 120 * @param lat0 latitude of false origin (in degrees) 121 * @param lat1 latitude of first standard parallel (in degrees) 122 * @param lat2 latitude of second standard parallel (in degrees) 123 */ 124 private void initialize2SP(double lat0, double lat1, double lat2) { 125 this.params = new Parameters2SP(lat0, lat1, lat2); 126 127 final double m1 = m(Utils.toRadians(lat1)); 128 final double m2 = m(Utils.toRadians(lat2)); 129 130 final double t1 = t(Utils.toRadians(lat1)); 131 final double t2 = t(Utils.toRadians(lat2)); 132 final double tf = t(Utils.toRadians(lat0)); 133 134 n = (log(m1) - log(m2)) / (log(t1) - log(t2)); 135 f = m1 / (n * pow(t1, n)); 136 r0 = f * pow(tf, n); 137 } 138 139 /** 140 * Initialize for LCC with 1 standard parallel. 141 * 142 * @param lat0 latitude of natural origin (in degrees) 143 */ 144 private void initialize1SP(double lat0) { 145 this.params = new Parameters1SP(lat0); 146 final double lat0rad = Utils.toRadians(lat0); 147 148 final double m0 = m(lat0rad); 149 final double t0 = t(lat0rad); 150 151 n = sin(lat0rad); 152 f = m0 / (n * pow(t0, n)); 153 r0 = f * pow(t0, n); 154 } 155 156 /** 157 * auxiliary function t 158 * @param latRad latitude in radians 159 * @return result 160 */ 161 protected double t(double latRad) { 162 return tan(PI/4 - latRad / 2.0) 163 / pow((1.0 - e * sin(latRad)) / (1.0 + e * sin(latRad)), e/2); 164 } 165 166 /** 167 * auxiliary function m 168 * @param latRad latitude in radians 169 * @return result 170 */ 171 protected double m(double latRad) { 172 return cos(latRad) / (sqrt(1 - e * e * pow(sin(latRad), 2))); 173 } 174 175 @Override 176 public String getName() { 177 return tr("Lambert Conformal Conic"); 178 } 179 180 @Override 181 public String getProj4Id() { 182 return "lcc"; 183 } 184 185 @Override 186 public double[] project(double phi, double lambda) { 187 double sinphi = sin(phi); 188 double l = (0.5*log((1+sinphi)/(1-sinphi))) - e/2*log((1+e*sinphi)/(1-e*sinphi)); 189 double r = f*exp(-n*l); 190 double gamma = n*lambda; 191 double x = r*sin(gamma); 192 double y = r0 - r*cos(gamma); 193 return new double[] {x, y}; 194 } 195 196 @Override 197 public double[] invproject(double east, double north) { 198 double r = sqrt(pow(east, 2) + pow(north-r0, 2)); 199 double gamma = atan(east / (r0-north)); 200 double lambda = gamma/n; 201 double latIso = (-1/n) * log(abs(r/f)); 202 double phi = ellps.latitude(latIso, e, epsilon); 203 return new double[] {phi, lambda}; 204 } 205 206 /** 207 * Returns projection parameters. 208 * @return projection parameters 209 */ 210 public final Parameters getParameters() { 211 return params; 212 } 213 214 @Override 215 public Bounds getAlgorithmBounds() { 216 double lat; 217 if (params instanceof Parameters2SP) { 218 Parameters2SP p2p = (Parameters2SP) params; 219 lat = (p2p.standardParallel1 + p2p.standardParallel2) / 2; 220 } else { 221 lat = params.latitudeOrigin; 222 } 223 double minlat = Math.max(lat - 60, -89); 224 double maxlat = Math.min(lat + 60, 89); 225 return new Bounds(minlat, -85, maxlat, 85, false); 226 } 227}