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