DistanceUtils.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You 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.
*/
// NOTE: we keep the header as it came from ASF; it did not originate in Spatial4j
package org.locationtech.spatial4j.distance;
import org.locationtech.spatial4j.context.SpatialContext;
import org.locationtech.spatial4j.shape.Point;
import org.locationtech.spatial4j.shape.Rectangle;
/**
* Various distance calculations and constants. To the extent possible, a {@link
* org.locationtech.spatial4j.distance.DistanceCalculator}, retrieved from {@link
* org.locationtech.spatial4j.context.SpatialContext#getDistCalc()} should be used in preference to calling
* these methods directly.
* <p>
* This code came from <a href="https://issues.apache.org/jira/browse/LUCENE-1387">Apache
* Lucene, LUCENE-1387</a>, which in turn came from "LocalLucene".
*/
public class DistanceUtils {
//pre-compute some angles that are commonly used
@Deprecated
public static final double DEG_45_AS_RADS = Math.PI / 4;
@Deprecated
public static final double SIN_45_AS_RADS = Math.sin(DEG_45_AS_RADS);
public static final double DEG_90_AS_RADS = Math.PI / 2;
public static final double DEG_180_AS_RADS = Math.PI;
@Deprecated
public static final double DEG_225_AS_RADS = 5 * DEG_45_AS_RADS;
@Deprecated
public static final double DEG_270_AS_RADS = 3 * DEG_90_AS_RADS;
public static final double DEGREES_TO_RADIANS = Math.PI / 180;
public static final double RADIANS_TO_DEGREES = 1 / DEGREES_TO_RADIANS;
public static final double KM_TO_MILES = 0.621371192;
public static final double MILES_TO_KM = 1 / KM_TO_MILES;//1.609
/**
* The International Union of Geodesy and Geophysics says the Earth's mean radius in KM is:
*
* [1] http://en.wikipedia.org/wiki/Earth_radius
*/
public static final double EARTH_MEAN_RADIUS_KM = 6371.0087714;
public static final double EARTH_EQUATORIAL_RADIUS_KM = 6378.1370;
/** Equivalent to degrees2Dist(1, EARTH_MEAN_RADIUS_KM) */
public static final double DEG_TO_KM = DEGREES_TO_RADIANS * EARTH_MEAN_RADIUS_KM;
public static final double KM_TO_DEG = 1 / DEG_TO_KM;
public static final double EARTH_MEAN_RADIUS_MI = EARTH_MEAN_RADIUS_KM * KM_TO_MILES;
public static final double EARTH_EQUATORIAL_RADIUS_MI = EARTH_EQUATORIAL_RADIUS_KM * KM_TO_MILES;
private DistanceUtils() {}
/**
* Calculate the p-norm (i.e. length) between two vectors.
* <p>
* See <a href="http://en.wikipedia.org/wiki/Lp_space">Lp space</a>
*
* @param vec1 The first vector
* @param vec2 The second vector
* @param power The power (2 for cartesian distance, 1 for manhattan, etc.)
* @return The length.
*
* @see #vectorDistance(double[], double[], double, double)
*
*/
@Deprecated
public static double vectorDistance(double[] vec1, double[] vec2, double power) {
//only calc oneOverPower if it's needed
double oneOverPower = (power == 0 || power == 1.0 || power == 2.0) ? Double.NaN : 1.0 / power;
return vectorDistance(vec1, vec2, power, oneOverPower);
}
/**
* Calculate the p-norm (i.e. length) between two vectors.
*
* @param vec1 The first vector
* @param vec2 The second vector
* @param power The power (2 for cartesian distance, 1 for manhattan, etc.)
* @param oneOverPower If you've pre-calculated oneOverPower and cached it, use this method to save
* one division operation over {@link #vectorDistance(double[], double[], double)}.
* @return The length.
*/
@Deprecated
public static double vectorDistance(double[] vec1, double[] vec2, double power, double oneOverPower) {
double result = 0;
if (power == 0) {
for (int i = 0; i < vec1.length; i++) {
result += vec1[i] - vec2[i] == 0 ? 0 : 1;
}
} else if (power == 1.0) { // Manhattan
for (int i = 0; i < vec1.length; i++) {
result += Math.abs(vec1[i] - vec2[i]);
}
} else if (power == 2.0) { // Cartesian
result = Math.sqrt(distSquaredCartesian(vec1, vec2));
} else if (power == Integer.MAX_VALUE || Double.isInfinite(power)) {//infinite norm?
for (int i = 0; i < vec1.length; i++) {
result = Math.max(result, Math.max(vec1[i], vec2[i]));
}
} else {
for (int i = 0; i < vec1.length; i++) {
result += Math.pow(vec1[i] - vec2[i], power);
}
result = Math.pow(result, oneOverPower);
}
return result;
}
/**
* Return the coordinates of a vector that is the corner of a box (upper right or lower left), assuming a Rectangular
* coordinate system. Note, this does not apply for points on a sphere or ellipse (although it could be used as an approximation).
*
* @param center The center point
* @param result Holds the result, potentially resizing if needed.
* @param distance The d from the center to the corner
* @param upperRight If true, return the coords for the upper right corner, else return the lower left.
* @return The point, either the upperLeft or the lower right
*/
@Deprecated
public static double[] vectorBoxCorner(double[] center, double[] result, double distance, boolean upperRight) {
if (result == null || result.length != center.length) {
result = new double[center.length];
}
if (!upperRight) {
distance = -distance;
}
//We don't care about the power here,
// b/c we are always in a rectangular coordinate system, so any norm can be used by
//using the definition of sine
distance = SIN_45_AS_RADS * distance; // sin(Pi/4) == (2^0.5)/2 == opp/hyp == opp/distance, solve for opp, similarly for cosine
for (int i = 0; i < center.length; i++) {
result[i] = center[i] + distance;
}
return result;
}
/**
* Given a start point (startLat, startLon), distance, and a bearing on a sphere, return the destination point.
*
* @param startLat The starting point latitude, in radians
* @param startLon The starting point longitude, in radians
* @param distanceRAD The distance to travel along the bearing in radians.
* @param bearingRAD The bearing, in radians. North is a 0, moving clockwise till radians(360).
* @param reuse A preallocated object to hold the results.
* @return The destination point, IN RADIANS.
*/
public static Point pointOnBearingRAD(double startLat, double startLon, double distanceRAD, double bearingRAD, SpatialContext ctx, Point reuse) {
/*
lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(��))
lon2 = lon1 + atan2(sin(��)*sin(d/R)*cos(lat1), cos(d/R)���sin(lat1)*sin(lat2))
*/
double cosAngDist = Math.cos(distanceRAD);
double cosStartLat = Math.cos(startLat);
double sinAngDist = Math.sin(distanceRAD);
double sinStartLat = Math.sin(startLat);
double sinLat2 = sinStartLat * cosAngDist +
cosStartLat * sinAngDist * Math.cos(bearingRAD);
double lat2 = Math.asin(sinLat2);
double lon2 = startLon + Math.atan2(Math.sin(bearingRAD) * sinAngDist * cosStartLat,
cosAngDist - sinStartLat * sinLat2);
// normalize lon first
if (lon2 > DEG_180_AS_RADS) {
lon2 = -1.0 * (DEG_180_AS_RADS - (lon2 - DEG_180_AS_RADS));
} else if (lon2 < -DEG_180_AS_RADS) {
lon2 = (lon2 + DEG_180_AS_RADS) + DEG_180_AS_RADS;
}
// normalize lat - could flip poles
if (lat2 > DEG_90_AS_RADS) {
lat2 = DEG_90_AS_RADS - (lat2 - DEG_90_AS_RADS);
if (lon2 < 0) {
lon2 = lon2 + DEG_180_AS_RADS;
} else {
lon2 = lon2 - DEG_180_AS_RADS;
}
} else if (lat2 < -DEG_90_AS_RADS) {
lat2 = -DEG_90_AS_RADS - (lat2 + DEG_90_AS_RADS);
if (lon2 < 0) {
lon2 = lon2 + DEG_180_AS_RADS;
} else {
lon2 = lon2 - DEG_180_AS_RADS;
}
}
if (reuse == null) {
return ctx.makePoint(lon2, lat2);
} else {
reuse.reset(lon2, lat2);//x y
return reuse;
}
}
/**
* Puts in range -180 <= lon_deg <= +180.
*/
public static double normLonDEG(double lon_deg) {
if (lon_deg >= -180 && lon_deg <= 180)
return lon_deg;//common case, and avoids slight double precision shifting
double off = (lon_deg + 180) % 360;
if (off < 0)
return 180 + off;
else if (off == 0 && lon_deg > 0)
return 180;
else
return -180 + off;
}
/**
* Puts in range -90 <= lat_deg <= 90.
*/
public static double normLatDEG(double lat_deg) {
if (lat_deg >= -90 && lat_deg <= 90)
return lat_deg;//common case, and avoids slight double precision shifting
double off = Math.abs((lat_deg + 90) % 360);
return (off <= 180 ? off : 360-off) - 90;
}
/**
* Calculates the bounding box of a circle, as specified by its center point
* and distance. <code>reuse</code> is an optional argument to store the
* results to avoid object creation.
*/
public static Rectangle calcBoxByDistFromPtDEG(double lat, double lon, double distDEG, SpatialContext ctx, Rectangle reuse) {
//See http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates Section 3.1, 3.2 and 3.3
double minX; double maxX; double minY; double maxY;
if (distDEG == 0) {
minX = lon; maxX = lon; minY = lat; maxY = lat;
} else if (distDEG >= 180) {//distance is >= opposite side of the globe
minX = -180; maxX = 180; minY = -90; maxY = 90;
} else {
//--calc latitude bounds
maxY = lat + distDEG;
minY = lat - distDEG;
if (maxY >= 90 || minY <= -90) {//touches either pole
//we have special logic for longitude
minX = -180; maxX = 180;//world wrap: 360 deg
if (maxY <= 90 && minY >= -90) {//doesn't pass either pole: 180 deg
minX = normLonDEG(lon - 90);
maxX = normLonDEG(lon + 90);
}
if (maxY > 90)
maxY = 90;
if (minY < -90)
minY = -90;
} else {
//--calc longitude bounds
double lon_delta_deg = calcBoxByDistFromPt_deltaLonDEG(lat, distDEG);
minX = normLonDEG(lon - lon_delta_deg);
maxX = normLonDEG(lon + lon_delta_deg);
}
}
if (reuse == null) {
return ctx.makeRectangle(minX, maxX, minY, maxY);
} else {
reuse.reset(minX, maxX, minY, maxY);
return reuse;
}
}
/**
* The delta longitude of a point-distance. In other words, half the width of
* the bounding box of a circle.
*/
public static double calcBoxByDistFromPt_deltaLonDEG(double lat, double distDEG) {
//http://gis.stackexchange.com/questions/19221/find-tangent-point-on-circle-furthest-east-or-west
if (distDEG == 0)
return 0;
double lat_rad = toRadians(lat);
double dist_rad = toRadians(distDEG);
double result_rad = Math.asin(Math.sin(dist_rad) / Math.cos(lat_rad));
if (!Double.isNaN(result_rad))
return toDegrees(result_rad);
return 90;
}
/**
* The latitude of the horizontal axis (e.g. left-right line)
* of a circle. The horizontal axis of a circle passes through its furthest
* left-most and right-most edges. On a 2D plane, this result is always
* <code>from.getY()</code> but, perhaps surprisingly, on a sphere it is going
* to be slightly different.
*/
public static double calcBoxByDistFromPt_latHorizAxisDEG(double lat, double lon, double distDEG) {
//http://gis.stackexchange.com/questions/19221/find-tangent-point-on-circle-furthest-east-or-west
if (distDEG == 0)
return lat;
// if we don't do this when == 90 or -90, computed result can be (+/-)89.9999 when at pole.
// No biggie but more accurate.
else if (lat + distDEG >= 90)
return 90;
else if (lat - distDEG <= -90)
return -90;
double lat_rad = toRadians(lat);
double dist_rad = toRadians(distDEG);
double result_rad = Math.asin( Math.sin(lat_rad) / Math.cos(dist_rad));
if (!Double.isNaN(result_rad))
return toDegrees(result_rad);
//handle NaN (shouldn't happen due to checks earlier)
if (lat > 0)
return 90;
if (lat < 0)
return -90;
return lat;//0
}
/**
* Calculates the degrees longitude distance at latitude {@code lat} to cover
* a distance {@code dist}.
* <p>
* Used to calculate a new expanded buffer distance to account for skewing
* effects for shapes that use the lat-lon space as a 2D plane instead of a
* sphere. The expanded buffer will be sure to cover the intended area, but
* the shape is still skewed and so it will cover a larger area. For latitude
* 0 (the equator) the result is the same buffer. At 60 (or -60) degrees, the
* result is twice the buffer, meaning that a shape at 60 degrees is twice as
* high as it is wide when projected onto a lat-lon plane even if in the real
* world it's equal all around.
* <p>
* If the result added to abs({@code lat}) is >= 90 degrees, then skewing is
* so severe that the caller should consider tossing the shape and
* substituting a spherical cap instead.
*
* @param lat latitude in degrees
* @param dist distance in degrees
* @return longitudinal degrees (x delta) at input latitude that is >= dist
* distance. Will be >= dist and <= 90.
*/
public static double calcLonDegreesAtLat(double lat, double dist) {
//This code was pulled out of DistanceUtils.pointOnBearingRAD() and
// optimized
// for bearing = 90 degrees, and so we can get an intermediate calculation.
double distanceRAD = DistanceUtils.toRadians(dist);
double startLat = DistanceUtils.toRadians(lat);
double cosAngDist = Math.cos(distanceRAD);
double cosStartLat = Math.cos(startLat);
double sinAngDist = Math.sin(distanceRAD);
double sinStartLat = Math.sin(startLat);
double lonDelta = Math.atan2(sinAngDist * cosStartLat,
cosAngDist * (1 - sinStartLat * sinStartLat));
return DistanceUtils.toDegrees(lonDelta);
}
/**
* The square of the cartesian Distance. Not really a distance, but useful if all that matters is
* comparing the result to another one.
*
* @param vec1 The first point
* @param vec2 The second point
* @return The squared cartesian distance
*/
@Deprecated
public static double distSquaredCartesian(double[] vec1, double[] vec2) {
double result = 0;
for (int i = 0; i < vec1.length; i++) {
double v = vec1[i] - vec2[i];
result += v * v;
}
return result;
}
/**
*
* @param lat1 The y coordinate of the first point, in radians
* @param lon1 The x coordinate of the first point, in radians
* @param lat2 The y coordinate of the second point, in radians
* @param lon2 The x coordinate of the second point, in radians
* @return The distance between the two points, as determined by the Haversine formula, in radians.
*/
public static double distHaversineRAD(double lat1, double lon1, double lat2, double lon2) {
//TODO investigate slightly different formula using asin() and min() http://www.movable-type.co.uk/scripts/gis-faq-5.1.html
// Check for same position
if (lat1 == lat2 && lon1 == lon2)
return 0.0;
double hsinX = Math.sin((lon1 - lon2) * 0.5);
double hsinY = Math.sin((lat1 - lat2) * 0.5);
double h = hsinY * hsinY +
(Math.cos(lat1) * Math.cos(lat2) * hsinX * hsinX);
if (h > 1)//numeric robustness issue. If we didn't check, the answer would be NaN!
h = 1;
return 2 * Math.atan2(Math.sqrt(h), Math.sqrt(1 - h));
}
/**
* Calculates the distance between two lat-lon's using the Law of Cosines. Due to numeric conditioning
* errors, it is not as accurate as the Haversine formula for small distances. But with
* double precision, it isn't that bad -- <a href="http://www.movable-type.co.uk/scripts/latlong.html">
* allegedly 1 meter</a>.
* <p>
* See <a href="http://gis.stackexchange.com/questions/4906/why-is-law-of-cosines-more-preferable-than-haversine-when-calculating-distance-b">
* Why is law of cosines more preferable than haversine when calculating distance between two latitude-longitude points?</a>
* <p>
* The arguments and return value are in radians.
*/
public static double distLawOfCosinesRAD(double lat1, double lon1, double lat2, double lon2) {
// Check for same position
if (lat1 == lat2 && lon1 == lon2)
return 0.0;
// Get the longitude difference. Don't need to worry about
// crossing dateline since cos(x) = cos(-x)
double dLon = lon2 - lon1;
double cosB = (Math.sin(lat1) * Math.sin(lat2))
+ (Math.cos(lat1) * Math.cos(lat2) * Math.cos(dLon));
// Find angle subtended (with some bounds checking) in radians
if (cosB < -1.0)
return Math.PI;
else if (cosB >= 1.0)
return 0;
else
return Math.acos(cosB);
}
/**
* Calculates the great circle distance using the Vincenty Formula, simplified for a spherical model. This formula
* is accurate for any pair of points. The equation
* was taken from <a href="http://en.wikipedia.org/wiki/Great-circle_distance">Wikipedia</a>.
* <p>
* The arguments are in radians, and the result is in radians.
*/
public static double distVincentyRAD(double lat1, double lon1, double lat2, double lon2) {
// Check for same position
if (lat1 == lat2 && lon1 == lon2)
return 0.0;
double cosLat1 = Math.cos(lat1);
double cosLat2 = Math.cos(lat2);
double sinLat1 = Math.sin(lat1);
double sinLat2 = Math.sin(lat2);
double dLon = lon2 - lon1;
double cosDLon = Math.cos(dLon);
double sinDLon = Math.sin(dLon);
double a = cosLat2 * sinDLon;
double b = cosLat1*sinLat2 - sinLat1*cosLat2*cosDLon;
double c = sinLat1*sinLat2 + cosLat1*cosLat2*cosDLon;
return Math.atan2(Math.sqrt(a*a+b*b),c);
}
/**
* Converts a distance in the units of the radius to degrees (360 degrees are
* in a circle). A spherical earth model is assumed.
*/
public static double dist2Degrees(double dist, double radius) {
return toDegrees(dist2Radians(dist, radius));
}
/**
* Converts <code>degrees</code> (1/360th of circumference of a circle) into a
* distance as measured by the units of the radius. A spherical earth model
* is assumed.
*/
public static double degrees2Dist(double degrees, double radius) {
return radians2Dist(toRadians(degrees), radius);
}
/**
* Converts a distance in the units of <code>radius</code> (e.g. kilometers)
* to radians (multiples of the radius). A spherical earth model is assumed.
*/
public static double dist2Radians(double dist, double radius) {
return dist / radius;
}
/**
* Converts <code>radians</code> (multiples of the <code>radius</code>) to
* distance in the units of the radius (e.g. kilometers).
*/
public static double radians2Dist(double radians, double radius) {
return radians * radius;
}
/**
* Same as {@link Math#toRadians(double)} but 3x faster (multiply vs. divide).
* See CompareRadiansSnippet.java in tests.
*/
public static double toRadians(double degrees) {
return degrees * DEGREES_TO_RADIANS;
}
/**
* Same as {@link Math#toDegrees(double)} but 3x faster (multiply vs. divide).
* See CompareRadiansSnippet.java in tests.
*/
public static double toDegrees(double radians) {
return radians * RADIANS_TO_DEGREES;
}
}