1 package org.opentrafficsim.core.gis;
2
3 import java.awt.geom.Point2D;
4 import java.awt.geom.Rectangle2D;
5 import java.util.Locale;
6
7 import org.opentrafficsim.base.logger.Cat;
8 import org.opentrafficsim.core.logger.SimLogger;
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 public final class TransformWGS84DutchRDNew
25 {
26
27
28 private static final double WGS84_WEST_LIMIT = 3.2;
29
30
31 private static final double WGS84_EAST_LIMIT = 7.3;
32
33
34 private static final double WGS84_SOUTH_LIMIT = 50.6;
35
36
37 private static final double WGS84_NORTH_LIMIT = 53.7;
38
39
40 private static final double RD_MINIMUM_X = 11000;
41
42
43 private static final double RD_MAXIMUM_X = 280000;
44
45
46 private static final double RD_MINIMUM_Y = 300000;
47
48
49 private static final double RD_MAXIMUM_Y = 630000;
50
51
52 private TransformWGS84DutchRDNew()
53 {
54
55 }
56
57
58
59
60
61
62
63 private static Point2D.Double ellipswgs842rd(final double wgs84East, final double wgs84North)
64 {
65 if (wgs84North > WGS84_NORTH_LIMIT || wgs84North < WGS84_SOUTH_LIMIT || wgs84East < WGS84_WEST_LIMIT
66 || (wgs84East > WGS84_EAST_LIMIT))
67 {
68 throw new IllegalArgumentException("ellipswgs842rd input out of range (" + wgs84East + ", " + wgs84North + ")");
69 }
70
71 @SuppressWarnings("checkstyle:nowhitespacebefore")
72
73 final double[][] r =
74 {
75 { 155000.00, 190094.945, -0.008, -32.391, 0.0 , },
76 { -0.705, -11832.228, 0.0 , 0.608, 0.0 , },
77 { 0.0 , -114.221, 0.0 , 0.148, 0.0 , },
78 { 0.0 , -2.340, 0.0 , 0.0 , 0.0 , },
79 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , }
80 };
81 @SuppressWarnings("checkstyle:nowhitespacebefore")
82 final double[][] s =
83 {
84 { 463000.00 , 0.433, 3638.893, 0.0 , 0.092, },
85 { 309056.544, -0.032, -157.984, 0.0 , -0.054, },
86 { 73.077, 0.0 , -6.439, 0.0 , 0.0 , },
87 { 59.788, 0.0 , 0.0 , 0.0 , 0.0 , },
88 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , }
89 };
90
91 double resultX = 0;
92 double resultY = 0;
93 double powNorth = 1;
94 double dNorth = 0.36 * (wgs84North - 52.15517440);
95 double dEast = 0.36 * (wgs84East - 5.38720621);
96
97 for (int p = 0; p < 5; p++)
98 {
99 double powEast = 1;
100 for (int q = 0; q < 5; q++)
101 {
102 resultX += r[p][q] * powEast * powNorth;
103 resultY += s[p][q] * powEast * powNorth;
104 powEast *= dEast;
105 }
106 powNorth *= dNorth;
107 }
108 return new Point2D.Double(resultX, resultY);
109 }
110
111
112
113
114
115
116
117
118 private static Point2D rd2ellipsWGS84(final double rdX, final double rdY) throws IllegalArgumentException
119 {
120 if (rdX < RD_MINIMUM_X || rdX > RD_MAXIMUM_X || rdY < RD_MINIMUM_Y || rdY > RD_MAXIMUM_Y)
121 {
122 throw new IllegalArgumentException(
123 "Location (" + rdX + "," + rdY + ") is not within the range " + "of the Dutch RD system");
124 }
125 final double dX = (rdX - 155000) / 100000;
126 final double dY = (rdY - 463000) / 100000;
127
128
129 @SuppressWarnings("checkstyle:nowhitespacebefore")
130 final double[][] k =
131 {
132 { 3600 * 52.15517440, 3235.65389, -0.24750, -0.06550, 0.0 , },
133 { -0.00738 , -0.00012, 0.0 , 0.0 , 0.0 , },
134 { -32.58297 , -0.84978, -0.01709, -0.00039, 0.0 , },
135 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , },
136 { 0.00530 , 0.00033, 0.0 , 0.0 , 0.0 , },
137 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , }
138 };
139 @SuppressWarnings("checkstyle:nowhitespacebefore")
140 final double[][] l = {
141 { 3600 * 5.38720621, 0.01199, 0.00022, 0.0 , 0.0 , },
142 { 5260.52916 , 105.94684, 2.45656, 0.05594, 0.00128, },
143 { -0.00022 , 0.0 , 0.0 , 0.0 , 0.0 , },
144 { -0.81885 , -0.05607, -0.00256, 0.0 , 0.0 , },
145 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , },
146 { 0.00026 , 0.0 , 0.0 , 0.0 , 0.0 , }
147 };
148
149 double resultNorth = 0;
150 double resultEast = 0;
151 double powX = 1;
152 for (int p = 0; p < 6; p++)
153 {
154 double powY = 1;
155 for (int q = 0; q < 5; q++)
156 {
157 resultNorth += k[p][q] * powX * powY / 3600;
158 resultEast += l[p][q] * powX * powY / 3600;
159 powY *= dY;
160 }
161 powX *= dX;
162 }
163 return new Point2D.Double(resultEast, resultNorth);
164 }
165
166
167
168
169
170
171
172 public static Point2D toWGS84(final Point2D local) throws IllegalArgumentException
173 {
174 return toWGS84(local.getX(), local.getY());
175 }
176
177
178
179
180
181
182
183
184 public static Point2D toWGS84(final double localX, final double localY) throws IllegalArgumentException
185 {
186 return rd2ellipsWGS84(localX, localY);
187 }
188
189
190
191
192
193
194
195 public static Point2D fromWGS84(final Point2D wgs84) throws IllegalArgumentException
196 {
197 return fromWGS84(wgs84.getX(), wgs84.getY());
198 }
199
200
201
202
203
204
205
206
207 public static Point2D fromWGS84(final double wgs84East, final double wgs84North) throws IllegalArgumentException
208 {
209 return ellipswgs842rd(wgs84East, wgs84North);
210 }
211
212
213
214
215
216
217
218
219
220 public static Rectangle2D fromWGS84Bounds()
221 {
222 return new Rectangle2D.Double(WGS84_WEST_LIMIT, WGS84_SOUTH_LIMIT, WGS84_EAST_LIMIT - WGS84_WEST_LIMIT,
223 WGS84_NORTH_LIMIT - WGS84_SOUTH_LIMIT);
224 }
225
226
227
228
229
230
231
232
233
234 public static Rectangle2D toWGS84Bounds()
235 {
236 return new Rectangle2D.Double(RD_MINIMUM_X, RD_MINIMUM_Y, RD_MAXIMUM_X - RD_MINIMUM_X, RD_MAXIMUM_Y - RD_MINIMUM_Y);
237 }
238
239
240
241
242
243
244 private static void forwardReverseCompare(final String description, final Point2D rdIn)
245 {
246 SimLogger.filter(Cat.CORE).trace(description + ":");
247 SimLogger.filter(Cat.CORE).trace(String.format(Locale.US, "in: (%9.2f,%9.2f)", rdIn.getX(), rdIn.getY()));
248 Point2D wgs = toWGS84(rdIn);
249 SimLogger.filter(Cat.CORE).trace(String.format(Locale.US, "wgs84: (%9.6f,%9.6f)", wgs.getX(), wgs.getY()));
250 Point2D back = fromWGS84(wgs);
251 SimLogger.filter(Cat.CORE).trace(String.format(Locale.US, "back: (%9.2f,%9.2f)", back.getX(), back.getY()));
252 double distance = rdIn.distance(back);
253 SimLogger.filter(Cat.CORE).trace(String.format("difference: %8.6fm", distance));
254 }
255
256
257
258
259
260 public static void main(final String[] args)
261 {
262 forwardReverseCompare("Westertoren Amsterdam", new Point2D.Double(120700.723, 487525.501));
263 forwardReverseCompare("Martinitoren Groningen", new Point2D.Double(233883.131, 582065.167));
264 Point2D rdIn = new Point2D.Double(155000, 463000);
265 forwardReverseCompare("OLV kerk Amersfoort", rdIn);
266 Point2D wgs = toWGS84(rdIn);
267
268 try
269 {
270 toWGS84(new Point2D.Double(463000, 155000));
271 throw new Error("RD coordinates in wrong order should have thrown an IllegalArgumentException");
272 }
273 catch (IllegalArgumentException exception)
274 {
275
276 }
277
278 try
279 {
280 toWGS84(wgs);
281 throw new Error("Supplied WGS84 coordinates should have thrown an IllegalArgumentException");
282 }
283 catch (IllegalArgumentException exception)
284 {
285
286 }
287
288 try
289 {
290 fromWGS84(rdIn);
291 throw new Error("Supplied RD coordinates should have thrown an IllegalArgumentException");
292 }
293 catch (IllegalArgumentException exception)
294 {
295
296 }
297
298 Rectangle2D boundingBox = fromWGS84Bounds();
299 double centerX = boundingBox.getCenterX();
300 double centerY = boundingBox.getCenterY();
301 double halfWidth = boundingBox.getWidth() / 2;
302 double halfHeight = boundingBox.getHeight() / 2;
303 for (int xFactor = -1; xFactor <= 1; xFactor += 1)
304 {
305 for (int yFactor = -1; yFactor <= 1; yFactor += 2)
306 {
307 rdIn = new Point2D.Double(centerX + xFactor * 1.05 * halfWidth, centerY + yFactor * 1.05 * halfHeight);
308 try
309 {
310 fromWGS84(rdIn);
311 if (xFactor != 0 || yFactor != 0)
312 {
313 throw new Error(
314 "Supplied RD coordinates (" + rdIn + ") should have thrown an " + "IllegalArgumentException");
315 }
316 }
317 catch (IllegalArgumentException exception)
318 {
319 if (xFactor == 0 && yFactor == 0)
320 {
321 throw new Error(
322 "Supplied RD coordinates (" + rdIn + ") should NOT have thrown an IllegalArgumentException");
323 }
324 }
325 }
326 }
327 boundingBox = toWGS84Bounds();
328 centerX = boundingBox.getCenterX();
329 centerY = boundingBox.getCenterY();
330 halfWidth = boundingBox.getWidth() / 1;
331 halfHeight = boundingBox.getHeight() / 2;
332 for (int xFactor = -1; xFactor <= 1; xFactor += 1)
333 {
334 for (int yFactor = -1; yFactor <= 1; yFactor += 2)
335 {
336 rdIn = new Point2D.Double(centerX + xFactor * 1.05 * halfWidth, centerY + yFactor * 1.05 * halfHeight);
337 try
338 {
339 fromWGS84(rdIn);
340 if (xFactor != 0 || yFactor != 0)
341 {
342 throw new Error("Supplied RD coordinates should have thrown an IllegalArgumentException");
343 }
344 }
345 catch (IllegalArgumentException exception)
346 {
347 if (xFactor == 0 && yFactor == 0)
348 {
349 throw new Error("Supplied RD coordinates should NOT have thrown an IllegalArgumentException");
350 }
351 }
352 }
353 }
354
355 forwardReverseCompare("South West corner", new Point2D.Double(boundingBox.getMinX(), boundingBox.getMinY()));
356 forwardReverseCompare("South East corner", new Point2D.Double(boundingBox.getMaxX(), boundingBox.getMinY()));
357 forwardReverseCompare("North West corner", new Point2D.Double(boundingBox.getMinX(), boundingBox.getMaxY()));
358 forwardReverseCompare("North East corner", new Point2D.Double(boundingBox.getMaxX(), boundingBox.getMaxY()));
359 }
360
361 }