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