Gnomonic.java
/*
* Copyright (C) 2018 Alberto Irurueta Carro (alberto@irurueta.com)
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package com.irurueta.navigation.geodesic;
/**
* Gnomonic projection.
* Gnomonic projection centered at an arbitrary position <i>C</i> on the ellipsoid. This projection
* is derived in Section 8 of
* <ul>
* <li>
* C. F. F. Karney, <a href="https://doi.org/10.1007/s00190-012-0578-z">Algorithms for
* geodesics</a>, J. Geodesy <b>87</b>, 43–55 (2013);
* DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">10.1007/s00190-012-0578-z</a>;
* </li>
* </ul>
* The gnomonic projection of a point <i>P</i> on the ellipsoid is defined as follows: compute the
* geodesic line from <i>C</i> to <i>P</i>; compute the reduced length <i>m12</i>, geodesic scale
* <i>M12</i>, and ρ = <i>m12</i>/<i>M12</i>; finally, this gives the coordinates <i>x</i> and
* <i>y</i> of <i>P</i> in gnomonic projection with <i>x</i> = ρ sin <i>azi1</i>; <i>y</i> = ρ
* cos <i>azi1</i>, where <i>azi1</i> is the azimuth of the geodesic at <i>C</i>. The method
* {@link Gnomonic#forward(double, double, double, double)} performs the forward projection and
* {@link Gnomonic#reverse(double, double, double, double)} is the inverse of the projection. The
* methods also return the azimuth <i>azi</i> of the geodesic at <i>P</i> and reciprocal scale
* <i>rk</i> in the azimuthal direction. The scale in the radial direction is 1/<i>rk</i><sup>2</sup>.
* For a sphere, ρ reduces to <i>a</i> tan (<i>s12</i>/<i>a</i>), where <i>s12</i> is the length
* of the geodesic from <i>C</i> to <i>P</i>, and the gnomonic projection has the property that all
* geodesics appear as straight lines. For an ellipsoid, this property holds only for geodesics
* interesting the centers. However, geodesic segments close to the center are approximately straight.
* Consider a geodesic segment of length <i>1</i>. Let <i>T</i> be the point on the geodesic
* (extended if necessary) closest to <i>C</i>, the center of the projection, and <i>t</i>, be the
* distance <i>CT</i>. To the lowest order, the maximum deviation (as a true distance) of the corresponding
* gnomonic line segment (i.e., with the same end points) from the geodesic is
* (<i>K</i>(<i>T</i> − <i>K</i>(<i>C</i>))
* <i>l</i><sup>2</sup> <i>t</i> / 32.
* where <i>K</i> is the Gaussian curvature.
* This result applies for any surface. For an ellipsoid of revolution, consider all geodesics whose end
* points are within a distance <i>r</i> of <i>C</i>. For a given <i>r</i>, the deviation is maximum
* when the latitude of <i>C</i> is 46°, when endpoints are a distance <i>r</i> away, and when
* their azimuths from the center are ± 45° or ± 135°. To the lowest order in
* <i>r</i> and the flattening <i>f</i>, the deviation is <i>f</i> (<i>r</i>/2<i>a</i>)<sup>3</sup>
* <i>r</i>.
* <b>CAUTION:</b> The definition of this projection for a sphere is standard. However, there is no
* standard for how it should be extended to an ellipsoid.
* The choices are:
* <ul>
* <li>
* Declare that the projection is undefined for an ellipsoid.
* </li>
* <li>
* Project to a tangent plane from the center of the ellipsoid. This causes great ellipses to
* appear as straight lines in the projection; i.e., it generalizes the spherical great circle
* to a great ellipse. This was proposed independently by Bowring and Williams in 1997.
* </li>
* <li>
* Project to the conformal sphere with the constant of integration chosen so that the values
* of the latitude match for the center point and perform a central projection onto the plane
* tangent to the conformal sphere at the center point. This causes normal section through the
* center point to appear as straight lines in the projection; i.e., it generalizes the
* spherical great circle to a normal section. This was proposed by I. G. Letoval'tsev,
* Generalization of the gnomonic projection for a spheroid and the principal geodetic
* problems involved in the alignment of surface routes, Geodesy and Aerophotography(5),
* 271–275 (1963)
* </li>
* <li>
* The projection given here. This causes geodesics close to the center point to appear as
* straight lines in the projection; i.e., it generalizes the spherical great circle to a
* geodesic.
* </li>
* </ul>
* Example of use:
* <pre>
* // Example of using the Gnomonic.java class
* import com.irurueta.navigation.geodesic.Geodesic;
* import com.irurueta.navigation.geodesic.Gnomonic;
* import com.irurueta.navigation.geodesic.GnomonicData;
* public class ExampleGnomonic {
* public static void main(String[] args) {
* Geodesic geod = Geodesic.WGS84;
* double lat0 = 48 + 50 / 60.0, lon0 = 2 + 20 / 60.0; // Paris
* Gnomonic gnom = new Gnomonic(geod);
* {
* // Sample forward calculation
* double lat = 50.9, lon = 1.8; // Calais
* GnomonicData proj = gnom.Forward(lat0, lon0, lat, lon);
* System.out.println(proj.x + " " + proj.y);
* }
* {
* // Sample reverse calculation
* double x = -38e3, y = 230e3;
* GnomonicData proj = gnom.Reverse(lat0, lon0, x, y);
* System.out.println(proj.lat + " " + proj.lon);
* }
* }
* }
* </pre>
*/
public class Gnomonic {
private static final double EPS = 0.01 * Math.sqrt(GeoMath.EPSILON);
private static final int NUMIT = 10;
/**
* Earth geodesic.
*/
private final Geodesic earth;
/**
* Major equatorial Earth radius.
*/
private final double a;
/**
* Earth flattening.
*/
private final double f;
/**
* Constructor for Gnomonic.
*
* @param earth the {@link Geodesic} object to use for geodesic calculations.
*/
public Gnomonic(final Geodesic earth) {
this.earth = earth;
a = this.earth.getMajorRadius();
f = this.earth.getFlattening();
}
/**
* Forward projection, from geographic to gnomonic.
* <i>lat0</i> and <i>lat</i> should be in the range [−90°, 90°] and <i>lon0</i>
* and <i>lon</i> should be in the range [−540°, 540°). The scale of the projection
* is 1/<i>rk<sup>2</sup></i> in the "radial" direction, <i>azi</i> clockwise from true north,
* and is 1/<i>rk</i> in the direction perpendicular to this. If the point lies "over the
* horizon", i.e., if <i>rk</i> ≤ 0, then NaNs are returned for <i>x</i> and <i>y</i> (the
* correct values are returned for <i>azi</i> and <i>rk</i>). A call to forward followed by a
* call to reverse will return the original (<i>lat</i>, <i>lon</i>) (to within roundoff)
* provided the point in not over the horizon.
*
* @param lat0 latitude of center point of projection (degrees).
* @param lon0 longitude of center point of projection (degrees).
* @param lat latitude of point (degrees).
* @param lon longitude of point (degrees).
* @return {@link GnomonicData} object with the following fields:
* <i>lat0</i>, <i>lon0</i>, <i>lat</i>, <i>lon</i>, <i>x</i>, <i>y</i>, <i>azi</i>, <i>rk</i>.
*/
public GnomonicData forward(final double lat0, final double lon0, final double lat, final double lon) {
final var inv = earth.inverse(lat0, lon0, lat, lon,
GeodesicMask.AZIMUTH | GeodesicMask.GEODESIC_SCALE | GeodesicMask.REDUCED_LENGTH);
final var fwd = new GnomonicData(lat0, lon0, lat, lon, Double.NaN, Double.NaN, inv.getAzi2(),
inv.getScaleM12());
if (inv.getScaleM12() > 0) {
final var rho = inv.getM12() / inv.getScaleM12();
final var p = GeoMath.sincosd(inv.getAzi1());
fwd.setX(rho * p.getFirst());
fwd.setY(rho * p.getSecond());
}
return fwd;
}
/**
* Reverse projection, from gnomonic to geographic.
* <i>lat</i> will be in the range [−90°, 90°] and <i>lon</i> will be in the
* range (−180°, 180°]. The scale of the projection is 1/<i>rk<sup>2</sup></i>
* in the "radial" direction, <i>azi</i> clockwise from true north, and is 1/<i>rk</i> in the
* direction perpendicular to this. Even though all inputs should return a valid <i>lat</i>
* and <i>lon</i>, it's possible that the procedure fails to converge for very large <i>x</i>
* or <i>y</i>; in this case NaNs are returned for very large <i>x</i> or <i>y</i>; in this
* case NaNs are returned for all the output arguments. A call to reverse followed by a call
* to forward will return the original (<i>x</i>, <i>y</i>) (to round-off).
*
* @param lat0 latitude of center point of projection (degrees). <i>lat0</i> should be in the
* range [−90°, 90°]
* @param lon0 longitude of center point of projection (degrees). <i>lon0</i> should be in the
* range [−540°, 540°).
* @param x easting of point (meters).
* @param y northing of point (meters).
* @return {@link GnomonicData} object with the following fields:
* <i>lat0</i>, <i>lon0</i>, <i>lat</i>, <i>lon</i>, <i>x</i>, <i>y</i>,
* <i>azi</i>, <i>rk</i>.
*/
public GnomonicData reverse(final double lat0, final double lon0, final double x, final double y) {
final var rev = new GnomonicData(lat0, lon0, Double.NaN, Double.NaN, x, y, Double.NaN, Double.NaN);
//noinspection all
final var azi0 = GeoMath.atan2d(x, y);
var rho = Math.hypot(x, y);
var s = a * Math.atan(rho / a);
final var little = rho <= a;
if (!little) {
rho = 1 / rho;
}
final var line = earth.line(lat0, lon0, azi0, GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE
| GeodesicMask.AZIMUTH | GeodesicMask.DISTANCE_IN | GeodesicMask.REDUCED_LENGTH
| GeodesicMask.GEODESIC_SCALE);
var count = NUMIT;
var trip = 0;
GeodesicData pos = null;
while (count-- > 0) {
pos = line.position(s, GeodesicMask.LONGITUDE | GeodesicMask.LATITUDE | GeodesicMask.AZIMUTH
| GeodesicMask.DISTANCE_IN | GeodesicMask.REDUCED_LENGTH | GeodesicMask.GEODESIC_SCALE);
if (trip > 0) {
break;
}
final var ds = little
? ((pos.getM12() / pos.getScaleM12()) - rho) * pos.getScaleM12() * pos.getScaleM12()
: (rho - (pos.getScaleM12() / pos.getM12())) * pos.getM12() * pos.getM12();
s -= ds;
if (Math.abs(ds) <= EPS * a) {
trip++;
}
}
if (trip == 0) {
return rev;
}
rev.setLat(pos.getLat2());
rev.setLon(pos.getLon2());
rev.setAzi(pos.getAzi2());
rev.setRk(pos.getScaleM12());
return rev;
}
/**
* Gets the equatorial radius of the ellipsoid (meters). This is the value inherited from the
* Geodesic object used in the constructor.
*
* @return <i>a</i> the equatorial radius of the ellipsoid (meters).
*/
public double getMajorRadius() {
return a;
}
/**
* Gets the flattening of the ellipsoid.
*
* @return <i>f</i> the flattening of the ellipsoid. This is the value inherited from the
* Geodesic object used in the constructor.
*/
public double getFlattening() {
return f;
}
}