/*
* Copyright (C) 2011 United States Government as represented by the Administrator of the
* National Aeronautics and Space Administration.
* All Rights Reserved.
*/
package gov.nasa.worldwind.geom;
import gov.nasa.worldwind.util.*;
import java.util.*;
/**
* Represents a point on the two-dimensional surface of a globe. Latitude is the degrees North and ranges between [-90,
* 90], while longitude refers to degrees East, and ranges between (-180, 180].
* <p/>
* Instances of <code>LatLon</code> are immutable.
*
* @author Tom Gaskins
* @version $Id$
*/
public class LatLon
{
public static final LatLon ZERO = new LatLon(Angle.ZERO, Angle.ZERO);
public final Angle latitude;
public final Angle longitude;
/**
* Factor method for obtaining a new <code>LatLon</code> from two angles expressed in radians.
*
* @param latitude in radians
* @param longitude in radians
*
* @return a new <code>LatLon</code> from the given angles, which are expressed as radians
*/
public static LatLon fromRadians(double latitude, double longitude)
{
return new LatLon(Math.toDegrees(latitude), Math.toDegrees(longitude));
}
/**
* Factory method for obtaining a new <code>LatLon</code> from two angles expressed in degrees.
*
* @param latitude in degrees
* @param longitude in degrees
*
* @return a new <code>LatLon</code> from the given angles, which are expressed as degrees
*/
public static LatLon fromDegrees(double latitude, double longitude)
{
return new LatLon(latitude, longitude);
}
private LatLon(double latitude, double longitude)
{
this.latitude = Angle.fromDegrees(latitude);
this.longitude = Angle.fromDegrees(longitude);
}
/**
* Constructs a new <code>LatLon</code> from two angles. Neither angle may be null.
*
* @param latitude latitude
* @param longitude longitude
*
* @throws IllegalArgumentException if <code>latitude</code> or <code>longitude</code> is null
*/
public LatLon(Angle latitude, Angle longitude)
{
if (latitude == null || longitude == null)
{
throw new IllegalArgumentException("Latitude Or Longitude Is Null");
}
this.latitude = latitude;
this.longitude = longitude;
}
public LatLon(LatLon latLon)
{
if (latLon == null)
{
throw new IllegalArgumentException("Lat Lon Is Null");
}
this.latitude = latLon.latitude;
this.longitude = latLon.longitude;
}
/**
* Obtains the latitude of this <code>LatLon</code>.
*
* @return this <code>LatLon</code>'s latitude
*/
public final Angle getLatitude()
{
return this.latitude;
}
/**
* Obtains the longitude of this <code>LatLon</code>.
*
* @return this <code>LatLon</code>'s longitude
*/
public final Angle getLongitude()
{
return this.longitude;
}
/**
* Returns an array of this object's latitude and longitude in degrees.
*
* @return the array of latitude and longitude, arranged in that order.
*/
public double[] asDegreesArray()
{
return new double[] {this.getLatitude().degrees, this.getLongitude().degrees};
}
/**
* Returns an array of this object's latitude and longitude in radians.
*
* @return the array of latitude and longitude, arranged in that order.
*/
public double[] asRadiansArray()
{
return new double[] {this.getLatitude().radians, this.getLongitude().radians};
}
/**
* Returns the linear interpolation of <code>value1</code> and <code>value2</code>, treating the geographic
* locations as simple 2D coordinate pairs.
*
* @param amount the interpolation factor
* @param value1 the first location.
* @param value2 the second location.
*
* @return the linear interpolation of <code>value1</code> and <code>value2</code>.
*
* @throws IllegalArgumentException if either location is null.
*/
public static LatLon interpolate(double amount, LatLon value1, LatLon value2)
{
if (value1 == null || value2 == null)
{
throw new IllegalArgumentException("Lat Lon Is Null");
}
if (LatLon.equals(value1, value2))
return value1;
Line line;
try
{
line = Line.fromSegment(
new Vec4(value1.getLongitude().radians, value1.getLatitude().radians, 0),
new Vec4(value2.getLongitude().radians, value2.getLatitude().radians, 0));
}
catch (IllegalArgumentException e)
{
// Locations became coincident after calculations.
return value1;
}
Vec4 p = line.getPointAt(amount);
return LatLon.fromRadians(p.y(), p.x);
}
/**
* Returns the an interpolated location along the great-arc between <code>value1</code> and <code>value2</code>. The
* interpolation factor <code>amount</code> defines the weight given to each value, and is clamped to the range [0,
* 1]. If <code>a</code> is 0 or less, this returns <code>value1</code>. If <code>amount</code> is 1 or more, this
* returns <code>value2</code>. Otherwise, this returns the location on the great-arc between <code>value1</code>
* and <code>value2</code> corresponding to the specified interpolation factor.
*
* @param amount the interpolation factor
* @param value1 the first location.
* @param value2 the second location.
*
* @return an interpolated location along the great-arc between <code>value1</code> and <code>value2</code>.
*
* @throws IllegalArgumentException if either location is null.
*/
public static LatLon interpolateGreatCircle(double amount, LatLon value1, LatLon value2)
{
if (value1 == null || value2 == null)
{
throw new IllegalArgumentException("Lat Lon Is Null");
}
if (LatLon.equals(value1, value2))
return value1;
double t = WWMath.clamp(amount, 0d, 1d);
Angle azimuth = LatLon.greatCircleAzimuth(value1, value2);
Angle distance = LatLon.greatCircleDistance(value1, value2);
Angle pathLength = Angle.fromDegrees(t * distance.degrees);
return LatLon.greatCircleEndPosition(value1, azimuth, pathLength);
}
/**
* Returns the an interpolated location along the rhumb line between <code>value1</code> and <code>value2</code>.
* The interpolation factor <code>amount</code> defines the weight given to each value, and is clamped to the range
* [0, 1]. If <code>a</code> is 0 or less, this returns <code>value1</code>. If <code>amount</code> is 1 or more,
* this returns <code>value2</code>. Otherwise, this returns the location on the rhumb line between
* <code>value1</code> and <code>value2</code> corresponding to the specified interpolation factor.
*
* @param amount the interpolation factor
* @param value1 the first location.
* @param value2 the second location.
*
* @return an interpolated location along the rhumb line between <code>value1</code> and <code>value2</code>
*
* @throws IllegalArgumentException if either location is null.
*/
public static LatLon interpolateRhumb(double amount, LatLon value1, LatLon value2)
{
if (value1 == null || value2 == null)
{
throw new IllegalArgumentException("Lat Lon Is Null");
}
if (LatLon.equals(value1, value2))
return value1;
double t = WWMath.clamp(amount, 0d, 1d);
Angle azimuth = LatLon.rhumbAzimuth(value1, value2);
Angle distance = LatLon.rhumbDistance(value1, value2);
Angle pathLength = Angle.fromDegrees(t * distance.degrees);
return LatLon.rhumbEndPosition(value1, azimuth, pathLength);
}
/**
* Computes the great circle angular distance between two locations. The return value gives the distance as the
* angle between the two positions on the pi radius circle. In radians, this angle is also the arc length of the
* segment between the two positions on that circle. To compute a distance in meters from this value, multiply it by
* the radius of the globe.
*
* @param p1 LatLon of the first location
* @param p2 LatLon of the second location
*
* @return the angular distance between the two locations. In radians, this value is the arc length on the radius pi
* circle.
*/
public static Angle greatCircleDistance(LatLon p1, LatLon p2)
{
if ((p1 == null) || (p2 == null))
{
throw new IllegalArgumentException("Lat Lon Is Null");
}
double lat1 = p1.getLatitude().radians;
double lon1 = p1.getLongitude().radians;
double lat2 = p2.getLatitude().radians;
double lon2 = p2.getLongitude().radians;
if (lat1 == lat2 && lon1 == lon2)
return Angle.ZERO;
// "Haversine formula," taken from http://en.wikipedia.org/wiki/Great-circle_distance#Formul.C3.A6
double a = Math.sin((lat2 - lat1) / 2.0);
double b = Math.sin((lon2 - lon1) / 2.0);
double c = a * a + +Math.cos(lat1) * Math.cos(lat2) * b * b;
double distanceRadians = 2.0 * Math.asin(Math.sqrt(c));
return Double.isNaN(distanceRadians) ? Angle.ZERO : Angle.fromRadians(distanceRadians);
}
/**
* Computes the azimuth angle (clockwise from North) that points from the first location to the second location.
* This angle can be used as the starting azimuth for a great circle arc that begins at the first location, and
* passes through the second location.
*
* @param p1 LatLon of the first location
* @param p2 LatLon of the second location
*
* @return Angle that points from the first location to the second location.
*/
public static Angle greatCircleAzimuth(LatLon p1, LatLon p2)
{
if ((p1 == null) || (p2 == null))
{
throw new IllegalArgumentException("Lat Lon Is Null");
}
double lat1 = p1.getLatitude().radians;
double lon1 = p1.getLongitude().radians;
double lat2 = p2.getLatitude().radians;
double lon2 = p2.getLongitude().radians;
if (lat1 == lat2 && lon1 == lon2)
return Angle.ZERO;
if (lon1 == lon2)
return lat1 > lat2 ? Angle.POS180 : Angle.ZERO;
// Taken from "Map Projections - A Working Manual", page 30, equation 5-4b.
// The atan2() function is used in place of the traditional atan(y/x) to simplify the case when x==0.
double y = Math.cos(lat2) * Math.sin(lon2 - lon1);
double x = Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1) * Math.cos(lat2) * Math.cos(lon2 - lon1);
double azimuthRadians = Math.atan2(y, x);
return Double.isNaN(azimuthRadians) ? Angle.ZERO : Angle.fromRadians(azimuthRadians);
}
/**
* Computes the location on a great circle arc with the given starting location, azimuth, and arc distance.
*
* @param p LatLon of the starting location
* @param greatCircleAzimuth great circle azimuth angle (clockwise from North)
* @param pathLength arc distance to travel
*
* @return LatLon location on the great circle arc.
*/
public static LatLon greatCircleEndPosition(LatLon p, Angle greatCircleAzimuth, Angle pathLength)
{
if (p == null)
{
throw new IllegalArgumentException("Lat Lon Is Null");
}
if (greatCircleAzimuth == null || pathLength == null)
{
throw new IllegalArgumentException("Angle Is Null");
}
double lat = p.getLatitude().radians;
double lon = p.getLongitude().radians;
double azimuth = greatCircleAzimuth.radians;
double distance = pathLength.radians;
if (distance == 0)
return p;
// Taken from "Map Projections - A Working Manual", page 31, equation 5-5 and 5-6.
double endLatRadians = Math.asin(Math.sin(lat) * Math.cos(distance)
+ Math.cos(lat) * Math.sin(distance) * Math.cos(azimuth));
double endLonRadians = lon + Math.atan2(
Math.sin(distance) * Math.sin(azimuth),
Math.cos(lat) * Math.cos(distance) - Math.sin(lat) * Math.sin(distance) * Math.cos(azimuth));
if (Double.isNaN(endLatRadians) || Double.isNaN(endLonRadians))
return p;
return new LatLon(
Angle.fromRadians(endLatRadians).normalizedLatitude(),
Angle.fromRadians(endLonRadians).normalizedLongitude());
}
/**
* Computes the location on a great circle arc with the given starting location, azimuth, and arc distance.
*
* @param p LatLon of the starting location
* @param greatCircleAzimuthRadians great circle azimuth angle (clockwise from North), in radians
* @param pathLengthRadians arc distance to travel, in radians
*
* @return LatLon location on the great circle arc.
*/
public static LatLon greatCircleEndPosition(LatLon p, double greatCircleAzimuthRadians, double pathLengthRadians)
{
if (p == null)
{
throw new IllegalArgumentException("Lat Lon Is Null");
}
return greatCircleEndPosition(p,
Angle.fromRadians(greatCircleAzimuthRadians), Angle.fromRadians(pathLengthRadians));
}
/**
* Returns two locations with the most extreme latitudes on the great circle with the given starting location and
* azimuth.
*
* @param location location on the great circle.
* @param azimuth great circle azimuth angle (clockwise from North).
*
* @return two locations where the great circle has its extreme latitudes.
*
* @throws IllegalArgumentException if either <code>location</code> or <code>azimuth</code> are null.
*/
public static LatLon[] greatCircleExtremeLocations(LatLon location, Angle azimuth)
{
if (location == null)
{
throw new IllegalArgumentException("Location Is Null");
}
if (azimuth == null)
{
throw new IllegalArgumentException("Azimuth Is Null");
}
double lat0 = location.getLatitude().radians;
double az = azimuth.radians;
// Derived by solving the function for longitude on a great circle against the desired longitude. We start with
// the equation in "Map Projections - A Working Manual", page 31, equation 5-5:
//
// lat = asin( sin(lat0) * cos(c) + cos(lat0) * sin(c) * cos(Az) )
//
// Where (lat0, lon) are the starting coordinates, c is the angular distance along the great circle from the
// starting coordinate, and Az is the azimuth. All values are in radians.
//
// Solving for angular distance gives distance to the equator:
//
// tan(c) = -tan(lat0) / cos(Az)
//
// The great circle is by definition centered about the Globe's origin. Therefore intersections with the
// equator will be antipodal (exactly 180 degrees opposite each other), as will be the extreme latitudes.
// My observing the symmetry of a great circle, it is also apparent that the extreme latitudes will be 90
// degrees from either interseciton with the equator.
//
// d1 = c + 90
// d2 = c - 90
double tanDistance = -Math.tan(lat0) / Math.cos(az);
double distance = Math.atan(tanDistance);
Angle extremeDistance1 = Angle.fromRadians(distance + (Math.PI / 2.0));
Angle extremeDistance2 = Angle.fromRadians(distance - (Math.PI / 2.0));
return new LatLon[]
{
greatCircleEndPosition(location, azimuth, extremeDistance1),
greatCircleEndPosition(location, azimuth, extremeDistance2)
};
}
/**
* Returns two locations with the most extreme latitudes on the great circle arc defined by, and limited to, the two
* locations.
*
* @param begin beginning location on the great circle arc.
* @param end ending location on the great circle arc.
*
* @return two locations with the most extreme latitudes on the great circle arc.
*
* @throws IllegalArgumentException if either <code>begin</code> or <code>end</code> are null.
*/
public static LatLon[] greatCircleArcExtremeLocations(LatLon begin, LatLon end)
{
if (begin == null)
{
throw new IllegalArgumentException("Begin Is Null");
}
if (end == null)
{
throw new IllegalArgumentException("End Is Null");
}
LatLon minLatLocation = null;
LatLon maxLatLocation = null;
double minLat = Angle.POS90.degrees;
double maxLat = Angle.NEG90.degrees;
// Compute the min and max latitude and assocated locations from the arc endpoints.
for (LatLon ll : java.util.Arrays.asList(begin, end))
{
if (minLat >= ll.getLatitude().degrees)
{
minLat = ll.getLatitude().degrees;
minLatLocation = ll;
}
if (maxLat <= ll.getLatitude().degrees)
{
maxLat = ll.getLatitude().degrees;
maxLatLocation = ll;
}
}
// Compute parameters for the great circle arc defined by begin and end. Then compute the locations of extreme
// latitude on entire the great circle which that arc is part of.
Angle greatArcAzimuth = greatCircleAzimuth(begin, end);
Angle greatArcDistance = greatCircleDistance(begin, end);
LatLon[] greatCircleExtremes = greatCircleExtremeLocations(begin, greatArcAzimuth);
// Determine whether either of the extreme locations are inside the arc defined by begin and end. If so,
// adjust the min and max latitude accordingly.
for (LatLon ll : greatCircleExtremes)
{
Angle az = LatLon.greatCircleAzimuth(begin, ll);
Angle d = LatLon.greatCircleDistance(begin, ll);
// The extreme location must be between the begin and end locations. Therefore its azimuth relative to
// the begin location should have the same signum, and its distance relative to the begin location should
// be between 0 and greatArcDistance, inclusive.
if (Math.signum(az.degrees) == Math.signum(greatArcAzimuth.degrees))
{
if (d.degrees >= 0 && d.degrees <= greatArcDistance.degrees)
{
if (minLat >= ll.getLatitude().degrees)
{
minLat = ll.getLatitude().degrees;
minLatLocation = ll;
}
if (maxLat <= ll.getLatitude().degrees)
{
maxLat = ll.getLatitude().degrees;
maxLatLocation = ll;
}
}
}
}
return new LatLon[] {minLatLocation, maxLatLocation};
}
/**
* Returns two locations with the most extreme latitudes on the sequence of great circle arcs defined by each pair
* of locations in the specified iterable.
*
* @param locations the pairs of locations defining a sequence of great circle arcs.
*
* @return two locations with the most extreme latitudes on the great circle arcs.
*
* @throws IllegalArgumentException if <code>locations</code> is null.
*/
public static LatLon[] greatCircleArcExtremeLocations(Iterable<? extends LatLon> locations)
{
if (locations == null)
{
throw new IllegalArgumentException("Locations List Is Null");
}
LatLon minLatLocation = null;
LatLon maxLatLocation = null;
LatLon lastLocation = null;
for (LatLon ll : locations)
{
if (lastLocation != null)
{
LatLon[] extremes = LatLon.greatCircleArcExtremeLocations(lastLocation, ll);
if (extremes == null)
continue;
if (minLatLocation == null || minLatLocation.getLatitude().degrees > extremes[0].getLatitude().degrees)
minLatLocation = extremes[0];
if (maxLatLocation == null || maxLatLocation.getLatitude().degrees < extremes[1].getLatitude().degrees)
maxLatLocation = extremes[1];
}
lastLocation = ll;
}
return new LatLon[] {minLatLocation, maxLatLocation};
}
/**
* Computes the length of the rhumb line between two locations. The return value gives the distance as the angular
* distance between the two positions on the pi radius circle. In radians, this angle is also the arc length of the
* segment between the two positions on that circle. To compute a distance in meters from this value, multiply it by
* the radius of the globe.
*
* @param p1 LatLon of the first location
* @param p2 LatLon of the second location
*
* @return the arc length of the rhumb line between the two locations. In radians, this value is the arc length on
* the radius pi circle.
*/
public static Angle rhumbDistance(LatLon p1, LatLon p2)
{
if (p1 == null || p2 == null)
{
throw new IllegalArgumentException("LatLon Is Null");
}
double lat1 = p1.getLatitude().radians;
double lon1 = p1.getLongitude().radians;
double lat2 = p2.getLatitude().radians;
double lon2 = p2.getLongitude().radians;
if (lat1 == lat2 && lon1 == lon2)
return Angle.ZERO;
// Taken from http://www.movable-type.co.uk/scripts/latlong.html
double dLat = lat2 - lat1;
double dLon = lon2 - lon1;
double dPhi = Math.log(Math.tan(lat2 / 2.0 + Math.PI / 4.0) / Math.tan(lat1 / 2.0 + Math.PI / 4.0));
double q = dLat / dPhi;
if (Double.isNaN(dPhi) || Double.isNaN(q))
{
q = Math.cos(lat1);
}
// If lonChange over 180 take shorter rhumb across 180 meridian.
if (Math.abs(dLon) > Math.PI)
{
dLon = dLon > 0 ? -(2 * Math.PI - dLon) : (2 * Math.PI + dLon);
}
double distanceRadians = Math.sqrt(dLat * dLat + q * q * dLon * dLon);
return Double.isNaN(distanceRadians) ? Angle.ZERO : Angle.fromRadians(distanceRadians);
}
/**
* Computes the azimuth angle (clockwise from North) of a rhumb line (a line of constant heading) between two
* locations.
*
* @param p1 LatLon of the first location
* @param p2 LatLon of the second location
*
* @return azimuth Angle of a rhumb line between the two locations.
*/
public static Angle rhumbAzimuth(LatLon p1, LatLon p2)
{
if (p1 == null || p2 == null)
{
throw new IllegalArgumentException("LatLon Is Null");
}
double lat1 = p1.getLatitude().radians;
double lon1 = p1.getLongitude().radians;
double lat2 = p2.getLatitude().radians;
double lon2 = p2.getLongitude().radians;
if (lat1 == lat2 && lon1 == lon2)
return Angle.ZERO;
// Taken from http://www.movable-type.co.uk/scripts/latlong.html
double dLon = lon2 - lon1;
double dPhi = Math.log(Math.tan(lat2 / 2.0 + Math.PI / 4.0) / Math.tan(lat1 / 2.0 + Math.PI / 4.0));
// If lonChange over 180 take shorter rhumb across 180 meridian.
if (Math.abs(dLon) > Math.PI)
{
dLon = dLon > 0 ? -(2 * Math.PI - dLon) : (2 * Math.PI + dLon);
}
double azimuthRadians = Math.atan2(dLon, dPhi);
return Double.isNaN(azimuthRadians) ? Angle.ZERO : Angle.fromRadians(azimuthRadians);
}
/**
* Computes the location on a rhumb line with the given starting location, rhumb azimuth, and arc distance along the
* line.
*
* @param p LatLon of the starting location
* @param rhumbAzimuth rhumb azimuth angle (clockwise from North)
* @param pathLength arc distance to travel
*
* @return LatLon location on the rhumb line.
*/
public static LatLon rhumbEndPosition(LatLon p, Angle rhumbAzimuth, Angle pathLength)
{
if (p == null)
{
throw new IllegalArgumentException("LatLon Is Null");
}
if (rhumbAzimuth == null || pathLength == null)
{
throw new IllegalArgumentException("Angle Is Null");
}
double lat1 = p.getLatitude().radians;
double lon1 = p.getLongitude().radians;
double azimuth = rhumbAzimuth.radians;
double distance = pathLength.radians;
if (distance == 0)
return p;
// Taken from http://www.movable-type.co.uk/scripts/latlong.html
double lat2 = lat1 + distance * Math.cos(azimuth);
double dPhi = Math.log(Math.tan(lat2 / 2.0 + Math.PI / 4.0) / Math.tan(lat1 / 2.0 + Math.PI / 4.0));
double q = (lat2 - lat1) / dPhi;
if (Double.isNaN(dPhi) || Double.isNaN(q) || Double.isInfinite(q))
{
q = Math.cos(lat1);
}
double dLon = distance * Math.sin(azimuth) / q;
// Handle latitude passing over either pole.
if (Math.abs(lat2) > Math.PI / 2.0)
{
lat2 = lat2 > 0 ? Math.PI - lat2 : -Math.PI - lat2;
}
double lon2 = (lon1 + dLon + Math.PI) % (2 * Math.PI) - Math.PI;
if (Double.isNaN(lat2) || Double.isNaN(lon2))
return p;
return new LatLon(
Angle.fromRadians(lat2).normalizedLatitude(),
Angle.fromRadians(lon2).normalizedLongitude());
}
/**
* Computes the location on a rhumb line with the given starting location, rhumb azimuth, and arc distance along the
* line.
*
* @param p LatLon of the starting location
* @param rhumbAzimuthRadians rhumb azimuth angle (clockwise from North), in radians
* @param pathLengthRadians arc distance to travel, in radians
*
* @return LatLon location on the rhumb line.
*/
public static LatLon rhumbEndPosition(LatLon p, double rhumbAzimuthRadians, double pathLengthRadians)
{
if (p == null)
{
throw new IllegalArgumentException("LatLon Is Null");
}
return rhumbEndPosition(p, Angle.fromRadians(rhumbAzimuthRadians), Angle.fromRadians(pathLengthRadians));
}
public static Angle getAverageDistance(Iterable<? extends LatLon> locations)
{
// Compute the average rhumb distance between locations.
if ((locations == null))
{
throw new IllegalArgumentException("Locations List Is Null");
}
double totalDistance = 0.0;
int count = 0;
for (LatLon p1 : locations)
{
for (LatLon p2 : locations)
{
if (p1 != p2)
{
double d = rhumbDistance(p1, p2).radians;
totalDistance += d;
count++;
}
}
}
return (count == 0) ? Angle.ZERO : Angle.fromRadians(totalDistance / (double) count);
}
public LatLon add(LatLon that)
{
if (that == null)
{
throw new IllegalArgumentException("Angle Is Null");
}
Angle lat = Angle.normalizedLatitude(this.latitude.add(that.latitude));
Angle lon = Angle.normalizedLongitude(this.longitude.add(that.longitude));
return new LatLon(lat, lon);
}
public LatLon subtract(LatLon that)
{
if (that == null)
{
throw new IllegalArgumentException("Angle Is Null");
}
Angle lat = Angle.normalizedLatitude(this.latitude.subtract(that.latitude));
Angle lon = Angle.normalizedLongitude(this.longitude.subtract(that.longitude));
return new LatLon(lat, lon);
}
public LatLon add(Position that)
{
if (that == null)
{
throw new IllegalArgumentException("Angle Is Null");
}
Angle lat = Angle.normalizedLatitude(this.latitude.add(that.getLatitude()));
Angle lon = Angle.normalizedLongitude(this.longitude.add(that.getLongitude()));
return new LatLon(lat, lon);
}
public LatLon subtract(Position that)
{
if (that == null)
{
throw new IllegalArgumentException("Angle Is Null");
}
Angle lat = Angle.normalizedLatitude(this.latitude.subtract(that.getLatitude()));
Angle lon = Angle.normalizedLongitude(this.longitude.subtract(that.getLongitude()));
return new LatLon(lat, lon);
}
public static boolean locationsCrossDateLine(Iterable<? extends LatLon> locations)
{
if (locations == null)
{
throw new IllegalArgumentException("Locations List Is Null");
}
LatLon pos = null;
for (LatLon posNext : locations)
{
if (pos != null)
{
// A segment cross the line if end pos have different longitude signs
// and are more than 180 degrees longitude apart
if (Math.signum(pos.getLongitude().degrees) != Math.signum(posNext.getLongitude().degrees))
{
double delta = Math.abs(pos.getLongitude().degrees - posNext.getLongitude().degrees);
if (delta > 180 && delta < 360)
return true;
}
}
pos = posNext;
}
return false;
}
public static boolean locationsCrossDateline(LatLon p1, LatLon p2)
{
if (p1 == null || p2 == null)
{
throw new IllegalArgumentException("Location Is Null");
}
// A segment cross the line if end pos have different longitude signs
// and are more than 180 degrees longitude apart
if (Math.signum(p1.getLongitude().degrees) != Math.signum(p2.getLongitude().degrees))
{
double delta = Math.abs(p1.getLongitude().degrees - p2.getLongitude().degrees);
if (delta > 180 && delta < 360)
return true;
}
return false;
}
/**
* Transform the negative longitudes of a dateline-spanning location list to positive values that maintain the
* relationship with the other locations in the list. Negative longitudes are transformed to values greater than 180
* degrees, as though longitude spanned [0, 360] rather than [-180, 180]. This enables arithmetic operations to be
* performed on the locations without having to take into account the longitude jump at the dateline.
*
* @param locations the locations to transform. This list is not modified.
*
* @return a new list of locations transformed as described above.
*
* @throws IllegalArgumentException if the location list is null.
*/
public static List<LatLon> makeDatelineCrossingLocationsPositive(Iterable<? extends LatLon> locations)
{
if (locations == null)
{
throw new IllegalArgumentException("Locations List Is Null");
}
Iterator<? extends LatLon> iter = locations.iterator();
if (!iter.hasNext())
return Collections.emptyList();
ArrayList<LatLon> newLocations = new ArrayList<LatLon>();
for (LatLon location : locations)
{
if (location == null)
continue;
if (location.getLongitude().degrees < 0)
{
newLocations.add(
LatLon.fromDegrees(location.getLatitude().degrees, location.getLongitude().degrees + 360));
}
else
{
newLocations.add(location);
}
}
return newLocations;
}
/**
* Parses a string containing latitude and longitude coordinates in either Degrees-minutes-seconds or decimal
* degrees. The latitude must precede the longitude and the angles must be separated by a comma.
*
* @param latLonString a string containing the comma separated latitude and longitude in either DMS or decimal
* degrees.
*
* @return a <code>LatLon</code> instance with the parsed angles.
*
* @throws IllegalArgumentException if <code>latLonString</code> is null.
* @throws NumberFormatException if the string does not form a latitude, longitude pair.
*/
public LatLon parseLatLon(String latLonString) // TODO
{
if (latLonString == null)
{
throw new IllegalArgumentException("String Is Null");
}
throw new UnsupportedOperationException(); // TODO: remove when implemented
}
@Override
public String toString()
{
String las = String.format("Lat %7.4f\u00B0", this.getLatitude().getDegrees());
String los = String.format("Lon %7.4f\u00B0", this.getLongitude().getDegrees());
return "(" + las + ", " + los + ")";
}
@Override
public boolean equals(Object o)
{
if (this == o)
return true;
if (o == null || getClass() != o.getClass())
return false;
final gov.nasa.worldwind.geom.LatLon latLon = (gov.nasa.worldwind.geom.LatLon) o;
if (!latitude.equals(latLon.latitude))
return false;
//noinspection RedundantIfStatement
if (!longitude.equals(latLon.longitude))
return false;
return true;
}
public static boolean equals(LatLon a, LatLon b)
{
return a.getLatitude().equals(b.getLatitude()) && a.getLongitude().equals(b.getLongitude());
}
@Override
public int hashCode()
{
int result;
result = latitude.hashCode();
result = 29 * result + longitude.hashCode();
return result;
}
/**
* Compute the forward azimuth between two positions
*
* @param p1 first position
* @param p2 second position
* @param equatorialRadius the equatorial radius of the globe in meters
* @param polarRadius the polar radius of the globe in meters
*
* @return the azimuth
*/
public static Angle ellipsoidalForwardAzimuth(LatLon p1, LatLon p2, double equatorialRadius, double polarRadius)
{
if (p1 == null || p2 == null)
{
throw new IllegalArgumentException("Position Is Null");
}
// TODO: What if polar radius is larger than equatorial radius?
// Calculate flattening
final double f = (equatorialRadius - polarRadius) / equatorialRadius; // flattening
// Calculate reduced latitudes and related sines/cosines
final double U1 = Math.atan((1.0 - f) * Math.tan(p1.latitude.radians));
final double cU1 = Math.cos(U1);
final double sU1 = Math.sin(U1);
final double U2 = Math.atan((1.0 - f) * Math.tan(p2.latitude.radians));
final double cU2 = Math.cos(U2);
final double sU2 = Math.sin(U2);
// Calculate difference in longitude
final double L = p2.longitude.subtract(p1.longitude).radians;
// Vincenty's Formula for Forward Azimuth
// iterate until change in lambda is negligible (e.g. 1e-12 ~= 0.06mm)
// first approximation
double lambda = L;
double sLambda = Math.sin(lambda);
double cLambda = Math.cos(lambda);
// dummy value to ensure
double lambda_prev = Double.MAX_VALUE;
int count = 0;
while (Math.abs(lambda - lambda_prev) > 1e-12 && count++ < 100)
{
// Store old lambda
lambda_prev = lambda;
// Calculate new lambda
double sSigma = Math.sqrt(Math.pow(cU2 * sLambda, 2)
+ Math.pow(cU1 * sU2 - sU1 * cU2 * cLambda, 2));
double cSigma = sU1 * sU2 + cU1 * cU2 * cLambda;
double sigma = Math.atan2(sSigma, cSigma);
double sAlpha = cU1 * cU2 * sLambda / sSigma;
double cAlpha2 = 1 - sAlpha * sAlpha; // trig identity
// As cAlpha2 approaches zeros, set cSigmam2 to zero to converge on a solution
double cSigmam2;
if (Math.abs(cAlpha2) < 1e-6)
{
cSigmam2 = 0;
}
else
{
cSigmam2 = cSigma - 2 * sU1 * sU2 / cAlpha2;
}
double c = f / 16 * cAlpha2 * (4 + f * (4 - 3 * cAlpha2));
lambda = L + (1 - c) * f * sAlpha * (sigma + c * sSigma * (cSigmam2 + c * cSigma * (-1 + 2 * cSigmam2)));
sLambda = Math.sin(lambda);
cLambda = Math.cos(lambda);
}
return Angle.fromRadians(Math.atan2(cU2 * sLambda, cU1 * sU2 - sU1 * cU2 * cLambda));
}
// TODO: Need method to compute end position from initial position, azimuth and distance. The companion to the
// spherical version, endPosition(), above.
/**
* Computes the distance between two points on an ellipsoid iteratively.
* <p/>
* NOTE: This method was copied from the UniData NetCDF Java library. http://www.unidata.ucar.edu/software/netcdf-java/
* <p/>
* Algorithm from U.S. National Geodetic Survey, FORTRAN program "inverse," subroutine "INVER1," by L. PFEIFER and
* JOHN G. GERGEN. See http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html
* <p/>
* Original documentation: SOLUTION OF THE GEODETIC INVERSE PROBLEM AFTER T.VINCENTY MODIFIED RAINSFORD'S METHOD
* WITH HELMERT'S ELLIPTICAL TERMS EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
* STANDPOINT/FOREPOINT MUST NOT BE THE GEOGRAPHIC POLE
* <p/>
* Requires close to 1.4 E-5 seconds wall clock time per call on a 550 MHz Pentium with Linux 7.2.
*
* @param p1 first position
* @param p2 second position
* @param equatorialRadius the equatorial radius of the globe in meters
* @param polarRadius the polar radius of the globe in meters
*
* @return distance in meters between the two points
*/
public static double ellipsoidalDistance(LatLon p1, LatLon p2, double equatorialRadius, double polarRadius)
{
// TODO: I think there is a non-iterative way to calculate the distance. Find it and compare with this one.
// TODO: What if polar radius is larger than equatorial radius?
final double F = (equatorialRadius - polarRadius) / equatorialRadius; // flattening = 1.0 / 298.257223563;
final double R = 1.0 - F;
final double EPS = 0.5E-13;
if (p1 == null || p2 == null)
{
throw new IllegalArgumentException("Position Is Null");
}
// Algorithm from National Geodetic Survey, FORTRAN program "inverse,"
// subroutine "INVER1," by L. PFEIFER and JOHN G. GERGEN.
// http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html
// Conversion to JAVA from FORTRAN was made with as few changes as possible
// to avoid errors made while recasting form, and to facilitate any future
// comparisons between the original code and the altered version in Java.
// Original documentation:
// SOLUTION OF THE GEODETIC INVERSE PROBLEM AFTER T.VINCENTY
// MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
// EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
// STANDPOINT/FOREPOINT MUST NOT BE THE GEOGRAPHIC POLE
// A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID
// F IS THE FLATTENING (NOT RECIPROCAL) OF THE REFERNECE ELLIPSOID
// LATITUDES GLAT1 AND GLAT2
// AND LONGITUDES GLON1 AND GLON2 ARE IN RADIANS POSITIVE NORTH AND EAST
// FORWARD AZIMUTHS AT BOTH POINTS RETURNED IN RADIANS FROM NORTH
//
// Reference ellipsoid is the WGS-84 ellipsoid.
// See http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
// FAZ is forward azimuth in radians from pt1 to pt2;
// BAZ is backward azimuth from point 2 to 1;
// S is distance in meters.
//
// Conversion to JAVA from FORTRAN was made with as few changes as possible
// to avoid errors made while recasting form, and to facilitate any future
// comparisons between the original code and the altered version in Java.
//
//IMPLICIT REAL*8 (A-H,O-Z)
// COMMON/CONST/PI,RAD
double GLAT1 = p1.getLatitude().radians;
double GLAT2 = p2.getLatitude().radians;
double TU1 = R * Math.sin(GLAT1) / Math.cos(GLAT1);
double TU2 = R * Math.sin(GLAT2) / Math.cos(GLAT2);
double CU1 = 1. / Math.sqrt(TU1 * TU1 + 1.);
double SU1 = CU1 * TU1;
double CU2 = 1. / Math.sqrt(TU2 * TU2 + 1.);
double S = CU1 * CU2;
double BAZ = S * TU2;
double FAZ = BAZ * TU1;
double GLON1 = p1.getLongitude().radians;
double GLON2 = p2.getLongitude().radians;
double X = GLON2 - GLON1;
double D, SX, CX, SY, CY, Y, SA, C2A, CZ, E, C;
do
{
SX = Math.sin(X);
CX = Math.cos(X);
TU1 = CU2 * SX;
TU2 = BAZ - SU1 * CU2 * CX;
SY = Math.sqrt(TU1 * TU1 + TU2 * TU2);
CY = S * CX + FAZ;
Y = Math.atan2(SY, CY);
SA = S * SX / SY;
C2A = -SA * SA + 1.;
CZ = FAZ + FAZ;
if (C2A > 0.)
{
CZ = -CZ / C2A + CY;
}
E = CZ * CZ * 2. - 1.;
C = ((-3. * C2A + 4.) * F + 4.) * C2A * F / 16.;
D = X;
X = ((E * CY * C + CZ) * SY * C + Y) * SA;
X = (1. - C) * X * F + GLON2 - GLON1;
//IF(DABS(D-X).GT.EPS) GO TO 100
}
while (Math.abs(D - X) > EPS);
//FAZ = Math.atan2(TU1, TU2);
//BAZ = Math.atan2(CU1 * SX, BAZ * CX - SU1 * CU2) + Math.PI;
X = Math.sqrt((1. / R / R - 1.) * C2A + 1.) + 1.;
X = (X - 2.) / X;
C = 1. - X;
C = (X * X / 4. + 1.) / C;
D = (0.375 * X * X - 1.) * X;
X = E * CY;
S = 1. - E - E;
S = ((((SY * SY * 4. - 3.) * S * CZ * D / 6. - X) * D / 4. + CZ) * SY
* D + Y) * C * equatorialRadius * R;
return S;
}
/**
* Computes a new set of locations translated from a specified location to a new location.
*
* @param oldLocation the original reference location.
* @param newLocation the new reference location.
* @param locations the locations to translate.
*
* @return the translated locations, or null if the locations could not be translated.
*
* @throws IllegalArgumentException if any argument is null.
*/
public static List<LatLon> computeShiftedLocations(Position oldLocation, Position newLocation,
Iterable<? extends LatLon> locations)
{
// TODO: Account for dateline spanning
if (oldLocation == null || newLocation == null)
{
throw new IllegalArgumentException("Position Is Null");
}
if (locations == null)
{
throw new IllegalArgumentException("Positions List Is Null");
}
ArrayList<LatLon> newPositions = new ArrayList<LatLon>();
for (LatLon location : locations)
{
Angle distance = LatLon.greatCircleDistance(oldLocation, location);
Angle azimuth = LatLon.greatCircleAzimuth(oldLocation, location);
newPositions.add(Position.greatCircleEndPosition(newLocation, azimuth, distance));
}
return newPositions;
}
}