Ticket #3214: GaussLaborde_Reunion.java

File GaussLaborde_Reunion.java, 8.8 KB (added by pieren, 17 years ago)
Line 
1// License: GPL. For details, see LICENSE file.
2package org.openstreetmap.josm.data.projection;
3
4import static org.openstreetmap.josm.tools.I18n.tr;
5
6import org.openstreetmap.josm.data.Bounds;
7import org.openstreetmap.josm.data.coor.EastNorth;
8import org.openstreetmap.josm.data.coor.LatLon;
9
10public class GaussLaborde_Reunion implements Projection {
11
12 private static final double lon0 = Math.toRadians(55.53333333333333333333);
13 private static final double lat0 = Math.toRadians(-21.11666666666666666666);
14 private static final double x0 = 160000.0;
15 private static final double y0 = 50000.0;
16 private static final double k0 = 1;
17
18 private static double sinLat0 = Math.sin(lat0);
19 private static double cosLat0 = Math.cos(lat0);
20 private static double n1 = Math.sqrt(1 + cosLat0*cosLat0*cosLat0*cosLat0*Ellipsoid.hayford.e2/(1-Ellipsoid.hayford.e2));
21 private static double phic = Math.asin(sinLat0/n1);
22 private static double c = Ellipsoid.hayford.latitudeIsometric(phic, 0) - n1*Ellipsoid.hayford.latitudeIsometric(lat0, Ellipsoid.hayford.e);
23 private static double n2 = k0*Ellipsoid.hayford.a*Math.sqrt(1-Ellipsoid.hayford.e2)/(1-Ellipsoid.hayford.e2*sinLat0*sinLat0);
24 private static double xs = x0;
25 private static double ys = y0 - n2*phic;
26 private static final double epsilon = 1e-11;
27 private static final double scaleDiff = -32.3241E-6;
28 private static final double Tx = 789.524;
29 private static final double Ty = -626.486;
30 private static final double Tz = -89.904;
31 private static final double rx = Math.toRadians(0.6006/3600);
32 private static final double ry = Math.toRadians(76.7946/3600);
33 private static final double rz = Math.toRadians(-10.5788/3600);
34 private static final double rx2 = rx*rx;
35 private static final double ry2 = ry*ry;
36 private static final double rz2 = rz*rz;
37
38 public LatLon eastNorth2latlon(EastNorth p) {
39 // plan Gauss-Laborde to geographic Piton-des-Neiges
40 LatLon geo = Geographic(p);
41
42 // geographic Piton-des-Neiges/Hayford to geographic WGS84/GRS80
43 LatLon wgs = PTN2GRS80(geo);
44 return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
45 }
46
47 /*
48 * Convert projected coordinates (Gauss-Laborde) to reference ellipsoid Hayforf geographic Piton-des-Neiges
49 * @return geographic coordinates in radian
50 */
51 private LatLon Geographic(EastNorth eastNorth) {
52 double dxn = (eastNorth.east()-xs)/n2;
53 double dyn = (eastNorth.north()-ys)/n2;
54 double lambda = Math.atan(sinh(dxn)/Math.cos(dyn));
55 double latIso = Ellipsoid.hayford.latitudeIsometric(Math.asin(Math.sin(dyn)/cosh(dxn)), 0);
56 double lon = lon0 + lambda/n1;
57 double lat = Ellipsoid.hayford.latitude((latIso-c)/n1, Ellipsoid.hayford.e, 1E-12);
58 return new LatLon(lat, lon);
59 }
60
61 /**
62 * Convert geographic Piton-des-Neiges ellipsoid Hayford to geographic WGS84/GRS80
63 * @param PTN in radian
64 * @return
65 */
66 private LatLon PTN2GRS80(LatLon PTN) {
67 double lat = PTN.lat();
68 double lon = PTN.lon();
69 double N = Ellipsoid.hayford.a / (Math.sqrt(1.0 - Ellipsoid.hayford.e2 * Math.sin(lat) * Math.sin(lat)));
70 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
71 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
72 double Z = (N * (1.0 - Ellipsoid.hayford.e2)/* + height */) * Math.sin(lat);
73
74 // translation
75 double coord[] = sevenParametersTransformation(X, Y, Z);
76
77 // WGS84 cartesian => WGS84 geographic
78 return cart2LatLon(coord[0], coord[1], coord[2], Ellipsoid.GRS80);
79 }
80
81 /**
82 * 7 parameters transformation
83 * @param coord X, Y, Z in array
84 * @return transformed X, Y, Z in array
85 */
86 private double[] sevenParametersTransformation(double Xa, double Ya, double Za){
87 double Xb = Tx + Xa*(1+scaleDiff) + Za*ry - Ya*rz;
88 double Yb = Ty + Ya*(1+scaleDiff) + Xa*rz - Za*rx;
89 double Zb = Tz + Za*(1+scaleDiff) + Ya*rx - Xa*ry;
90 return new double[]{Xb, Yb, Zb};
91 }
92
93 public EastNorth latlon2eastNorth(LatLon wgs) {
94 // translate ellipsoid GRS80 (WGS83) => reference ellipsoid geographic Réunion
95 LatLon geo = GRS802Hayford(wgs);
96
97 // reference ellipsoid geographic => GaussLaborde plan
98 return GaussLabordeProjection(geo);
99 }
100
101 private LatLon GRS802Hayford(LatLon wgs) {
102 double lat = Math.toRadians(wgs.lat()); // degree to radian
103 double lon = Math.toRadians(wgs.lon());
104 // WGS84 geographic => WGS84 cartesian
105 double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
106 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
107 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
108 double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
109 // translation
110 double coord[] = invSevenParametersTransformation(X, Y, Z);
111 // Gauss cartesian => Gauss geographic
112 return Geographic(coord[0], coord[1], coord[2], Ellipsoid.hayford);
113 }
114
115 /**
116 * 7 parameters inverse transformation
117 * @param coord X, Y, Z in array
118 * @return transformed X, Y, Z in array
119 */
120 private double[] invSevenParametersTransformation(double Vx, double Vy, double Vz){
121 Vx = Vx - Tx;
122 Vy = Vy - Ty;
123 Vz = Vz - Tz;
124 double e = 1 + scaleDiff;
125 double e2 = e*e;
126 double det = e*(e2 + rx2 + ry2 + rz2);
127 double Ux = ((e2 + rx2)*Vx + (e*rz + rx*ry)*Vy + (rx*rz - e*ry)*Vz)/det;
128 double Uy = ((-e*rz + rx*ry)*Vx + (e2 + ry2)*Vy + (e*rx + ry*rz)*Vz)/det;
129 double Uz = ((e*ry + rx*rz)*Vx + (-e*rx + ry*rz)*Vy + (e2 + rz2)*Vz)/det;
130 return new double[]{Ux, Uy, Uz};
131 }
132
133 private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
134 double norm = Math.sqrt(X * X + Y * Y);
135 double lg = 2.0 * Math.atan(Y / (X + norm));
136 double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
137 double delta = 1.0;
138 while (delta > epsilon) {
139 double s2 = Math.sin(lt);
140 s2 *= s2;
141 double l = Math.atan((Z / norm)
142 / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
143 delta = Math.abs(l - lt);
144 lt = l;
145 }
146 double s2 = Math.sin(lt);
147 s2 *= s2;
148 // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
149 return new LatLon(lt, lg);
150 }
151
152 private EastNorth GaussLabordeProjection(LatLon geo) {
153 double lambda = n1*(geo.lon()-lon0);
154 double latIso = c + n1*Ellipsoid.hayford.latitudeIsometric(geo.lat());
155 double x = xs + n2*Ellipsoid.hayford.latitudeIsometric(Math.asin(Math.sin(lambda)/((Math.exp(latIso) + Math.exp(-latIso))/2)),0);
156 double y = ys + n2*Math.atan(((Math.exp(latIso) - Math.exp(-latIso))/2)/(Math.cos(lambda)));
157 return new EastNorth(x, y);
158 }
159
160 /**
161 * initializes from cartesian coordinates
162 *
163 * @param X 1st coordinate in meters
164 * @param Y 2nd coordinate in meters
165 * @param Z 3rd coordinate in meters
166 * @param ell reference ellipsoid
167 */
168 private LatLon cart2LatLon(double X, double Y, double Z, Ellipsoid ell) {
169 double norm = Math.sqrt(X * X + Y * Y);
170 double lg = 2.0 * Math.atan(Y / (X + norm));
171 double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
172 double delta = 1.0;
173 while (delta > epsilon) {
174 double s2 = Math.sin(lt);
175 s2 *= s2;
176 double l = Math.atan((Z / norm)
177 / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
178 delta = Math.abs(l - lt);
179 lt = l;
180 }
181 double s2 = Math.sin(lt);
182 s2 *= s2;
183 // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
184 return new LatLon(lt, lg);
185 }
186
187 /*
188 * hyperbolic sine
189 */
190 public static final double sinh(double x) {
191 return (Math.exp(x) - Math.exp(-x))/2;
192 }
193 /*
194 * hyperbolic cosine
195 */
196 public static final double cosh(double x) {
197 return (Math.exp(x) + Math.exp(-x))/2;
198 }
199
200 public String getCacheDirectoryName() {
201 return this.toString();
202 }
203
204 public Bounds getWorldBoundsLatLon() {
205 return new Bounds(
206 new LatLon(-21.5, 55.14),
207 new LatLon(-20.76, 55.94));
208 }
209
210 public String toCode() {
211 return "EPSG::3727";
212 }
213
214 @Override public String toString() {
215 return tr("Gauss-Laborde Réunion 1947");
216 }
217
218}