/*******************************************************************************
* Copyright (c) 2015 MITRE and VoyagerSearch
* All rights reserved. This program and the accompanying materials
* are made available under the terms of the Apache License, Version 2.0 which
* accompanies this distribution and is available at
* http://www.apache.org/licenses/LICENSE-2.0.txt
******************************************************************************/
package org.locationtech.spatial4j.distance;
import com.carrotsearch.randomizedtesting.RandomizedTest;
import org.locationtech.spatial4j.context.SpatialContext;
import org.locationtech.spatial4j.shape.Circle;
import org.locationtech.spatial4j.shape.Point;
import org.locationtech.spatial4j.shape.Rectangle;
import org.locationtech.spatial4j.shape.SpatialRelation;
import org.locationtech.spatial4j.shape.impl.PointImpl;
import org.junit.Before;
import org.junit.Test;
import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertTrue;
import static org.locationtech.spatial4j.distance.DistanceUtils.DEG_TO_KM;
import static org.locationtech.spatial4j.distance.DistanceUtils.KM_TO_DEG;
public class TestDistances extends RandomizedTest {
//NOTE! These are sometimes modified by tests.
private SpatialContext ctx;
private double EPS;
@Before
public void beforeTest() {
ctx = SpatialContext.GEO;
EPS = 10e-4;//delta when doing double assertions. Geo eps is not that small.
}
private DistanceCalculator dc() {
return ctx.getDistCalc();
}
@Test
public void testSomeDistances() {
//See to verify: from http://www.movable-type.co.uk/scripts/latlong.html
Point ctr = pLL(0,100);
assertEquals(11100, dc().distance(ctr, pLL(10, 0)) * DEG_TO_KM, 3);
double deg = dc().distance(ctr, pLL(10, -160));
assertEquals(11100, deg * DEG_TO_KM, 3);
assertEquals(314.40338, dc().distance(pLL(1, 2), pLL(3, 4)) * DEG_TO_KM, EPS);
}
@Test
public void testDistancesAgainstVincenty() {
DistanceCalculator vincenty = new GeodesicSphereDistCalc.Vincenty();
DistanceCalculator haversine = new GeodesicSphereDistCalc.Haversine();
DistanceCalculator lawOfCos = new GeodesicSphereDistCalc.LawOfCosines();
final int TRIES = 100000 * (int)multiplier();
for (int i = 0; i < TRIES; i++) {
Point p1 = randomGeoPoint();
Point p2 = randomGeoPointFrom(p1);
double distV = vincenty.distance(p1, p2);
//Haversine: accurate to a centimeter if on same side of globe,
// otherwise possibly 1m apart (antipodal)
double havV = haversine.distance(p1, p2);
assertEquals(distV, havV, (distV <= 90) ? DistanceUtils.KM_TO_DEG * 0.00001 : DistanceUtils.KM_TO_DEG * 0.001);
// Fractionally compared to truth, also favorably accurate.
if (distV != 0 && distV > 0.0000001)
assertEquals(1.0, havV/distV, 0.001);//0.1%
//LawOfCosines: accurate to within 1 meter (or better?)
double locV = lawOfCos.distance(p1, p2);
assertEquals(distV, locV, DistanceUtils.KM_TO_DEG * 0.001);
}
}
private Point randomGeoPoint() {
//not uniformly distributed but that's ok
return ctx.makePoint(randomDouble()*360 + -180, randomDouble()*180 + -90);
}
private Point randomGeoPointFrom(Point p1) {
int which = randomInt(10);//inclusive
double distDEG;
if (which <= 2) {
distDEG = 180 - randomDouble() * 0.001 / Math.pow(10, which);
} else if (which >= 8) {
distDEG = randomDouble() * 0.001 / Math.pow(10, 10-which);
} else {
distDEG = randomDouble()*180;
}
double bearingDEG = randomDouble() * 360;
Point p2RAD = DistanceUtils.pointOnBearingRAD(DistanceUtils.toRadians(p1.getY()), DistanceUtils.toRadians(p1.getX()),
DistanceUtils.toRadians(distDEG), DistanceUtils.toRadians(bearingDEG), ctx, null);
p2RAD.reset(DistanceUtils.toDegrees(p2RAD.getX()), DistanceUtils.toDegrees(p2RAD.getY()));
return p2RAD;//now it's in degrees
}
@Test /** See #81 */
public void testHaversineNaN() {
assertEquals(180, new GeodesicSphereDistCalc.Haversine().distance(
ctx.makePoint(-81.05206968336057, 71.82629271026536),
98.9479297952497, -71.82629264390964),
0.00001);
}
@Test
public void testCalcBoxByDistFromPt() {
//first test regression
{
double d = 6894.1 * KM_TO_DEG;
Point pCtr = pLL(-20, 84);
Point pTgt = pLL(-42, 15);
assertTrue(dc().distance(pCtr, pTgt) < d);
//since the pairwise distance is less than d, a bounding box from ctr with d should contain pTgt.
Rectangle r = dc().calcBoxByDistFromPt(pCtr, d, ctx, null);
assertEquals(SpatialRelation.CONTAINS,r.relate(pTgt));
checkBBox(pCtr,d);
}
assertEquals("0 dist, horiz line",
-45,dc().calcBoxByDistFromPt_yHorizAxisDEG(ctx.makePoint(-180, -45), 0, ctx),0);
double MAXDIST = (double) 180 * DEG_TO_KM;
checkBBox(ctx.makePoint(0,0), MAXDIST);
checkBBox(ctx.makePoint(0,0), MAXDIST *0.999999);
checkBBox(ctx.makePoint(0,0),0);
checkBBox(ctx.makePoint(0,0),0.000001);
checkBBox(ctx.makePoint(0,90),0.000001);
checkBBox(ctx.makePoint(-32.7,-5.42),9829);
checkBBox(ctx.makePoint(0,90-20), (double) 20 * DEG_TO_KM);
{
double d = 0.010;//10m
checkBBox(ctx.makePoint(0,90- (d + 0.001) * KM_TO_DEG),d);
}
for (int T = 0; T < 100; T++) {
double lat = -90 + randomDouble()*180;
double lon = -180 + randomDouble()*360;
Point ctr = ctx.makePoint(lon, lat);
double dist = MAXDIST*randomDouble();
checkBBox(ctr, dist);
}
}
private void checkBBox(Point ctr, double distKm) {
String msg = "ctr: "+ctr+" distKm: "+distKm;
double dist = distKm * KM_TO_DEG;
Rectangle r = dc().calcBoxByDistFromPt(ctr, dist, ctx, null);
double horizAxisLat = dc().calcBoxByDistFromPt_yHorizAxisDEG(ctr, dist, ctx);
if (!Double.isNaN(horizAxisLat))
assertTrue(r.relateYRange(horizAxisLat, horizAxisLat).intersects());
//horizontal
if (r.getWidth() >= 180) {
double deg = dc().distance(ctr, r.getMinX(), r.getMaxY() == 90 ? 90 : -90);
double calcDistKm = deg * DEG_TO_KM;
assertTrue(msg, calcDistKm <= distKm + EPS);
//horizAxisLat is meaningless in this context
} else {
Point tPt = findClosestPointOnVertToPoint(r.getMinX(), r.getMinY(), r.getMaxY(), ctr);
double calcDistKm = dc().distance(ctr, tPt) * DEG_TO_KM;
assertEquals(msg, distKm, calcDistKm, EPS);
assertEquals(msg, tPt.getY(), horizAxisLat, EPS);
}
//vertical
double topDistKm = dc().distance(ctr, ctr.getX(), r.getMaxY()) * DEG_TO_KM;
if (r.getMaxY() == 90)
assertTrue(msg, topDistKm <= distKm + EPS);
else
assertEquals(msg, distKm, topDistKm, EPS);
double botDistKm = dc().distance(ctr, ctr.getX(), r.getMinY()) * DEG_TO_KM;
if (r.getMinY() == -90)
assertTrue(msg, botDistKm <= distKm + EPS);
else
assertEquals(msg, distKm, botDistKm, EPS);
}
private Point findClosestPointOnVertToPoint(double lon, double lowLat, double highLat, Point ctr) {
//A binary search algorithm to find the point along the vertical lon between lowLat & highLat that is closest
// to ctr, and returns the distance.
double midLat = (highLat - lowLat)/2 + lowLat;
double midLatDist = ctx.getDistCalc().distance(ctr,lon,midLat);
for(int L = 0; L < 100 && (highLat - lowLat > 0.001|| L < 20); L++) {
boolean bottom = (midLat - lowLat > highLat - midLat);
double newMid = bottom ? (midLat - lowLat)/2 + lowLat : (highLat - midLat)/2 + midLat;
double newMidDist = ctx.getDistCalc().distance(ctr,lon,newMid);
if (newMidDist < midLatDist) {
if (bottom) {
highLat = midLat;
} else {
lowLat = midLat;
}
midLat = newMid;
midLatDist = newMidDist;
} else {
if (bottom) {
lowLat = newMid;
} else {
highLat = newMid;
}
}
}
return ctx.makePoint(lon,midLat);
}
@Test
public void testDistCalcPointOnBearing_cartesian() {
ctx = new SpatialContext(false);
EPS = 10e-6;//tighter epsilon (aka delta)
for(int i = 0; i < 1000; i++) {
testDistCalcPointOnBearing(randomInt(100));
}
}
@Test
public void testDistCalcPointOnBearing_geo() {
//The haversine formula has a higher error if the points are near antipodal. We adjust EPS tolerance for this case.
//TODO Eventually we should add the Vincenty formula for improved accuracy, or try some other cleverness.
//test known high delta
// {
// Point c = ctx.makePoint(-103,-79);
// double angRAD = Math.toRadians(236);
// double dist = 20025;
// Point p2 = dc().pointOnBearingRAD(c, dist, angRAD, ctx);
// //Pt(x=76.61200011750923,y=79.04946929870962)
// double calcDist = dc().distance(c, p2);
// assertEqualsRatio(dist, calcDist);
// }
double maxDistKm = (double) 180 * DEG_TO_KM;
for(int i = 0; i < 1000; i++) {
int distKm = randomInt((int) maxDistKm);
EPS = (distKm < maxDistKm*0.75 ? 10e-6 : 10e-3);
testDistCalcPointOnBearing(distKm);
}
}
private void testDistCalcPointOnBearing(double distKm) {
for(int angDEG = 0; angDEG < 360; angDEG += randomIntBetween(1,20)) {
Point c = ctx.makePoint(
DistanceUtils.normLonDEG(randomInt(359)),
randomIntBetween(-90,90));
//0 distance means same point
Point p2 = dc().pointOnBearing(c, 0, angDEG, ctx, null);
assertEquals(c,p2);
p2 = dc().pointOnBearing(c, distKm * KM_TO_DEG, angDEG, ctx, null);
double calcDistKm = dc().distance(c, p2) * DEG_TO_KM;
assertEqualsRatio(distKm, calcDistKm);
}
}
private void assertEqualsRatio(double expected, double actual) {
double delta = Math.abs(actual - expected);
double base = Math.min(actual, expected);
double deltaRatio = base==0 ? delta : Math.min(delta,delta / base);
assertEquals(0,deltaRatio, EPS);
}
@Test
public void testNormLat() {
double[][] lats = new double[][] {
{1.23,1.23},//1.23 might become 1.2299999 after some math and we want to ensure that doesn't happen
{-90,-90},{90,90},{0,0}, {-100,-80},
{-90-180,90},{-90-360,-90},{90+180,-90},{90+360,90},
{-12+180,12}};
for (double[] pair : lats) {
assertEquals("input "+pair[0], pair[1], DistanceUtils.normLatDEG(pair[0]), 0);
}
for(int i = -1000; i < 1000; i += randomInt(9)*10) {
double d = DistanceUtils.normLatDEG(i);
assertTrue(i + " " + d, d >= -90 && d <= 90);
}
}
@Test
public void testNormLon() {
double[][] lons = new double[][] {
{1.23,1.23},//1.23 might become 1.2299999 after some math and we want to ensure that doesn't happen
{-180,-180},{180,+180},{0,0}, {-190,170},{181,-179},
{-180-360,-180},{-180-720,-180},
{180+360,+180},{180+720,+180}};
for (double[] pair : lons) {
assertEquals("input " + pair[0], pair[1], DistanceUtils.normLonDEG(pair[0]), 0);
}
for(int i = -1000; i < 1000; i += randomInt(9)*10) {
double d = DistanceUtils.normLonDEG(i);
assertTrue(i + " " + d, d >= -180 && d <= 180);
}
}
@Test
public void assertDistanceConversion() {
assertDistanceConversion(0);
assertDistanceConversion(500);
assertDistanceConversion(DistanceUtils.EARTH_MEAN_RADIUS_KM);
}
private void assertDistanceConversion(double dist) {
double radius = DistanceUtils.EARTH_MEAN_RADIUS_KM;
//test back & forth conversion for both
double distRAD = DistanceUtils.dist2Radians(dist, radius);
assertEquals(dist, DistanceUtils.radians2Dist(distRAD, radius), EPS);
double distDEG = DistanceUtils.dist2Degrees(dist, radius);
assertEquals(dist, DistanceUtils.degrees2Dist(distDEG, radius), EPS);
//test across rad & deg
assertEquals(distDEG,DistanceUtils.toDegrees(distRAD),EPS);
//test point on bearing
assertEquals(
DistanceUtils.pointOnBearingRAD(0, 0, DistanceUtils.dist2Radians(dist, radius), DistanceUtils.DEG_90_AS_RADS, ctx, new PointImpl(0, 0, ctx)).getX(),
distRAD, 10e-5);
}
private Point pLL(double lat, double lon) {
return ctx.makePoint(lon,lat);
}
@Test
public void testArea() {
double radius = DistanceUtils.EARTH_MEAN_RADIUS_KM * KM_TO_DEG;
//surface of a sphere is 4 * pi * r^2
final double earthArea = 4 * Math.PI * radius * radius;
Circle c = ctx.makeCircle(randomIntBetween(-180,180), randomIntBetween(-90,90),
180);//180 means whole earth
assertEquals(earthArea, c.getArea(ctx), 1.0);
assertEquals(earthArea, ctx.getWorldBounds().getArea(ctx), 1.0);
//now check half earth
Circle cHalf = ctx.makeCircle(c.getCenter(), 90);
assertEquals(earthArea/2, cHalf.getArea(ctx), 1.0);
//circle with same radius at +20 lat with one at -20 lat should have same area as well as bbox with same area
Circle c2 = ctx.makeCircle(c.getCenter(), 30);
Circle c3 = ctx.makeCircle(c.getCenter().getX(), 20, 30);
assertEquals(c2.getArea(ctx), c3.getArea(ctx), 0.01);
Circle c3Opposite = ctx.makeCircle(c.getCenter().getX(), -20, 30);
assertEquals(c3.getArea(ctx), c3Opposite.getArea(ctx), 0.01);
assertEquals(c3.getBoundingBox().getArea(ctx), c3Opposite.getBoundingBox().getArea(ctx), 0.01);
//small shapes near the equator should have similar areas to euclidean rectangle
Rectangle smallRect = ctx.makeRectangle(0, 1, 0, 1);
assertEquals(1.0, smallRect.getArea(null), 0.0);
double smallDelta = smallRect.getArea(null) - smallRect.getArea(ctx);
assertTrue(smallDelta > 0 && smallDelta < 0.0001);
Circle smallCircle = ctx.makeCircle(0,0,1);
smallDelta = smallCircle.getArea(null) - smallCircle.getArea(ctx);
assertTrue(smallDelta > 0 && smallDelta < 0.0001);
//bigger, but still fairly similar
//c2 = ctx.makeCircle(c.getCenter(), 30);
double areaRatio = c2.getArea(null) / c2.getArea(ctx);
assertTrue(areaRatio > 1 && areaRatio < 1.1);
}
}