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