1 package org.opentrafficsim.core.gis;
2
3 import java.awt.geom.Point2D;
4 import java.awt.geom.Rectangle2D;
5
6 /**
7 * Convert geographical coordinates between WGS84 and the Dutch RD (Rijksdriehoek) system. <br>
8 * Specific MathTransform for WGS84 (EPSG:4326) to RD_new (EPSG:28992) conversions. Code based on C code by Peter Knoppers as
9 * applied <a href="http://www.regiolab-delft.nl/?q=node/36">here</a>, which is based on
10 * <a href="http://home.solcon.nl/pvanmanen/Download/Transformatieformules.pdf">this</a> paper.
11 * <p>
12 * Copyright (c) ~2000-2019 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
13 * BSD-style license. See <a href="http://opentrafficsim.org/docs/current/license.html">OpenTrafficSim License</a>.
14 * <p>
15 * @version $Revision$, $LastChangedDate$, by $Author$, initial version Aug 4, 2016 <br>
16 * @author <a href="http://www.tbm.tudelft.nl/averbraeck">Alexander Verbraeck</a>
17 * @author <a href="http://www.tudelft.nl/pknoppers">Peter Knoppers</a>
18 * @author Gert-Jan Stolk
19 **/
20 public final class TransformWGS84DutchRDNew
21 {
22
23 /** Western boundary of the Dutch RD system. */
24 private static final double WGS84_WEST_LIMIT = 3.2;
25
26 /** Eastern boundary of the Dutch RD system. */
27 private static final double WGS84_EAST_LIMIT = 7.3;
28
29 /** Northern boundary of the Dutch RD system. */
30 private static final double WGS84_SOUTH_LIMIT = 50.6;
31
32 /** Southern boundary of the Dutch RD system. */
33 private static final double WGS84_NORTH_LIMIT = 53.7;
34
35 /** Western boundary of the Dutch RD system. */
36 private static final double RD_MINIMUM_X = 11000;
37
38 /** Eastern boundary of the Dutch RD system. */
39 private static final double RD_MAXIMUM_X = 280000;
40
41 /** Southern boundary of the Dutch RD system. */
42 private static final double RD_MINIMUM_Y = 300000;
43
44 /** Northern boundary of the Dutch RD system. */
45 private static final double RD_MAXIMUM_Y = 630000;
46
47 /** This class is a utility class and instances cannot be constructed. */
48 private TransformWGS84DutchRDNew()
49 {
50 // This class is a utility class and instances cannot be constructed.
51 }
52
53 /**
54 * Convert from WGS84 to RD coordinates.
55 * @param wgs84East double; Degrees East of Greenwich
56 * @param wgs84North double; Degrees North of the equator
57 * @return Point2D; equivalent location in the Dutch RD system
58 */
59 private static Point2D.Double ellipswgs842rd(final double wgs84East, final double wgs84North)
60 {
61 if (wgs84North > WGS84_NORTH_LIMIT || wgs84North < WGS84_SOUTH_LIMIT || wgs84East < WGS84_WEST_LIMIT
62 || (wgs84East > WGS84_EAST_LIMIT))
63 {
64 throw new IllegalArgumentException("ellipswgs842rd input out of range (" + wgs84East + ", " + wgs84North + ")");
65 }
66 //@formatter:off
67 @SuppressWarnings("checkstyle:nowhitespacebefore")
68 /** Coefficients for ellipswgs842rd. */
69 final double[][] r =
70 { /* p down, q right */
71 { 155000.00, 190094.945, -0.008, -32.391, 0.0 , },
72 { -0.705, -11832.228, 0.0 , 0.608, 0.0 , },
73 { 0.0 , -114.221, 0.0 , 0.148, 0.0 , },
74 { 0.0 , -2.340, 0.0 , 0.0 , 0.0 , },
75 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , }
76 };
77 @SuppressWarnings("checkstyle:nowhitespacebefore")
78 final double[][] s =
79 { /* p down, q right */
80 { 463000.00 , 0.433, 3638.893, 0.0 , 0.092, },
81 { 309056.544, -0.032, -157.984, 0.0 , -0.054, },
82 { 73.077, 0.0 , -6.439, 0.0 , 0.0 , },
83 { 59.788, 0.0 , 0.0 , 0.0 , 0.0 , },
84 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , }
85 };
86 //@formatter:on
87 double resultX = 0;
88 double resultY = 0;
89 double powNorth = 1;
90 double dNorth = 0.36 * (wgs84North - 52.15517440);
91 double dEast = 0.36 * (wgs84East - 5.38720621);
92
93 for (int p = 0; p < 5; p++)
94 {
95 double powEast = 1;
96 for (int q = 0; q < 5; q++)
97 {
98 resultX += r[p][q] * powEast * powNorth;
99 resultY += s[p][q] * powEast * powNorth;
100 powEast *= dEast;
101 }
102 powNorth *= dNorth;
103 }
104 return new Point2D.Double(resultX, resultY);
105 }
106
107 /**
108 * Convert coordinates from WGS84 to the Dutch RD system.
109 * @param rdX double; X coordinate in the Dutch RD system
110 * @param rdY double; Y coordinate in the Dutch RD system
111 * @return Point2D; equivalent location in the WGS84 system
112 * @throws IllegalArgumentException when the given coordinates are not within the area of the Dutch RD system
113 */
114 private static Point2D rd2ellipsWGS84(final double rdX, final double rdY) throws IllegalArgumentException
115 {
116 if (rdX < RD_MINIMUM_X || rdX > RD_MAXIMUM_X || rdY < RD_MINIMUM_Y || rdY > RD_MAXIMUM_Y)
117 {
118 throw new IllegalArgumentException(
119 "Location (" + rdX + "," + rdY + ") is not within the range " + "of the Dutch RD system");
120 }
121 final double dX = (rdX - 155000) / 100000;
122 final double dY = (rdY - 463000) / 100000;
123 /* Coefficients are for zone 31 (0E .. 6E); roughly west of Apeldoorn */
124 //@formatter:off
125 @SuppressWarnings("checkstyle:nowhitespacebefore")
126 final double[][] k =
127 { /* p down, q right */
128 { 3600 * 52.15517440, 3235.65389, -0.24750, -0.06550, 0.0 , },
129 { -0.00738 , -0.00012, 0.0 , 0.0 , 0.0 , },
130 { -32.58297 , -0.84978, -0.01709, -0.00039, 0.0 , },
131 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , },
132 { 0.00530 , 0.00033, 0.0 , 0.0 , 0.0 , },
133 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , }
134 };
135 @SuppressWarnings("checkstyle:nowhitespacebefore")
136 final double[][] l = { /* p down, q right */
137 { 3600 * 5.38720621, 0.01199, 0.00022, 0.0 , 0.0 , },
138 { 5260.52916 , 105.94684, 2.45656, 0.05594, 0.00128, },
139 { -0.00022 , 0.0 , 0.0 , 0.0 , 0.0 , },
140 { -0.81885 , -0.05607, -0.00256, 0.0 , 0.0 , },
141 { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , },
142 { 0.00026 , 0.0 , 0.0 , 0.0 , 0.0 , }
143 };
144 //@formatter:on
145 double resultNorth = 0;
146 double resultEast = 0;
147 double powX = 1;
148 for (int p = 0; p < 6; p++)
149 {
150 double powY = 1;
151 for (int q = 0; q < 5; q++)
152 {
153 resultNorth += k[p][q] * powX * powY / 3600;
154 resultEast += l[p][q] * powX * powY / 3600;
155 powY *= dY;
156 }
157 powX *= dX;
158 }
159 return new Point2D.Double(resultEast, resultNorth);
160 }
161
162 /**
163 * Convert a coordinate pair in the local system to WGS84 coordinates.
164 * @param local Point2D; coordinates in the local system.
165 * @return Point2D; the equivalent location in degrees in the WGS84 coordinate system
166 * @throws IllegalArgumentException when <cite>local</cite> is not valid in the local system
167 */
168 public static Point2D toWGS84(final Point2D local) throws IllegalArgumentException
169 {
170 return toWGS84(local.getX(), local.getY());
171 }
172
173 /**
174 * Convert a coordinate pair in the local system to WGS84 coordinates.
175 * @param localX double; X-coordinate in the local system
176 * @param localY double; Y-coordinate in the local system
177 * @return Point2D; the equivalent location in degrees in the WGS84 coordinate system
178 * @throws IllegalArgumentException when <cite>localX</cite>, <cite>localY</cite> is not valid in the local system
179 */
180 public static Point2D toWGS84(final double localX, final double localY) throws IllegalArgumentException
181 {
182 return rd2ellipsWGS84(localX, localY);
183 }
184
185 /**
186 * Convert a coordinate pair in WGS84 coordinates to local coordinates.
187 * @param wgs84 Point2D; coordinates in degrees in the WGS84 coordinate system
188 * @return Point2D; the equivalent location in the local coordinate system
189 * @throws IllegalArgumentException when <cite>wgs84</cite> is not valid in the local system
190 */
191 public static Point2D fromWGS84(final Point2D wgs84) throws IllegalArgumentException
192 {
193 return fromWGS84(wgs84.getX(), wgs84.getY());
194 }
195
196 /**
197 * Convert a coordinate pair in WGS84 coordinates to local coordinates.
198 * @param wgs84East double; East coordinate in degrees in the WGS84 system (negative value indicates West)
199 * @param wgs84North double; North coordinate in degrees in the WGS84 system (negative value indicates South)
200 * @return Point2D; the equivalent location in the local coordinate system
201 * @throws IllegalArgumentException when <cite>wgs84</cite> is not valid in the local system
202 */
203 public static Point2D fromWGS84(final double wgs84East, final double wgs84North) throws IllegalArgumentException
204 {
205 return ellipswgs842rd(wgs84East, wgs84North);
206 }
207
208 /**
209 * Report the bounding box for conversion to the local coordinate system.<br>
210 * Conversions from WGS84 to the local coordinate system should fail for locations outside this bounding box. If the valid
211 * range is not adequately described by a rectangular bounding box, conversions for some areas within this bounding box may
212 * also fail (with an IllegalArgumentException). There is no guarantee that the result of a conversion lies within the
213 * bounding box for the reverse conversion.
214 * @return Rectangle2D; bounding box in WGS84 degrees
215 */
216 public static Rectangle2D fromWGS84Bounds()
217 {
218 return new Rectangle2D.Double(WGS84_WEST_LIMIT, WGS84_SOUTH_LIMIT, WGS84_EAST_LIMIT - WGS84_WEST_LIMIT,
219 WGS84_NORTH_LIMIT - WGS84_SOUTH_LIMIT);
220 }
221
222 /**
223 * Report the bounding box for conversions from the local coordinate system. <br>
224 * Conversions from the local coordinate system to WGS84 should fail for locations outside this bounding box. If the valid
225 * range is not adequately described by a rectangular bounding box, conversions for some areas within this bounding box may
226 * also fail (with an IllegalArgumentException). There is no guarantee that the result of a conversion lies within the
227 * bounding box for the reverse conversion.
228 * @return Rectangle2D; bounding box of the local coordinate system
229 */
230 public static Rectangle2D toWGS84Bounds()
231 {
232 return new Rectangle2D.Double(RD_MINIMUM_X, RD_MINIMUM_Y, RD_MAXIMUM_X - RD_MINIMUM_X, RD_MAXIMUM_Y - RD_MINIMUM_Y);
233 }
234
235 }