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