net.sf.geographiclib.Gnomonic Maven / Gradle / Ivy
Show all versions of GeographicLib-Java Show documentation
/**
* Implementation of the net.sf.geographiclib.Gnomonic class
*
* Copyright (c) BMW Car IT GmbH (2014-2016)
* and licensed under the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/
**********************************************************************/
package net.sf.geographiclib;
/**
* Gnomonic projection.
*
* Note: Gnomonic.java has been ported to Java from its C++ equivalent
* Gnomonic.cpp, authored by C. F. F. Karney and licensed under MIT/X11
* license. The following documentation is mostly the same as for its C++
* equivalent, but has been adopted to apply to this Java implementation.
*
* Gnomonic projection centered at an arbitrary position C on the
* ellipsoid. This projection is derived in Section 8 of
*
* -
* C. F. F. Karney,
* Algorithms for geodesics, J. Geodesy 87, 43–55 (2013);
* DOI:
* 10.1007/s00190-012-0578-z; addenda:
*
* geod-addenda.html.
*
*
*
* The gnomonic projection of a point P on the ellipsoid is defined as
* follows: compute the geodesic line from C to P; compute the
* reduced length m12, geodesic scale M12, and ρ =
* m12/M12; finally, this gives the coordinates x and
* y of P in gnomonic projection with x = ρ sin
* azi1; y = ρ cos azi1, where azi1 is the
* azimuth of the geodesic at C. 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
* azi of the geodesic at P and reciprocal scale
* rk in the azimuthal direction. The scale in the radial
* direction is 1/rk2.
*
* For a sphere, ρ reduces to a tan(s12/a), where
* s12 is the length of the geodesic from C to P, 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 l. Let T be the point on
* the geodesic (extended if necessary) closest to C, the center of the
* projection, and t, be the distance CT. To 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
*
* (K(T) - K(C)) l2 t / 32.
*
*
* where K 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 r of
* C. For a given r, the deviation is maximum when the latitude
* of C is 45°, when endpoints are a distance r away, and
* when their azimuths from the center are ± 45° or ±
* 135°. To lowest order in r and the flattening f, the
* deviation is f (r/2a)3 r.
*
* CAUTION: 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:
*
* -
* Declare that the projection is undefined for an ellipsoid.
*
* -
* 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
* by independently by Bowring and Williams in 1997.
*
* -
* 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 sections 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–274 (1963).
*
* -
* 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.
*
*
*
* Example of use:
*
*
* // Example of using the Gnomonic.java class
* import net.sf.geographiclib.Geodesic;
* import net.sf.geographiclib.Gnomonic;
* import net.sf.geographiclib.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);
* }
* }
* }
*
*/
public class Gnomonic {
private static final double eps_ = 0.01 * Math.sqrt(GeoMath.epsilon);
private static final int numit_ = 10;
private Geodesic _earth;
private double _a, _f;
/**
* Constructor for Gnomonic.
*
* @param earth the {@link Geodesic} object to use for geodesic
* calculations.
*/
public Gnomonic(Geodesic earth) {
_earth = earth;
_a = _earth.MajorRadius();
_f = _earth.Flattening();
}
/**
* Forward projection, from geographic to gnomonic.
*
* @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:
* lat0, lon0, lat, lon, x, y,
* azi, rk.
*
* lat0 and lat should be in the range [−90°,
* 90°] and lon0 and lon should be in the range
* [−540°, 540°). The scale of the projection is
* 1/rk2 in the "radial" direction, azi clockwise
* from true north, and is 1/rk in the direction perpendicular to
* this. If the point lies "over the horizon", i.e., if rk ≤ 0,
* then NaNs are returned for x and y (the correct values are
* returned for azi and rk). A call to Forward followed by a
* call to Reverse will return the original (lat, lon) (to
* within roundoff) provided the point in not over the horizon.
*/
public GnomonicData Forward(double lat0, double lon0, double lat, double lon)
{
GeodesicData inv =
_earth.Inverse(lat0, lon0, lat, lon, GeodesicMask.AZIMUTH
| GeodesicMask.GEODESICSCALE | GeodesicMask.REDUCEDLENGTH);
GnomonicData fwd =
new GnomonicData(lat0, lon0, lat, lon, Double.NaN, Double.NaN,
inv.azi2, inv.M12);
if (inv.M12 > 0) {
double rho = inv.m12 / inv.M12;
Pair p = GeoMath.sincosd(inv.azi1);
fwd.x = rho * p.first;
fwd.y = rho * p.second;
}
return fwd;
}
/**
* Reverse projection, from gnomonic to geographic.
*
* @param lat0 latitude of center point of projection (degrees).
* @param lon0 longitude of center point of projection (degrees).
* @param x easting of point (meters).
* @param y northing of point (meters).
* @return {@link GnomonicData} object with the following fields:
* lat0, lon0, lat, lon, x, y,
* azi, rk.
*
* lat0 should be in the range [−90°, 90°] and
* lon0 should be in the range [−540°, 540°).
* lat will be in the range [−90°, 90°] and lon
* will be in the range [−180°, 180°). The scale of the
* projection is 1/rk2 in the "radial" direction,
* azi clockwise from true north, and is 1/rk in the direction
* perpendicular to this. Even though all inputs should return a valid
* lat and lon, it's possible that the procedure fails to
* converge for very large x or y; 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 (x, y) (to
* roundoff).
*/
public GnomonicData Reverse(double lat0, double lon0, double x, double y) {
GnomonicData rev =
new GnomonicData(lat0, lon0, Double.NaN, Double.NaN, x, y, Double.NaN,
Double.NaN);
double azi0 = GeoMath.atan2d(x, y);
double rho = Math.hypot(x, y);
double s = _a * Math.atan(rho / _a);
boolean little = rho <= _a;
if (!little)
rho = 1 / rho;
GeodesicLine line =
_earth.Line(lat0, lon0, azi0, GeodesicMask.LATITUDE
| GeodesicMask.LONGITUDE | GeodesicMask.AZIMUTH
| GeodesicMask.DISTANCE_IN | GeodesicMask.REDUCEDLENGTH
| GeodesicMask.GEODESICSCALE);
int count = numit_, trip = 0;
GeodesicData pos = null;
while (count-- > 0) {
pos =
line.Position(s, GeodesicMask.LONGITUDE | GeodesicMask.LATITUDE
| GeodesicMask.AZIMUTH | GeodesicMask.DISTANCE_IN
| GeodesicMask.REDUCEDLENGTH
| GeodesicMask.GEODESICSCALE);
if (trip > 0)
break;
double ds =
little ? ((pos.m12 / pos.M12) - rho) * pos.M12 * pos.M12
: (rho - (pos.M12 / pos.m12)) * pos.m12 * pos.m12;
s -= ds;
if (Math.abs(ds) <= eps_ * _a)
trip++;
}
if (trip == 0)
return rev;
rev.lat = pos.lat2;
rev.lon = pos.lon2;
rev.azi = pos.azi2;
rev.rk = pos.M12;
return rev;
}
/**
* @return a the equatorial radius of the ellipsoid (meters). This is
* the value inherited from the Geodesic object used in the constructor.
**********************************************************************/
public double MajorRadius() { return _a; }
/**
* @return f the flattening of the ellipsoid. This is
* the value inherited from the Geodesic object used in the constructor.
**********************************************************************/
public double Flattening() { return _f; }
}