001 /*
002 * Java GPX Library (jpx-3.1.0).
003 * Copyright (c) 2016-2023 Franz Wilhelmstötter
004 *
005 * Licensed under the Apache License, Version 2.0 (the "License");
006 * you may not use this file except in compliance with the License.
007 * You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 *
017 * Author:
018 * Franz Wilhelmstötter (franz.wilhelmstoetter@gmail.com)
019 */
020 package io.jenetics.jpx.geom;
021
022 import static java.lang.Math.abs;
023 import static java.lang.Math.asin;
024 import static java.lang.Math.atan;
025 import static java.lang.Math.atan2;
026 import static java.lang.Math.cos;
027 import static java.lang.Math.sin;
028 import static java.lang.Math.sqrt;
029 import static java.lang.Math.tan;
030 import static java.lang.String.format;
031 import static java.util.Objects.requireNonNull;
032 import static io.jenetics.jpx.geom.MathUtils.equal;
033
034 import java.util.stream.Collector;
035
036 import io.jenetics.jpx.Length;
037 import io.jenetics.jpx.Length.Unit;
038 import io.jenetics.jpx.Point;
039
040 /**
041 * Implementation of <em>geodetic</em> functions.
042 *
043 * @see <a href="https://en.wikipedia.org/wiki/Geoid">Wikipedia: Geoid</a>
044 * @see Ellipsoid
045 *
046 * @author <a href="mailto:franz.wilhelmstoetter@gmail.com">Franz Wilhelmstötter</a>
047 * @version 1.0
048 * @since 1.0
049 */
050 public final class Geoid {
051
052 private static final int EPSILON_ULP = 10_000;
053
054 /**
055 * {@link Geoid} using of the <em>World Geodetic System: WGS 84</em>
056 *
057 * @see <a href="https://en.wikipedia.org/wiki/World_Geodetic_System#A_new_World_Geodetic_System:_WGS_84">
058 * WGS-84</a>
059 */
060 public static final Geoid WGS84 = of(Ellipsoid.WGS84);
061
062 /**
063 * {@link Geoid} using the <em>International Earth Rotation and Reference
064 * Systems Service (1989)</em>
065 *
066 * @see <a href="https://en.wikipedia.org/wiki/IERS">IERS-89</a>
067 */
068 public static final Geoid IERS_1989 = of(Ellipsoid.IERS_1989);
069
070 /**
071 * {@link Geoid} using the <em>International Earth Rotation and Reference
072 * Systems Service (2003)</em>
073 *
074 * @see <a href="https://en.wikipedia.org/wiki/IERS">IERS-89</a>
075 */
076 public static final Geoid IERS_2003 = of(Ellipsoid.IERS_2003);
077
078 /**
079 * {@link Geoid} using the {@link Ellipsoid#DEFAULT} ellipsoid.
080 */
081 public static final Geoid DEFAULT = of(Ellipsoid.DEFAULT);
082
083 private final Ellipsoid _ellipsoid;
084
085 // Minor semi-axes of the ellipsoid.
086 private final double B;
087
088 private final double AABBBB;
089
090 // Flattening (A - B)/A
091 private final double F;
092
093 // The maximal iteration of the 'distance'
094 private static final int DISTANCE_ITERATION_MAX = 1000;
095
096 // The epsilon of the result, when to stop iteration.
097 private static final double DISTANCE_ITERATION_EPSILON = 1E-12;
098
099 /**
100 * Create a new {@code Geoid} object with the given ellipsoid.
101 *
102 * @param ellipsoid the ellipsoid used by the geoid
103 * @throws NullPointerException if the given {@code ellipsoid} is {@code null}
104 */
105 private Geoid(final Ellipsoid ellipsoid) {
106 _ellipsoid = requireNonNull(ellipsoid);
107
108 final double a = ellipsoid.A();
109 final double aa = a*a;
110
111 B = ellipsoid.B();
112 final double bb = B*B;
113
114 AABBBB = (aa - bb)/bb;
115 F = 1.0/ellipsoid.F();
116 }
117
118 /**
119 * Return the ellipsoid the {@code Geom} object is using.
120 *
121 * @return the ellipsoid the {@code Geom} object is using
122 */
123 public Ellipsoid ellipsoid() {
124 return _ellipsoid;
125 }
126
127 /**
128 * Calculate the distance between points on an ellipsoidal earth model. This
129 * method will throw an {@link ArithmeticException} if the algorithm doesn't
130 * converge while calculating the distance, which is the case for a point
131 * and its (near) antidote.
132 *
133 * @see <a href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">DIRECT AND
134 * INVERSE SOLUTIONS OF GEODESICS 0 THE ELLIPSOID
135 * WITH APPLICATION OF NESTED EQUATIONS</a>
136 * @see <a href="http://www.movable-type.co.uk/scripts/latlong-vincenty.html">
137 * Vincenty solutions of geodesics on the ellipsoid</a>
138 *
139 * @param start the start point
140 * @param end the end point
141 * @return the distance between {@code start} and {@code end} in meters
142 * @throws NullPointerException if one of the points is {@code null}
143 * @throws ArithmeticException if the algorithm used for calculating the
144 * distance between {@code start} and {@code end} didn't converge,
145 * which is the case for a point and its (near) antidote.
146 */
147 public Length distance(final Point start, final Point end) {
148 final double lat1 = start.getLatitude().toRadians();
149 final double lon1 = start.getLongitude().toRadians();
150 final double lat2 = end.getLatitude().toRadians();
151 final double lon2 = end.getLongitude().toRadians();
152
153 final double omega = lon2 - lon1;
154
155 final double tanphi1 = tan(lat1);
156 final double tanU1 = (1.0 - F)*tanphi1;
157 final double U1 = atan(tanU1);
158 final double sinU1 = sin(U1);
159 final double cosU1 = cos(U1);
160
161 final double tanphi2 = tan(lat2);
162 final double tanU2 = (1.0 - F)*tanphi2;
163 final double U2 = atan(tanU2);
164 final double sinU2 = sin(U2);
165 final double cosU2 = cos(U2);
166
167 final double sinU1sinU2 = sinU1*sinU2;
168 final double cosU1sinU2 = cosU1*sinU2;
169 final double sinU1cosU2 = sinU1*cosU2;
170 final double cosU1cosU2 = cosU1*cosU2;
171
172 // Eq. 13
173 double lambda = omega;
174
175 // Intermediates we'll need to compute distance 's'
176 double a;
177 double b;
178 double sigma;
179 double deltasigma;
180 double lambda0;
181
182 int iteration = 0;
183 do {
184 lambda0 = lambda;
185
186 final double sinlambda = sin(lambda);
187 final double coslambda = cos(lambda);
188
189 // Eq. 14
190 final double sin2sigma =
191 (cosU2*sinlambda*cosU2*sinlambda) +
192 (cosU1sinU2 - sinU1cosU2*coslambda)*
193 (cosU1sinU2 - sinU1cosU2*coslambda);
194 final double sinsigma = sqrt(sin2sigma);
195
196 // Eq. 15
197 final double cossigma = sinU1sinU2 + (cosU1cosU2*coslambda);
198
199 // Eq. 16
200 sigma = atan2(sinsigma, cossigma);
201
202 // Eq. 17 Careful! sin2sigma might be almost 0!
203 final double sinalpha = equal(sin2sigma, 0.0, EPSILON_ULP)
204 ? 0.0
205 : cosU1cosU2*sinlambda/sinsigma;
206 final double alpha = asin(sinalpha);
207 final double cosalpha = cos(alpha);
208 double cos2alpha = cosalpha*cosalpha;
209
210 // Eq. 18 Careful! cos2alpha might be almost 0!
211 final double cos2sigmam = equal(cos2alpha, 0.0, EPSILON_ULP)
212 ? 0.0
213 : cossigma - 2*sinU1sinU2/cos2alpha;
214 final double u2 = cos2alpha*AABBBB;
215
216 final double cos2sigmam2 = cos2sigmam*cos2sigmam;
217
218 // Eq. 3
219 a = 1.0 + u2/16384*(4096 + u2*(-768 + u2*(320 - 175*u2)));
220
221 // Eq. 4
222 b = u2/1024*(256 + u2*(-128 + u2*(74 - 47*u2)));
223
224 // Eq. 6
225 deltasigma = b*sinsigma*(cos2sigmam +
226 b/4*(cossigma*(-1 + 2 * cos2sigmam2) -
227 b/6*cos2sigmam*(-3 + 4*sin2sigma)*(-3 + 4*cos2sigmam2)));
228
229 // Eq. 10
230 final double C = F/16*cos2alpha*(4 + F*(4 - 3*cos2alpha));
231
232 // Eq. 11
233 lambda = omega + (1 - C)*F*sinalpha*
234 (sigma + C*sinsigma*(cos2sigmam +
235 C*cossigma*(-1 + 2*cos2sigmam2)));
236
237 } while (iteration++ < DISTANCE_ITERATION_MAX &&
238 (abs((lambda - lambda0)/lambda) > DISTANCE_ITERATION_EPSILON));
239
240 if (iteration >= DISTANCE_ITERATION_MAX) {
241 throw new ArithmeticException(format(
242 "Calculating distance between %s and %s didn't converge.",
243 start, end
244 ));
245 }
246
247 // Eq. 19
248 final double s = B*a*(sigma - deltasigma);
249
250 return Length.of(s, Unit.METER);
251 }
252
253 /**
254 * Return a collector which calculates the length of the (open) path which
255 * is defined by the {@code Point} stream.
256 *
257 * <pre>{@code
258 * final Length length = gpx.tracks()
259 * .flatMap(Track::segments)
260 * .flatMap(TrackSegment::points)
261 * .collect(Geoid.WGSC_84.toPathLength());
262 * }</pre>
263 *
264 * <b>The returned {@code Collector} doesn't work for <em>parallel</em>
265 * stream. Using it for a <em>parallel</em> point stream will throw an
266 * {@link UnsupportedOperationException} at runtime.</b>
267 *
268 * @see #toTourLength()
269 *
270 * @return a new path length collector
271 */
272 public Collector<Point, ?, Length> toPathLength() {
273 return Collector.of(
274 () -> new LengthCollector(this),
275 LengthCollector::add,
276 LengthCollector::combine,
277 LengthCollector::pathLength
278 );
279 }
280
281 /**
282 * Return a collector which calculates the length of the (closed) tour which
283 * is defined by the {@code Point} stream. The <em>tour</em> length
284 * additionally adds the distance of the last point back to the first point.
285 *
286 * <pre>{@code
287 * final Length length = gpx.tracks()
288 * .flatMap(Track::segments)
289 * .flatMap(TrackSegment::points)
290 * .collect(Geoid.WGSC_84.toTourLength());
291 * }</pre>
292 *
293 * <b>The returned {@code Collector} doesn't work for <em>parallel</em>
294 * stream. Using it for a <em>parallel</em> point stream will throw an
295 * {@link UnsupportedOperationException} at runtime.</b>
296 *
297 * @see #toPathLength()
298 *
299 * @return a new path length collector
300 */
301 public Collector<Point, ?, Length> toTourLength() {
302 return Collector.of(
303 () -> new LengthCollector(this),
304 LengthCollector::add,
305 LengthCollector::combine,
306 LengthCollector::tourLength
307 );
308 }
309
310
311
312 /**
313 * Create a new {@code Geoid} object with the given ellipsoid.
314 *
315 * @param ellipsoid the ellipsoid used by the geoid
316 * @return a new {@code Geoid} object with the given ellipsoid
317 * @throws NullPointerException if the given {@code ellipsoid} is {@code null}
318 */
319 public static Geoid of(final Ellipsoid ellipsoid) {
320 return new Geoid(ellipsoid);
321 }
322
323 }
|