1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
3 |
|
---|
4 | import static java.lang.Math.*;
|
---|
5 |
|
---|
6 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
7 |
|
---|
8 | import org.openstreetmap.josm.data.projection.Ellipsoid;
|
---|
9 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
10 |
|
---|
11 | /**
|
---|
12 | * Implementation of the Lambert Conformal Conic projection.
|
---|
13 | *
|
---|
14 | * @author Pieren
|
---|
15 | */
|
---|
16 | public class LambertConformalConic implements Proj {
|
---|
17 |
|
---|
18 | protected Ellipsoid ellps;
|
---|
19 | protected double e;
|
---|
20 |
|
---|
21 | public static abstract class Parameters {
|
---|
22 | public final double latitudeOrigin;
|
---|
23 | public Parameters(double latitudeOrigin) {
|
---|
24 | this.latitudeOrigin = latitudeOrigin;
|
---|
25 | }
|
---|
26 | };
|
---|
27 |
|
---|
28 | public static class Parameters1SP extends Parameters {
|
---|
29 | public Parameters1SP(double latitudeOrigin) {
|
---|
30 | super(latitudeOrigin);
|
---|
31 | }
|
---|
32 | }
|
---|
33 |
|
---|
34 | public static class Parameters2SP extends Parameters {
|
---|
35 | public final double standardParallel1;
|
---|
36 | public final double standardParallel2;
|
---|
37 | public Parameters2SP(double latitudeOrigin, double standardParallel1, double standardParallel2) {
|
---|
38 | super(latitudeOrigin);
|
---|
39 | this.standardParallel1 = standardParallel1;
|
---|
40 | this.standardParallel2 = standardParallel2;
|
---|
41 | }
|
---|
42 | }
|
---|
43 |
|
---|
44 | private Parameters params;
|
---|
45 |
|
---|
46 | /**
|
---|
47 | * projection exponent
|
---|
48 | */
|
---|
49 | protected double n;
|
---|
50 | /**
|
---|
51 | * projection factor
|
---|
52 | */
|
---|
53 | protected double F;
|
---|
54 | /**
|
---|
55 | * radius of the parallel of latitude of the false origin (2SP) or at
|
---|
56 | * natural origin (1SP)
|
---|
57 | */
|
---|
58 | protected double r0;
|
---|
59 |
|
---|
60 | /**
|
---|
61 | * precision in iterative schema
|
---|
62 | */
|
---|
63 | protected static final double epsilon = 1e-12;
|
---|
64 |
|
---|
65 | @Override
|
---|
66 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
67 | ellps = params.ellps;
|
---|
68 | e = ellps.e;
|
---|
69 | if (params.lcc_n != null && params.lcc_F != null && params.lcc_r0 != null) {
|
---|
70 | initializeDirect(params.lcc_n, params.lcc_F, params.lcc_r0);
|
---|
71 | } else {
|
---|
72 | if (params.lat_0 == null)
|
---|
73 | throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
|
---|
74 | if (params.lat_1 != null && params.lat_2 != null) {
|
---|
75 | initialize2SP(params.lat_0, params.lat_1, params.lat_2);
|
---|
76 | } else {
|
---|
77 | initialize1SP(params.lat_0);
|
---|
78 | }
|
---|
79 | }
|
---|
80 | }
|
---|
81 |
|
---|
82 | /**
|
---|
83 | * Initialize for LCC with 2 standard parallels.
|
---|
84 | *
|
---|
85 | * @param lat_0 latitude of false origin (in degrees)
|
---|
86 | * @param lat_1 latitude of first standard parallel (in degrees)
|
---|
87 | * @param lat_2 latitude of second standard parallel (in degrees)
|
---|
88 | */
|
---|
89 | private void initialize2SP(double lat_0, double lat_1, double lat_2) {
|
---|
90 | this.params = new Parameters2SP(lat_0, lat_1, lat_2);
|
---|
91 |
|
---|
92 | final double m1 = m(toRadians(lat_1));
|
---|
93 | final double m2 = m(toRadians(lat_2));
|
---|
94 |
|
---|
95 | final double t1 = t(toRadians(lat_1));
|
---|
96 | final double t2 = t(toRadians(lat_2));
|
---|
97 | final double tf = t(toRadians(lat_0));
|
---|
98 |
|
---|
99 | n = (log(m1) - log(m2)) / (log(t1) - log(t2));
|
---|
100 | F = m1 / (n * pow(t1, n));
|
---|
101 | r0 = F * pow(tf, n);
|
---|
102 | }
|
---|
103 |
|
---|
104 | /**
|
---|
105 | * Initialize for LCC with 1 standard parallel.
|
---|
106 | *
|
---|
107 | * @param lat_0 latitude of natural origin (in degrees)
|
---|
108 | */
|
---|
109 | private void initialize1SP(double lat_0) {
|
---|
110 | this.params = new Parameters1SP(lat_0);
|
---|
111 | final double lat_0_rad = toRadians(lat_0);
|
---|
112 |
|
---|
113 | final double m0 = m(lat_0_rad);
|
---|
114 | final double t0 = t(lat_0_rad);
|
---|
115 |
|
---|
116 | n = sin(lat_0_rad);
|
---|
117 | F = m0 / (n * pow(t0, n));
|
---|
118 | r0 = F * pow(t0, n);
|
---|
119 | }
|
---|
120 |
|
---|
121 | /**
|
---|
122 | * Initialize LCC by providing the projection parameters directly.
|
---|
123 | *
|
---|
124 | * @param n see field n
|
---|
125 | * @param F see field F
|
---|
126 | * @param r0 see field r0
|
---|
127 | */
|
---|
128 | private void initializeDirect(double n, double F, double r0) {
|
---|
129 | this.params = null;
|
---|
130 | this.n = n;
|
---|
131 | this.F = F;
|
---|
132 | this.r0 = r0;
|
---|
133 | }
|
---|
134 |
|
---|
135 | /**
|
---|
136 | * auxiliary function t
|
---|
137 | */
|
---|
138 | protected double t(double lat_rad) {
|
---|
139 | return tan(PI/4 - lat_rad / 2.0)
|
---|
140 | / pow(( (1.0 - e * sin(lat_rad)) / (1.0 + e * sin(lat_rad))) , e/2);
|
---|
141 | }
|
---|
142 |
|
---|
143 | /**
|
---|
144 | * auxiliary function m
|
---|
145 | */
|
---|
146 | protected double m(double lat_rad) {
|
---|
147 | return cos(lat_rad) / (sqrt(1 - e * e * pow(sin(lat_rad), 2)));
|
---|
148 | }
|
---|
149 |
|
---|
150 | @Override
|
---|
151 | public String getName() {
|
---|
152 | return tr("Lambert Conformal Conic");
|
---|
153 | }
|
---|
154 |
|
---|
155 | @Override
|
---|
156 | public String getProj4Id() {
|
---|
157 | return "lcc";
|
---|
158 | }
|
---|
159 |
|
---|
160 | @Override
|
---|
161 | public double[] project(double phi, double 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 | }
|
---|