SOLR-2125: change interpretation of bounding box to better fit user expectations on accuracy. fix issues with radians. add more tests

git-svn-id: https://svn.apache.org/repos/asf/lucene/dev/trunk@1000428 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Grant Ingersoll 2010-09-23 12:29:40 +00:00
parent f68b4b4adb
commit 3f547a822e
4 changed files with 151 additions and 48 deletions

View File

@ -34,11 +34,15 @@ public class DistanceUtils {
public static final double DEGREES_TO_RADIANS = Math.PI / 180.0;
public static final double RADIANS_TO_DEGREES = 180.0 / Math.PI;
public static final double DEG_45 = Math.PI / 4.0;
public static final double DEG_225 = 5 * DEG_45;
public static final double DEG_90 = Math.PI / 2;
public static final double DEG_180 = Math.PI;
public static final double SIN_45 = Math.sin(DEG_45);
//pre-compute some angles that are commonly used
public static final double DEG_45_AS_RADS = Math.PI / 4.0;
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;
public static final double DEG_225_AS_RADS = 5 * DEG_45_AS_RADS;
public static final double DEG_270_AS_RADS = 3*DEG_90_AS_RADS;
public static final double KM_TO_MILES = 0.621371192;
public static final double MILES_TO_KM = 1.609344;
/**
@ -162,7 +166,7 @@ public class DistanceUtils {
//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 * distance; // sin(Pi/4) == (2^0.5)/2 == opp/hyp == opp/distance, solve for opp, similarily for cosine
distance = SIN_45_AS_RADS * distance; // sin(Pi/4) == (2^0.5)/2 == opp/hyp == opp/distance, solve for opp, similarily for cosine
for (int i = 0; i < center.length; i++) {
result[i] = center[i] + distance;
}
@ -175,41 +179,68 @@ public class DistanceUtils {
* @param distance The distance
* @param result A preallocated array to hold the results. If null, a new one is constructed.
* @param upperRight If true, calculate the upper right corner, else the lower left
* @param radius The radius of the sphere to use.
* @param sphereRadius The radius of the sphere to use.
* @return The Lat/Lon in degrees
*
* @see #latLonCorner(double, double, double, double[], boolean, double)
*/
public static double[] latLonCornerDegs(double latCenter, double lonCenter,
double distance, double [] result,
boolean upperRight, double radius) {
boolean upperRight, double sphereRadius) {
result = latLonCorner(latCenter * DEGREES_TO_RADIANS,
lonCenter * DEGREES_TO_RADIANS, distance, result, upperRight, radius);
lonCenter * DEGREES_TO_RADIANS, distance, result, upperRight, sphereRadius);
result[0] = result[0] * RADIANS_TO_DEGREES;
result[1] = result[1] * RADIANS_TO_DEGREES;
return result;
}
/**
* Uses Haversine to calculate the corner
* Uses Haversine to calculate the corner of a box (upper right or lower left) that is the <i>distance</i> away, given a sphere of the specified <i>radius</i>.
*
* NOTE: This is not the same as calculating a box that transcribes a circle of the given distance.
*
* @param latCenter In radians
* @param lonCenter In radians
* @param distance The distance
* @param result A preallocated array to hold the results. If null, a new one is constructed.
* @param upperRight If true, give lat/lon for the upper right corner, else lower left
* @param radius The radius to use for the calculation
* @param sphereRadius The radius to use for the calculation
* @return The Lat/Lon in Radians
*/
public static double[] latLonCorner(double latCenter, double lonCenter,
double distance, double [] result, boolean upperRight, double radius) {
double distance, double [] result, boolean upperRight, double sphereRadius) {
// Haversine formula
double brng = upperRight ? DEG_45 : DEG_225;
double lat2 = Math.asin(Math.sin(latCenter) * Math.cos(distance / radius) +
Math.cos(latCenter) * Math.sin(distance / radius) * Math.cos(brng));
double lon2 = lonCenter + Math.atan2(Math.sin(brng) * Math.sin(distance / radius) * Math.cos(latCenter),
Math.cos(distance / radius) - Math.sin(latCenter) * Math.sin(lat2));
double brng = upperRight ? DEG_45_AS_RADS : DEG_225_AS_RADS;
result = pointOnBearing(latCenter, lonCenter, distance, brng, result, sphereRadius);
return result;
}
/**
* Given a start point (startLat, startLon) and a bearing on a sphere of radius <i>sphereRadius</i>, return the destination point.
* @param startLat The starting point latitude, in radians
* @param startLon The starting point longitude, in radians
* @param distance The distance to travel along the bearing. The units are assumed to be the same as the sphereRadius units, both of which is up to the caller to know
* @param bearing The bearing, in radians. North is a 0 deg. bearing, east is 90 deg, south is 180 deg, west is 270 deg.
* @param result A preallocated array to hold the results. If null, a new one is constructed.
* @param sphereRadius The radius of the sphere to use for the calculation.
* @return The destination point, in radians. First entry is latitude, second is longitude
*/
public static double[] pointOnBearing(double startLat, double startLon, double distance, double bearing, double[] result, double sphereRadius) {
/*
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(distance / sphereRadius);
double cosStartLat = Math.cos(startLat);
double sinAngDist = Math.sin(distance / sphereRadius);
double lat2 = Math.asin(Math.sin(startLat) * cosAngDist +
cosStartLat * sinAngDist * Math.cos(bearing));
double lon2 = startLon + Math.atan2(Math.sin(bearing) * sinAngDist * cosStartLat,
cosAngDist - Math.sin(startLat) * Math.sin(lat2));
/*lat2 = (lat2*180)/Math.PI;
lon2 = (lon2*180)/Math.PI;*/
@ -224,7 +255,6 @@ public class DistanceUtils {
// normalize lat - could flip poles
normLat(result);
return result;
}
@ -233,19 +263,19 @@ public class DistanceUtils {
*/
public static void normLat(double[] latLng) {
if (latLng[0] > DEG_90) {
latLng[0] = DEG_90 - (latLng[0] - DEG_90);
if (latLng[0] > DEG_90_AS_RADS) {
latLng[0] = DEG_90_AS_RADS - (latLng[0] - DEG_90_AS_RADS);
if (latLng[1] < 0) {
latLng[1] = latLng[1] + DEG_180;
latLng[1] = latLng[1] + DEG_180_AS_RADS;
} else {
latLng[1] = latLng[1] - DEG_180;
latLng[1] = latLng[1] - DEG_180_AS_RADS;
}
} else if (latLng[0] < -DEG_90) {
latLng[0] = -DEG_90 - (latLng[0] + DEG_90);
} else if (latLng[0] < -DEG_90_AS_RADS) {
latLng[0] = -DEG_90_AS_RADS - (latLng[0] + DEG_90_AS_RADS);
if (latLng[1] < 0) {
latLng[1] = latLng[1] + DEG_180;
latLng[1] = latLng[1] + DEG_180_AS_RADS;
} else {
latLng[1] = latLng[1] - DEG_180;
latLng[1] = latLng[1] - DEG_180_AS_RADS;
}
}
@ -257,10 +287,10 @@ public class DistanceUtils {
* @param latLng The lat/lon, in radians, lat in position 0, long in position 1
*/
public static void normLng(double[] latLng) {
if (latLng[1] > DEG_180) {
latLng[1] = -1.0 * (DEG_180 - (latLng[1] - DEG_180));
} else if (latLng[1] < -DEG_180) {
latLng[1] = (latLng[1] + DEG_180) + DEG_180;
if (latLng[1] > DEG_180_AS_RADS) {
latLng[1] = -1.0 * (DEG_180_AS_RADS - (latLng[1] - DEG_180_AS_RADS));
} else if (latLng[1] < -DEG_180_AS_RADS) {
latLng[1] = (latLng[1] + DEG_180_AS_RADS) + DEG_180_AS_RADS;
}
}

View File

@ -64,7 +64,7 @@ public class DistanceUtilsTest extends LuceneTestCase {
public void testLatLonCorner() throws Exception {
double[] zero = new double[]{0, 0};
double[] zero45 = new double[]{0, DistanceUtils.DEG_45};
double[] zero45 = new double[]{0, DistanceUtils.DEG_45_AS_RADS};
double[] result;
// 00°3809N, 000°3809E
//Verify at http://www.movable-type.co.uk/scripts/latlong.html
@ -94,8 +94,59 @@ public class DistanceUtilsTest extends LuceneTestCase {
}
public void testPointBearing() throws Exception {
double[] zero = new double[]{0, 0};
double[] zero45 = new double[]{40 * DistanceUtils.DEGREES_TO_RADIANS, DistanceUtils.DEG_45_AS_RADS};
double[] result;
// 00°3809N, 000°3809E
//Verify at http://www.movable-type.co.uk/scripts/latlong.html
result = DistanceUtils.pointOnBearing(zero[0], zero[1], 100, DistanceUtils.DEG_45_AS_RADS, null, DistanceUtils.EARTH_MEAN_RADIUS_KM);
assertEquals(0.63583 * DistanceUtils.DEGREES_TO_RADIANS, result[0], 0.001);
assertEquals(0.63583 * DistanceUtils.DEGREES_TO_RADIANS, result[1], 0.001);
//should be above the current point at 0.8994°,0.0000°
result = DistanceUtils.pointOnBearing(zero[0], zero[1], 100, 0, null, DistanceUtils.EARTH_MEAN_RADIUS_KM);
assertEquals(0.8994 * DistanceUtils.DEGREES_TO_RADIANS, result[0], 0.001);
assertEquals(0, result[1], 0.001);
//directly below
result = DistanceUtils.pointOnBearing(zero[0], zero[1], 100, DistanceUtils.DEG_180_AS_RADS, null, DistanceUtils.EARTH_MEAN_RADIUS_KM);
assertEquals(-0.8994 * DistanceUtils.DEGREES_TO_RADIANS, result[0], 0.001);
assertEquals(0, result[1], 0.001);
//0.7183°,0.5414° -- 37 deg bearing
result = DistanceUtils.pointOnBearing(zero[0], zero[1], 100, 37 * DistanceUtils.DEGREES_TO_RADIANS, null, DistanceUtils.EARTH_MEAN_RADIUS_KM);
assertEquals(0.7183 * DistanceUtils.DEGREES_TO_RADIANS, result[0], 0.001);
assertEquals(0.5414 * DistanceUtils.DEGREES_TO_RADIANS, result[1], 0.001);
result = DistanceUtils.pointOnBearing(zero45[0], zero45[1], 100, DistanceUtils.DEG_45_AS_RADS, null, DistanceUtils.EARTH_MEAN_RADIUS_KM);
//40.6328°,45.8381°
assertEquals(40.6328 * DistanceUtils.DEGREES_TO_RADIANS, result[0], 0.001);
assertEquals(45.8381 * DistanceUtils.DEGREES_TO_RADIANS, result[1], 0.001);
result = DistanceUtils.pointOnBearing(1 * DistanceUtils.DEGREES_TO_RADIANS, 1 * DistanceUtils.DEGREES_TO_RADIANS, 100, DistanceUtils.DEG_90_AS_RADS, null, DistanceUtils.EARTH_MEAN_RADIUS_KM);
//0.9997°,1.8994°
assertEquals(0.9997 * DistanceUtils.DEGREES_TO_RADIANS, result[0], 0.001);
assertEquals(1.8994 * DistanceUtils.DEGREES_TO_RADIANS, result[1], 0.001);
result = DistanceUtils.pointOnBearing(-10 * DistanceUtils.DEGREES_TO_RADIANS, -150 * DistanceUtils.DEGREES_TO_RADIANS, 15, 205*DistanceUtils.DEGREES_TO_RADIANS, null, DistanceUtils.EARTH_MEAN_RADIUS_KM);
//-10.1222°,-150.0578°
assertEquals(-10.1222 * DistanceUtils.DEGREES_TO_RADIANS, result[0], 0.001);
assertEquals(-150.0578 * DistanceUtils.DEGREES_TO_RADIANS, result[1], 0.001);
result = DistanceUtils.pointOnBearing(-10 * DistanceUtils.DEGREES_TO_RADIANS, -150 * DistanceUtils.DEGREES_TO_RADIANS, 200, 63*DistanceUtils.DEGREES_TO_RADIANS, null, DistanceUtils.EARTH_MEAN_RADIUS_KM);
//-9.1797°,-148.3767°
assertEquals(-9.1797 * DistanceUtils.DEGREES_TO_RADIANS, result[0], 0.001);
assertEquals(-148.3767 * DistanceUtils.DEGREES_TO_RADIANS, result[1], 0.001);
result = DistanceUtils.pointOnBearing(-10 * DistanceUtils.DEGREES_TO_RADIANS, -150 * DistanceUtils.DEGREES_TO_RADIANS, 3000, 63*DistanceUtils.DEGREES_TO_RADIANS, null, DistanceUtils.EARTH_MEAN_RADIUS_KM);
//2.7561°,-126.1281°
assertEquals(2.7561 * DistanceUtils.DEGREES_TO_RADIANS, result[0], 0.001);
assertEquals(-126.1281 * DistanceUtils.DEGREES_TO_RADIANS, result[1], 0.001);
}
public void testVectorDistance() throws Exception {
double[] zero = new double[]{0, 0};
double[] zeroOne = new double[]{0, 1};
double[] oneZero = new double[]{1, 0};
double[] oneOne = new double[]{1, 1};

View File

@ -91,26 +91,34 @@ public class LatLonType extends AbstractSubTypeFieldType implements SpatialQuery
} catch (InvalidGeoException e) {
throw new SolrException(SolrException.ErrorCode.BAD_REQUEST, e);
}
point[0] = point[0] * DistanceUtils.DEGREES_TO_RADIANS;
point[1] = point[1] * DistanceUtils.DEGREES_TO_RADIANS;
//Get the distance
double[] ur;
double[] ll;
if (options.measStr == null || options.measStr.equals("hsin")) {
ur = DistanceUtils.latLonCornerDegs(point[LAT], point[LONG], options.distance, null, true, options.radius);
ll = DistanceUtils.latLonCornerDegs(point[LAT], point[LONG], options.distance, null, false, options.radius);
} else {
ur = DistanceUtils.vectorBoxCorner(point, null, options.distance, true);
ll = DistanceUtils.vectorBoxCorner(point, null, options.distance, false);
}
double[] ur = new double[2];
double[] ll = new double[2];
double[] tmp = new double[2];
//these calculations aren't totally accurate, but it should be good enough
//TODO: Optimize to do in single calculations. Would need to deal with poles, prime meridian, etc.
double [] north = DistanceUtils.pointOnBearing(point[LAT], point[LONG], options.distance, 0, tmp, options.radius);
//This returns the point as radians, but we need degrees b/c that is what the field is stored as
ur[LAT] = north[LAT] * DistanceUtils.RADIANS_TO_DEGREES;//get it now, as we are going to reuse tmp
double [] east = DistanceUtils.pointOnBearing(point[LAT], point[LONG], options.distance, DistanceUtils.DEG_90_AS_RADS, tmp, options.radius);
ur[LONG] = east[LONG] * DistanceUtils.RADIANS_TO_DEGREES;
double [] south = DistanceUtils.pointOnBearing(point[LAT], point[LONG], options.distance, DistanceUtils.DEG_180_AS_RADS, tmp, options.radius);
ll[LAT] = south[LAT] * DistanceUtils.RADIANS_TO_DEGREES;
double [] west = DistanceUtils.pointOnBearing(point[LAT], point[LONG], options.distance, DistanceUtils.DEG_270_AS_RADS, tmp, options.radius);
ll[LONG] = west[LONG] * DistanceUtils.RADIANS_TO_DEGREES;
SchemaField subSF;
Query range;
//TODO: can we reuse our bearing calculations?
double angDist = DistanceUtils.angularDistance(options.distance,
options.radius);//in radians
double angDistDegs = DistanceUtils.angularDistance(options.distance,
options.radius) * DistanceUtils.RADIANS_TO_DEGREES;
//for the poles, do something slightly different
if (point[LAT] + angDistDegs > 90.0) { //we cross the north pole
//Also, note point[LAT] is in radians, but ur and ll are in degrees
if (point[LAT] + angDist > DistanceUtils.DEG_90_AS_RADS) { //we cross the north pole
//we don't need a longitude boundary at all
double minLat = Math.min(ll[LAT], ur[LAT]);
@ -119,7 +127,7 @@ public class LatLonType extends AbstractSubTypeFieldType implements SpatialQuery
String.valueOf(minLat),
"90", true, true);
result.add(range, BooleanClause.Occur.MUST);
} else if (point[LAT] - angDistDegs < -90.0) {//we cross the south pole
} else if (point[LAT] - angDist < -DistanceUtils.DEG_90_AS_RADS) {//we cross the south pole
subSF = subField(options.field, LAT);
double maxLat = Math.max(ll[LAT], ur[LAT]);
range = subSF.getType().getRangeQuery(parser, subSF,

View File

@ -93,9 +93,23 @@ public class SpatialFilterTest extends SolrTestCaseJ4 {
checkHits(fieldName, "33.0,-80.0", 300, 2);
//large distance
checkHits(fieldName, "1,1", 5000, 3, 5, 6, 7);
//Try alternate distance
checkHits(fieldName, "0.1,0.1", 15, 1, 6);
//Because we are generating a box based on the west/east longitudes and the south/north latitudes, which then
//translates to a range query, which is slightly more inclusive. Thus, even though 0.0 is 15.725 kms away,
//it will be included, b/c of the box calculation.
checkHits(fieldName, "0.1,0.1", 15, 2, 5, 6);
//try some more
clearIndex();
assertU(adoc("id", "14", fieldName, "0,5"));
assertU(adoc("id", "15", fieldName, "0,15"));
//3000KM from 0,0, see http://www.movable-type.co.uk/scripts/latlong.html
assertU(adoc("id", "16", fieldName, "18.71111,19.79750"));
assertU(commit());
checkHits(fieldName, "0,0", 1000, 1, 14);
checkHits(fieldName, "0,0", 2000, 2, 14, 15);
checkHits(fieldName, "0,0", 3000, 3, 14, 15, 16);
checkHits(fieldName, "0,0", 3001, 3, 14, 15, 16);
checkHits(fieldName, "0,0", 3000.1, 3, 14, 15, 16);
}
private void checkHits(String fieldName, String pt, double distance, int count, int ... docIds) {