1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection;
|
---|
3 |
|
---|
4 | /**
|
---|
5 | * This class is not a direct implementation of Projection. It collects all methods used
|
---|
6 | * for the projection of the French departements in the Caribbean Sea UTM zone 20N
|
---|
7 | * but using each time different local geodesic settings for the 7 parameters transformation algorithm.
|
---|
8 | */
|
---|
9 | import org.openstreetmap.josm.data.coor.EastNorth;
|
---|
10 | import org.openstreetmap.josm.data.coor.LatLon;
|
---|
11 |
|
---|
12 | public class UTM_20N_France_DOM {
|
---|
13 |
|
---|
14 | /**
|
---|
15 | * false east in meters (constant)
|
---|
16 | */
|
---|
17 | private static final double Xs = 500000.0;
|
---|
18 | /**
|
---|
19 | * false north in meters (0 in northern hemisphere, 10000000 in southern
|
---|
20 | * hemisphere)
|
---|
21 | */
|
---|
22 | private static double Ys = 0;
|
---|
23 | /**
|
---|
24 | * origin meridian longitude
|
---|
25 | */
|
---|
26 | protected double lg0;
|
---|
27 | /**
|
---|
28 | * UTM zone (from 1 to 60)
|
---|
29 | */
|
---|
30 | int zone = 20;
|
---|
31 | /**
|
---|
32 | * whether north or south hemisphere
|
---|
33 | */
|
---|
34 | private boolean isNorth = true;
|
---|
35 | /**
|
---|
36 | * 7 parameters transformation
|
---|
37 | */
|
---|
38 | double tx = 0.0;
|
---|
39 | double ty = 0.0;
|
---|
40 | double tz = 0.0;
|
---|
41 | double rx = 0;
|
---|
42 | double ry = 0;
|
---|
43 | double rz = 0;
|
---|
44 | double scaleDiff = 0;
|
---|
45 | /**
|
---|
46 | * precision in iterative schema
|
---|
47 | */
|
---|
48 | public static final double epsilon = 1e-11;
|
---|
49 | public final static double DEG_TO_RAD = Math.PI / 180;
|
---|
50 | public final static double RAD_TO_DEG = 180 / Math.PI;
|
---|
51 |
|
---|
52 | public UTM_20N_France_DOM(double[] translation, double[] rotation, double scaleDiff) {
|
---|
53 | this.tx = translation[0];
|
---|
54 | this.ty = translation[1];
|
---|
55 | this.tz = translation[2];
|
---|
56 | this.rx = rotation[0]/206264.806247096355; // seconds to radian
|
---|
57 | this.ry = rotation[1]/206264.806247096355;
|
---|
58 | this.rz = rotation[2]/206264.806247096355;
|
---|
59 | this.scaleDiff = scaleDiff;
|
---|
60 | }
|
---|
61 | public EastNorth latlon2eastNorth(LatLon p) {
|
---|
62 | // translate ellipsoid GRS80 (WGS83) => reference ellipsoid geographic
|
---|
63 | LatLon geo = GRS802Hayford(p);
|
---|
64 |
|
---|
65 | // reference ellipsoid geographic => UTM projection
|
---|
66 | return MTProjection(geo, Ellipsoid.hayford.a, Ellipsoid.hayford.e);
|
---|
67 | }
|
---|
68 |
|
---|
69 | /**
|
---|
70 | * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to UTM
|
---|
71 | * geographic, (ellipsoid Hayford)
|
---|
72 | */
|
---|
73 | private LatLon GRS802Hayford(LatLon wgs) {
|
---|
74 | double lat = Math.toRadians(wgs.lat()); // degree to radian
|
---|
75 | double lon = Math.toRadians(wgs.lon());
|
---|
76 | // WGS84 geographic => WGS84 cartesian
|
---|
77 | double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
|
---|
78 | double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
|
---|
79 | double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
|
---|
80 | double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
|
---|
81 | // translation
|
---|
82 | double coord[] = invSevenParametersTransformation(X, Y, Z);
|
---|
83 | // UTM cartesian => UTM geographic
|
---|
84 | return Geographic(coord[0], coord[1], coord[2], Ellipsoid.hayford);
|
---|
85 | }
|
---|
86 |
|
---|
87 | /**
|
---|
88 | * initializes from cartesian coordinates
|
---|
89 | *
|
---|
90 | * @param X
|
---|
91 | * 1st coordinate in meters
|
---|
92 | * @param Y
|
---|
93 | * 2nd coordinate in meters
|
---|
94 | * @param Z
|
---|
95 | * 3rd coordinate in meters
|
---|
96 | * @param ell
|
---|
97 | * reference ellipsoid
|
---|
98 | */
|
---|
99 | private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
|
---|
100 | double norm = Math.sqrt(X * X + Y * Y);
|
---|
101 | double lg = 2.0 * Math.atan(Y / (X + norm));
|
---|
102 | double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
|
---|
103 | double delta = 1.0;
|
---|
104 | while (delta > epsilon) {
|
---|
105 | double s2 = Math.sin(lt);
|
---|
106 | s2 *= s2;
|
---|
107 | double l = Math.atan((Z / norm)
|
---|
108 | / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
|
---|
109 | delta = Math.abs(l - lt);
|
---|
110 | lt = l;
|
---|
111 | }
|
---|
112 | double s2 = Math.sin(lt);
|
---|
113 | s2 *= s2;
|
---|
114 | // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
|
---|
115 | return new LatLon(lt, lg);
|
---|
116 | }
|
---|
117 |
|
---|
118 | /**
|
---|
119 | * initalizes from geographic coordinates
|
---|
120 | *
|
---|
121 | * @param coord geographic coordinates triplet
|
---|
122 | * @param a reference ellipsoid long axis
|
---|
123 | * @param e reference ellipsoid excentricity
|
---|
124 | */
|
---|
125 | private EastNorth MTProjection(LatLon coord, double a, double e) {
|
---|
126 | double n = 0.9996 * a;
|
---|
127 | Ys = (coord.lat() >= 0.0) ? 0.0 : 10000000.0;
|
---|
128 | double r6d = Math.PI / 30.0;
|
---|
129 | //zone = (int) Math.floor((coord.lon() + Math.PI) / r6d) + 1;
|
---|
130 | lg0 = r6d * (zone - 0.5) - Math.PI;
|
---|
131 | double e2 = e * e;
|
---|
132 | double e4 = e2 * e2;
|
---|
133 | double e6 = e4 * e2;
|
---|
134 | double e8 = e4 * e4;
|
---|
135 | double C[] = {
|
---|
136 | 1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
|
---|
137 | e2/8.0 - e4/96.0 - 9.0*e6/1024.0 - 901.0*e8/184320.0,
|
---|
138 | 13.0*e4/768.0 + 17.0*e6/5120.0 - 311.0*e8/737280.0,
|
---|
139 | 61.0*e6/15360.0 + 899.0*e8/430080.0,
|
---|
140 | 49561.0*e8/41287680.0
|
---|
141 | };
|
---|
142 | double s = e * Math.sin(coord.lat());
|
---|
143 | double l = Math.log(Math.tan(Math.PI/4.0 + coord.lat()/2.0) *
|
---|
144 | Math.pow((1.0 - s) / (1.0 + s), e/2.0));
|
---|
145 | double phi = Math.asin(Math.sin(coord.lon() - lg0) /
|
---|
146 | ((Math.exp(l) + Math.exp(-l)) / 2.0));
|
---|
147 | double ls = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
|
---|
148 | double lambda = Math.atan(((Math.exp(l) - Math.exp(-l)) / 2.0) /
|
---|
149 | Math.cos(coord.lon() - lg0));
|
---|
150 |
|
---|
151 | double north = C[0] * lambda;
|
---|
152 | double east = C[0] * ls;
|
---|
153 | for(int k = 1; k < 5; k++) {
|
---|
154 | double r = 2.0 * k * lambda;
|
---|
155 | double m = 2.0 * k * ls;
|
---|
156 | double em = Math.exp(m);
|
---|
157 | double en = Math.exp(-m);
|
---|
158 | double sr = Math.sin(r)/2.0 * (em + en);
|
---|
159 | double sm = Math.cos(r)/2.0 * (em - en);
|
---|
160 | north += C[k] * sr;
|
---|
161 | east += C[k] * sm;
|
---|
162 | }
|
---|
163 | east *= n;
|
---|
164 | east += Xs;
|
---|
165 | north *= n;
|
---|
166 | north += Ys;
|
---|
167 | return new EastNorth(east, north);
|
---|
168 | }
|
---|
169 |
|
---|
170 | public LatLon eastNorth2latlon(EastNorth p) {
|
---|
171 | MTProjection(p.east(), p.north(), zone, isNorth);
|
---|
172 | LatLon geo = Geographic(p, Ellipsoid.hayford.a, Ellipsoid.hayford.e, 0.0 /* z */);
|
---|
173 |
|
---|
174 | // reference ellipsoid geographic => reference ellipsoid cartesian
|
---|
175 | //Hayford2GRS80(ellipsoid, geo);
|
---|
176 | double N = Ellipsoid.hayford.a / (Math.sqrt(1.0 - Ellipsoid.hayford.e2 * Math.sin(geo.lat()) * Math.sin(geo.lat())));
|
---|
177 | double X = (N /*+ h*/) * Math.cos(geo.lat()) * Math.cos(geo.lon());
|
---|
178 | double Y = (N /*+ h*/) * Math.cos(geo.lat()) * Math.sin(geo.lon());
|
---|
179 | double Z = (N * (1.0-Ellipsoid.hayford.e2) /*+ h*/) * Math.sin(geo.lat());
|
---|
180 |
|
---|
181 | // translation
|
---|
182 | double coord[] = sevenParametersTransformation(X, Y, Z);
|
---|
183 |
|
---|
184 | // WGS84 cartesian => WGS84 geographic
|
---|
185 | LatLon wgs = cart2LatLon(coord[0], coord[1], coord[2], Ellipsoid.GRS80);
|
---|
186 | return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
|
---|
187 | }
|
---|
188 |
|
---|
189 | /**
|
---|
190 | * initializes new projection coordinates (in north hemisphere)
|
---|
191 | *
|
---|
192 | * @param east east from origin in meters
|
---|
193 | * @param north north from origin in meters
|
---|
194 | * @param zone zone number (from 1 to 60)
|
---|
195 | * @param isNorth true in north hemisphere, false in south hemisphere
|
---|
196 | */
|
---|
197 | private void MTProjection(double east, double north, int zone, boolean isNorth) {
|
---|
198 | Ys = isNorth ? 0.0 : 10000000.0;
|
---|
199 | double r6d = Math.PI / 30.0;
|
---|
200 | lg0 = r6d * (zone - 0.5) - Math.PI;
|
---|
201 | }
|
---|
202 |
|
---|
203 | public double scaleFactor() {
|
---|
204 | return 1/Math.PI/2;
|
---|
205 | }
|
---|
206 |
|
---|
207 | /**
|
---|
208 | * initalizes from projected coordinates (Mercator transverse projection)
|
---|
209 | *
|
---|
210 | * @param coord projected coordinates pair
|
---|
211 | * @param e reference ellipsoid excentricity
|
---|
212 | * @param a reference ellipsoid long axis
|
---|
213 | * @param z altitude in meters
|
---|
214 | */
|
---|
215 | private LatLon Geographic(EastNorth coord, double a, double e, double z) {
|
---|
216 | double n = 0.9996 * a;
|
---|
217 | double e2 = e * e;
|
---|
218 | double e4 = e2 * e2;
|
---|
219 | double e6 = e4 * e2;
|
---|
220 | double e8 = e4 * e4;
|
---|
221 | double C[] = {
|
---|
222 | 1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
|
---|
223 | e2/8.0 + e4/48.0 + 7.0*e6/2048.0 + e8/61440.0,
|
---|
224 | e4/768.0 + 3.0*e6/1280.0 + 559.0*e8/368640.0,
|
---|
225 | 17.0*e6/30720.0 + 283.0*e8/430080.0,
|
---|
226 | 4397.0*e8/41287680.0
|
---|
227 | };
|
---|
228 | double l = (coord.north() - Ys) / (n * C[0]);
|
---|
229 | double ls = (coord.east() - Xs) / (n * C[0]);
|
---|
230 | double l0 = l;
|
---|
231 | double ls0 = ls;
|
---|
232 | for(int k = 1; k < 5; k++) {
|
---|
233 | double r = 2.0 * k * l0;
|
---|
234 | double m = 2.0 * k * ls0;
|
---|
235 | double em = Math.exp(m);
|
---|
236 | double en = Math.exp(-m);
|
---|
237 | double sr = Math.sin(r)/2.0 * (em + en);
|
---|
238 | double sm = Math.cos(r)/2.0 * (em - en);
|
---|
239 | l -= C[k] * sr;
|
---|
240 | ls -= C[k] * sm;
|
---|
241 | }
|
---|
242 | double lon = lg0 + Math.atan(((Math.exp(ls) - Math.exp(-ls)) / 2.0) /
|
---|
243 | Math.cos(l));
|
---|
244 | double phi = Math.asin(Math.sin(l) /
|
---|
245 | ((Math.exp(ls) + Math.exp(-ls)) / 2.0));
|
---|
246 | l = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
|
---|
247 | double lat = 2.0 * Math.atan(Math.exp(l)) - Math.PI / 2.0;
|
---|
248 | double lt0;
|
---|
249 | do {
|
---|
250 | lt0 = lat;
|
---|
251 | double s = e * Math.sin(lat);
|
---|
252 | lat = 2.0 * Math.atan(Math.pow((1.0 + s) / (1.0 - s), e/2.0) *
|
---|
253 | Math.exp(l)) - Math.PI / 2.0;
|
---|
254 | }
|
---|
255 | while(Math.abs(lat-lt0) >= epsilon);
|
---|
256 | //h = z;
|
---|
257 |
|
---|
258 | return new LatLon(lat, lon);
|
---|
259 | }
|
---|
260 |
|
---|
261 | /**
|
---|
262 | * initializes from cartesian coordinates
|
---|
263 | *
|
---|
264 | * @param X 1st coordinate in meters
|
---|
265 | * @param Y 2nd coordinate in meters
|
---|
266 | * @param Z 3rd coordinate in meters
|
---|
267 | * @param ell reference ellipsoid
|
---|
268 | */
|
---|
269 | private LatLon cart2LatLon(double X, double Y, double Z, Ellipsoid ell) {
|
---|
270 | double norm = Math.sqrt(X * X + Y * Y);
|
---|
271 | double lg = 2.0 * Math.atan(Y / (X + norm));
|
---|
272 | double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
|
---|
273 | double delta = 1.0;
|
---|
274 | while (delta > epsilon) {
|
---|
275 | double s2 = Math.sin(lt);
|
---|
276 | s2 *= s2;
|
---|
277 | double l = Math.atan((Z / norm)
|
---|
278 | / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
|
---|
279 | delta = Math.abs(l - lt);
|
---|
280 | lt = l;
|
---|
281 | }
|
---|
282 | double s2 = Math.sin(lt);
|
---|
283 | s2 *= s2;
|
---|
284 | // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
|
---|
285 | return new LatLon(lt, lg);
|
---|
286 | }
|
---|
287 |
|
---|
288 | /**
|
---|
289 | * 7 parameters transformation
|
---|
290 | * @param coord X, Y, Z in array
|
---|
291 | * @return transformed X, Y, Z in array
|
---|
292 | */
|
---|
293 | private double[] sevenParametersTransformation(double Xa, double Ya, double Za){
|
---|
294 | double Xb = tx + Xa*(1+scaleDiff) + Za*ry - Ya*rz;
|
---|
295 | double Yb = ty + Ya*(1+scaleDiff) + Xa*rz - Za*rx;
|
---|
296 | double Zb = tz + Za*(1+scaleDiff) + Ya*rx - Xa*ry;
|
---|
297 | return new double[]{Xb, Yb, Zb};
|
---|
298 | }
|
---|
299 |
|
---|
300 | /**
|
---|
301 | * 7 parameters inverse transformation
|
---|
302 | * @param coord X, Y, Z in array
|
---|
303 | * @return transformed X, Y, Z in array
|
---|
304 | */
|
---|
305 | private double[] invSevenParametersTransformation(double Xa, double Ya, double Za){
|
---|
306 | double Xb = (1-scaleDiff)*(-tx + Xa + ((-tz+Za)*(-ry) - (-ty+Ya)*(-rz)));
|
---|
307 | double Yb = (1-scaleDiff)*(-ty + Ya + ((-tx+Xa)*(-rz) - (-tz+Za)*(-rx)));
|
---|
308 | double Zb = (1-scaleDiff)*(-tz + Za + ((-ty+Ya)*(-rx) - (-tx+Xa)*(-ry)));
|
---|
309 | return new double[]{Xb, Yb, Zb};
|
---|
310 | }
|
---|
311 |
|
---|
312 | }
|
---|