Ticket #3214: UTM_20N_France_DOM.java

File UTM_20N_France_DOM.java, 11.7 KB (added by pieren, 17 years ago)
Line 
1// License: GPL. For details, see LICENSE file.
2package 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 */
9import org.openstreetmap.josm.data.coor.EastNorth;
10import org.openstreetmap.josm.data.coor.LatLon;
11
12public 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}