1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection;
|
---|
3 |
|
---|
4 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
5 |
|
---|
6 | import org.openstreetmap.josm.data.Bounds;
|
---|
7 | import org.openstreetmap.josm.data.coor.EastNorth;
|
---|
8 | import org.openstreetmap.josm.data.coor.LatLon;
|
---|
9 |
|
---|
10 | public 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\u00E9union
|
---|
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\u00E9union 1947");
|
---|
216 | }
|
---|
217 |
|
---|
218 | public double getDefaultZoomInPPD() {
|
---|
219 | // this will set the map scaler to about 1000 m
|
---|
220 | return 10.02;
|
---|
221 | }
|
---|
222 | }
|
---|