LUCENE-7158: use the same value (from WGS84) for earth's mean radius when we approximate it as a sphere

This commit is contained in:
Mike McCandless 2016-04-01 05:40:50 -04:00
parent cc75be53f9
commit cf7967cc46
14 changed files with 42 additions and 51 deletions

View File

@ -251,6 +251,10 @@ Bug Fixes
* LUCENE-7126: Remove GeoPointDistanceRangeQuery. This query was implemented
with boolean NOT, and incorrect for multi-valued documents. (Robert Muir)
* LUCENE-7158: Consistently use earth's WGS84 mean radius wherever our
geo search implementations approximate the earth as a sphere (Karl
Wright via Mike McCandless)
Other
* LUCENE-7035: Upgrade icu4j to 56.1/unicode 8. (Robert Muir)

View File

@ -176,8 +176,10 @@ public class SloppyMath {
// TODO: remove these for java 9, they fixed Math.toDegrees()/toRadians() to work just like this.
public static final double TO_RADIANS = Math.PI / 180D;
public static final double TO_DEGREES = 180D / Math.PI;
private static final double TO_METERS = 6_378_137D; // equatorial radius
private static final double TO_KILOMETERS = 6_378.137D; // equatorial radius
// Earth's mean radius, in meters and kilometers; see http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
private static final double TO_METERS = 6_371_008.7714D; // equatorial radius
private static final double TO_KILOMETERS = 6_371.0087714D; // equatorial radius
// cos/asin
private static final double ONE_DIV_F2 = 1/2.0;

View File

@ -99,8 +99,8 @@ public class TestSloppyMath extends LuceneTestCase {
assertEquals(0, haversinMeters(90, -180, 90, 180), 0D);
assertEquals(0, haversinMeters(90, 180, 90, 180), 0D);
// Test half a circle on the equator, using WGS84 equatorial earth radius
double earthRadiusMs = 6_378_137D;
// Test half a circle on the equator, using WGS84 mean earth radius in meters
double earthRadiusMs = 6_371_008.7714;
double halfCircle = earthRadiusMs * Math.PI;
assertEquals(halfCircle, haversinMeters(0, 0, 0, 180), 0D);
@ -111,17 +111,17 @@ public class TestSloppyMath extends LuceneTestCase {
double randomLat2 = 40.65 + (r.nextInt(10) - 5) * 360;
double randomLon2 = -73.95 + (r.nextInt(10) - 5) * 360;
assertEquals(8_581.7047, haversinMeters(randomLat1, randomLon1, randomLat2, randomLon2), 0.01D);
assertEquals(8_572.1137, haversinMeters(randomLat1, randomLon1, randomLat2, randomLon2), 0.01D);
// from solr and ES tests (with their respective epsilons)
assertEquals(0, haversinMeters(40.7143528, -74.0059731, 40.7143528, -74.0059731), 0D);
assertEquals(5_291.80, haversinMeters(40.7143528, -74.0059731, 40.759011, -73.9844722), 0.01D);
assertEquals(462.62, haversinMeters(40.7143528, -74.0059731, 40.718266, -74.007819), 0.01D);
assertEquals(1_056.16, haversinMeters(40.7143528, -74.0059731, 40.7051157, -74.0088305), 0.01D);
assertEquals(1_259.53, haversinMeters(40.7143528, -74.0059731, 40.7247222, -74), 0.01D);
assertEquals(2_030.79, haversinMeters(40.7143528, -74.0059731, 40.731033, -73.9962255), 0.01D);
assertEquals(8_581.70, haversinMeters(40.7143528, -74.0059731, 40.65, -73.95), 0.01D);
assertEquals(5_285.89, haversinMeters(40.7143528, -74.0059731, 40.759011, -73.9844722), 0.01D);
assertEquals(462.10, haversinMeters(40.7143528, -74.0059731, 40.718266, -74.007819), 0.01D);
assertEquals(1_054.98, haversinMeters(40.7143528, -74.0059731, 40.7051157, -74.0088305), 0.01D);
assertEquals(1_258.12, haversinMeters(40.7143528, -74.0059731, 40.7247222, -74), 0.01D);
assertEquals(2_028.52, haversinMeters(40.7143528, -74.0059731, 40.731033, -73.9962255), 0.01D);
assertEquals(8_572.11, haversinMeters(40.7143528, -74.0059731, 40.65, -73.95), 0.01D);
}
/** Test this method sorts the same way as real haversin */
@ -164,6 +164,6 @@ public class TestSloppyMath extends LuceneTestCase {
double h1 = (1 - StrictMath.cos(StrictMath.toRadians(lat2) - StrictMath.toRadians(lat1))) / 2;
double h2 = (1 - StrictMath.cos(StrictMath.toRadians(lon2) - StrictMath.toRadians(lon1))) / 2;
double h = h1 + StrictMath.cos(StrictMath.toRadians(lat1)) * StrictMath.cos(StrictMath.toRadians(lat2)) * h2;
return 2 * 6378137 * StrictMath.asin(Math.min(1, Math.sqrt(h)));
return 2 * 6371008.7714 * StrictMath.asin(Math.min(1, Math.sqrt(h)));
}
}

View File

@ -70,12 +70,13 @@ public class DistanceFacetsExample implements Closeable {
/** The "home" longitude. */
public final static double ORIGIN_LONGITUDE = -74.0059731;
/** Radius of the Earth in KM
/** Mean radius of the Earth in KM
*
* NOTE: this is approximate, because the earth is a bit
* wider at the equator than the poles. See
* http://en.wikipedia.org/wiki/Earth_radius */
public final static double EARTH_RADIUS_KM = 6371.01;
// see http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
public final static double EARTH_RADIUS_KM = 6_371.0087714;
/** Empty constructor */
public DistanceFacetsExample() {}

View File

@ -224,13 +224,13 @@ public class TestDemoExpressions extends LuceneTestCase {
TopFieldDocs td = searcher.search(new MatchAllDocsQuery(), 3, sort);
FieldDoc d = (FieldDoc) td.scoreDocs[0];
assertEquals(0.4626D, (Double)d.fields[0], 1E-4);
assertEquals(0.4621D, (Double)d.fields[0], 1E-4);
d = (FieldDoc) td.scoreDocs[1];
assertEquals(1.0562D, (Double)d.fields[0], 1E-4);
assertEquals(1.055D, (Double)d.fields[0], 1E-4);
d = (FieldDoc) td.scoreDocs[2];
assertEquals(5.2918D, (Double)d.fields[0], 1E-4);
assertEquals(5.2859D, (Double)d.fields[0], 1E-4);
}
public void testStaticExtendedVariableExample() throws Exception {

View File

@ -158,7 +158,7 @@ public class TestJavascriptFunction extends LuceneTestCase {
}
public void testHaversinMethod() throws Exception {
assertEvaluatesTo("haversin(40.7143528,-74.0059731,40.759011,-73.9844722)", 5.291799723323441);
assertEvaluatesTo("haversin(40.7143528,-74.0059731,40.759011,-73.9844722)", 5.285885589128259);
}
public void testLnMethod() throws Exception {

View File

@ -62,13 +62,13 @@ public class TestLatLonPointDistanceSort extends LuceneTestCase {
TopDocs td = searcher.search(new MatchAllDocsQuery(), 3, sort);
FieldDoc d = (FieldDoc) td.scoreDocs[0];
assertEquals(462.6174876948475D, (Double)d.fields[0], 0.0D);
assertEquals(462.1004647449412, (Double)d.fields[0], 0.0D);
d = (FieldDoc) td.scoreDocs[1];
assertEquals(1056.163041670945D, (Double)d.fields[0], 0.0D);
assertEquals(1054.9826700985088, (Double)d.fields[0], 0.0D);
d = (FieldDoc) td.scoreDocs[2];
assertEquals(5291.798081190281D, (Double)d.fields[0], 0.0D);
assertEquals(5285.883948830351, (Double)d.fields[0], 0.0D);
reader.close();
dir.close();
@ -99,10 +99,10 @@ public class TestLatLonPointDistanceSort extends LuceneTestCase {
TopDocs td = searcher.search(new MatchAllDocsQuery(), 3, sort);
FieldDoc d = (FieldDoc) td.scoreDocs[0];
assertEquals(462.6174876948475D, (Double)d.fields[0], 0.0D);
assertEquals(462.1004647449412D, (Double)d.fields[0], 0.0D);
d = (FieldDoc) td.scoreDocs[1];
assertEquals(1056.163041670945D, (Double)d.fields[0], 0.0D);
assertEquals(1054.9826700985088, (Double)d.fields[0], 0.0D);
d = (FieldDoc) td.scoreDocs[2];
assertEquals(Double.POSITIVE_INFINITY, (Double)d.fields[0], 0.0D);

View File

@ -53,8 +53,9 @@ public final class GeoUtils {
public static final double MAX_LAT_RADIANS = TO_RADIANS * MAX_LAT_INCL;
// WGS84 earth-ellipsoid parameters
/** major (a) axis in meters */
public static final double SEMIMAJOR_AXIS = 6_378_137; // [m]
/** mean earth axis in meters */
// see http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
public static final double EARTH_MEAN_RADIUS_METERS = 6_371_008.7714;
// No instance:
private GeoUtils() {
@ -79,7 +80,7 @@ public final class GeoUtils {
final double radLat = TO_RADIANS * centerLat;
final double radLon = TO_RADIANS * centerLon;
// LUCENE-7143
double radDistance = (radiusMeters + 7E-2) / SEMIMAJOR_AXIS;
double radDistance = (radiusMeters + 7E-2) / EARTH_MEAN_RADIUS_METERS;
double minLat = radLat - radDistance;
double maxLat = radLat + radDistance;
double minLon;
@ -129,7 +130,7 @@ public final class GeoUtils {
}
/** maximum error from {@link #axisLat(double, double)}. logic must be prepared to handle this */
public static final double AXISLAT_ERROR = 0.1D / SEMIMAJOR_AXIS * TO_DEGREES;
public static final double AXISLAT_ERROR = 0.1D / EARTH_MEAN_RADIUS_METERS * TO_DEGREES;
/**
* Calculate the latitude of a circle's intersections with its bbox meridians.
@ -156,7 +157,7 @@ public final class GeoUtils {
// the argument to arc cosine, resulting in a range (0, PI/2].
double l1 = TO_RADIANS * centerLat;
double r = (radiusMeters + 7E-2) / SEMIMAJOR_AXIS;
double r = (radiusMeters + 7E-2) / EARTH_MEAN_RADIUS_METERS;
// if we are within radius range of a pole, the lat is the pole itself
if (Math.abs(l1) + r >= MAX_LAT_RADIANS) {

View File

@ -957,7 +957,7 @@ public abstract class BaseGeoPointTestCase extends LuceneTestCase {
radiusMeters = random().nextDouble() * 333000 + 1.0;
} else {
// So the query can cover at most 50% of the earth's surface:
radiusMeters = random().nextDouble() * GeoUtils.SEMIMAJOR_AXIS * Math.PI / 2.0 + 1.0;
radiusMeters = random().nextDouble() * GeoUtils.EARTH_MEAN_RADIUS_METERS * Math.PI / 2.0 + 1.0;
}
if (VERBOSE) {

View File

@ -265,7 +265,7 @@ public class TestGeoUtils extends LuceneTestCase {
}
public void testAxisLat() {
double earthCircumference = 2D * Math.PI * GeoUtils.SEMIMAJOR_AXIS;
double earthCircumference = 2D * Math.PI * GeoUtils.EARTH_MEAN_RADIUS_METERS;
assertEquals(90, GeoUtils.axisLat(0, earthCircumference / 4), 0.0D);
for (int i = 0; i < 100; ++i) {

View File

@ -46,11 +46,8 @@ import org.apache.lucene.util.NumericUtils;
* @lucene.experimental */
public final class Geo3DPoint extends Field {
/** Mean radius of the earth, in meters */
protected final static double MEAN_EARTH_RADIUS_METERS = 6371008.7714;
/** How many radians are in one earth surface meter */
protected final static double RADIANS_PER_METER = 1.0 / MEAN_EARTH_RADIUS_METERS;
protected final static double RADIANS_PER_METER = 1.0 / PlanetModel.WGS84_MEAN;
/** Indexing {@link FieldType}. */
public static final FieldType TYPE = new FieldType();

View File

@ -26,7 +26,8 @@ public class PlanetModel {
public static final PlanetModel SPHERE = new PlanetModel(1.0,1.0);
/** Mean radius */
public static final double WGS84_MEAN = 6371009.0;
// see http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
public static final double WGS84_MEAN = 6371008.7714;
/** Polar radius */
public static final double WGS84_POLAR = 6356752.314245;
/** Equatorial radius */

View File

@ -796,7 +796,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
public void testToString() {
Geo3DPoint point = new Geo3DPoint("point", 44.244272, 7.769736);
assertEquals("Geo3DPoint <point: x=0.709426287693908 y=0.09679758561541502 z=0.6973564369288621>", point.toString());
assertEquals("Geo3DPoint <point: x=0.709426313149037 y=0.09679758908863707 z=0.6973564619509093>", point.toString());
}
public void testShapeQueryToString() {

View File

@ -211,21 +211,6 @@ public class GeoCircleTest extends LuceneTestCase {
xyzb.getMinimumX(), xyzb.getMaximumX(), xyzb.getMinimumY(), xyzb.getMaximumY(), xyzb.getMinimumZ(), xyzb.getMaximumZ());
assertTrue(GeoArea.WITHIN == area.getRelationship(c) || GeoArea.OVERLAPS == area.getRelationship(c));
// Yet another test case from BKD
c = GeoCircleFactory.makeGeoCircle(PlanetModel.WGS84, 0.006229478708446979, 0.005570196723795424, 3.840276763694387E-5);
xyzb = new XYZBounds();
c.getBounds(xyzb);
area = GeoAreaFactory.makeGeoArea(PlanetModel.WGS84,
xyzb.getMinimumX(), xyzb.getMaximumX(), xyzb.getMinimumY(), xyzb.getMaximumY(), xyzb.getMinimumZ(), xyzb.getMaximumZ());
p1 = new GeoPoint(PlanetModel.WGS84, 0.006224927111830945, 0.005597367237251763);
p2 = new GeoPoint(1.0010836083810235, 0.005603490759433942, 0.006231850560862502);
assertTrue(PlanetModel.WGS84.pointOnSurface(p1));
//assertTrue(PlanetModel.WGS84.pointOnSurface(p2));
assertTrue(c.isWithin(p1));
assertTrue(c.isWithin(p2));
assertTrue(area.isWithin(p1));
assertTrue(area.isWithin(p2));
// Another test case from BKD
c = GeoCircleFactory.makeGeoCircle(PlanetModel.SPHERE, -0.005955031040627789, -0.0029274772647399153, 1.601488279374338E-5);
xyzb = new XYZBounds();