mirror of https://github.com/apache/lucene.git
LUCENE-6908: Fix TestGeoUtils.testGeoRelations to handle irregular rectangles. Refactors relational methods to new GeoRelationUtils class.
git-svn-id: https://svn.apache.org/repos/asf/lucene/dev/trunk@1718481 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Nick Knize
parent
6d8a195822
commit
92b659ab9e
|
|
@ -116,6 +116,10 @@ New Features
|
|||
|
||||
API Changes
|
||||
|
||||
* LUCENE-6908: GeoUtils static relational methods have been refactored to new
|
||||
GeoRelationUtils and now correctly handle large irregular rectangles, and
|
||||
pole crossing distance queries. (Nick Knize)
|
||||
|
||||
* LUCENE-6900: Grouping sortWithinGroup variables used to allow null to mean
|
||||
Sort.RELEVANCE. Null is no longer permitted. (David Smiley)
|
||||
|
||||
|
|
|
|||
|
|
@ -87,7 +87,41 @@ public class SloppyMath {
|
|||
double indexSin = sinTab[index];
|
||||
return indexCos + delta * (-indexSin + delta * (-indexCos * ONE_DIV_F2 + delta * (indexSin * ONE_DIV_F3 + delta * indexCos * ONE_DIV_F4)));
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Returns the trigonometric sine of an angle converted as a cos operation.
|
||||
* <p>
|
||||
* Error is around 1E-15.
|
||||
* <p>
|
||||
* Special cases:
|
||||
* <ul>
|
||||
* <li>If the argument is {@code NaN} or an infinity, then the result is {@code NaN}.
|
||||
* </ul>
|
||||
* @param a an angle, in radians.
|
||||
* @return the sine of the argument.
|
||||
* @see Math#cos(double)
|
||||
*/
|
||||
public static double sin(double a) {
|
||||
return cos(a - PIO2);
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the trigonometric tangent of an angle converted in terms of a sin/cos operation.
|
||||
* <p>
|
||||
* Error is around 1E-15.
|
||||
* <p>
|
||||
* Special cases:
|
||||
* <ul>
|
||||
* <li>If the argument is {@code NaN} or an infinity, then the result is {@code NaN}.
|
||||
* </ul>
|
||||
* @param a an angle, in radians.
|
||||
* @return the tangent of the argument.
|
||||
* @see Math#sin(double) aand Math#cos(double)
|
||||
*/
|
||||
public static double tan(double a) {
|
||||
return sin(a) / cos(a);
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the arc sine of a value.
|
||||
* <p>
|
||||
|
|
@ -146,13 +180,15 @@ public class SloppyMath {
|
|||
}
|
||||
|
||||
// haversin
|
||||
private static final double TO_RADIANS = Math.PI / 180D;
|
||||
static final double TO_RADIANS = Math.PI / 180D;
|
||||
static final double TO_DEGREES = 180D / Math.PI;
|
||||
|
||||
// cos/asin
|
||||
private static final double ONE_DIV_F2 = 1/2.0;
|
||||
private static final double ONE_DIV_F3 = 1/6.0;
|
||||
private static final double ONE_DIV_F4 = 1/24.0;
|
||||
|
||||
|
||||
static final double PIO2 = Math.PI / 2D;
|
||||
private static final double PIO2_HI = Double.longBitsToDouble(0x3FF921FB54400000L); // 1.57079632673412561417e+00 first 33 bits of pi/2
|
||||
private static final double PIO2_LO = Double.longBitsToDouble(0x3DD0B4611A626331L); // 6.07710050650619224932e-11 pi/2 - PIO2_HI
|
||||
private static final double TWOPI_HI = 4*PIO2_HI;
|
||||
|
|
|
|||
|
|
@ -27,6 +27,7 @@ import org.apache.lucene.index.DimensionalValues.Relation;
|
|||
import org.apache.lucene.index.LeafReader;
|
||||
import org.apache.lucene.index.LeafReaderContext;
|
||||
import org.apache.lucene.util.DocIdSetBuilder;
|
||||
import org.apache.lucene.util.GeoRelationUtils;
|
||||
import org.apache.lucene.util.GeoUtils;
|
||||
import org.apache.lucene.util.bkd.BKDUtil;
|
||||
|
||||
|
|
@ -125,7 +126,7 @@ public class DimensionalPointInPolygonQuery extends Query {
|
|||
assert packedValue.length == 8;
|
||||
double lat = DimensionalLatLonField.decodeLat(BKDUtil.bytesToInt(packedValue, 0));
|
||||
double lon = DimensionalLatLonField.decodeLon(BKDUtil.bytesToInt(packedValue, 1));
|
||||
if (GeoUtils.pointInPolygon(polyLons, polyLats, lat, lon)) {
|
||||
if (GeoRelationUtils.pointInPolygon(polyLons, polyLats, lat, lon)) {
|
||||
hitCount[0]++;
|
||||
result.add(docID);
|
||||
}
|
||||
|
|
@ -141,11 +142,11 @@ public class DimensionalPointInPolygonQuery extends Query {
|
|||
if (cellMinLat <= minLat && cellMaxLat >= maxLat && cellMinLon <= minLon && cellMaxLon >= maxLon) {
|
||||
// Cell fully encloses the query
|
||||
return Relation.CELL_CROSSES_QUERY;
|
||||
} else if (GeoUtils.rectWithinPoly(cellMinLon, cellMinLat, cellMaxLon, cellMaxLat,
|
||||
} else if (GeoRelationUtils.rectWithinPoly(cellMinLon, cellMinLat, cellMaxLon, cellMaxLat,
|
||||
polyLons, polyLats,
|
||||
minLon, minLat, maxLon, maxLat)) {
|
||||
return Relation.CELL_INSIDE_QUERY;
|
||||
} else if (GeoUtils.rectCrossesPoly(cellMinLon, cellMinLat, cellMaxLon, cellMaxLat,
|
||||
} else if (GeoRelationUtils.rectCrossesPoly(cellMinLon, cellMinLat, cellMaxLon, cellMaxLat,
|
||||
polyLons, polyLats,
|
||||
minLon, minLat, maxLon, maxLat)) {
|
||||
return Relation.CELL_CROSSES_QUERY;
|
||||
|
|
|
|||
|
|
@ -24,6 +24,7 @@ import org.apache.lucene.index.Terms;
|
|||
import org.apache.lucene.index.TermsEnum;
|
||||
import org.apache.lucene.util.AttributeSource;
|
||||
import org.apache.lucene.util.GeoRect;
|
||||
import org.apache.lucene.util.GeoRelationUtils;
|
||||
import org.apache.lucene.util.GeoUtils;
|
||||
import org.apache.lucene.util.SloppyMath;
|
||||
|
||||
|
|
@ -77,12 +78,14 @@ final class GeoPointDistanceQueryImpl extends GeoPointInBBoxQueryImpl {
|
|||
|
||||
@Override
|
||||
protected boolean cellCrosses(final double minLon, final double minLat, final double maxLon, final double maxLat) {
|
||||
return GeoUtils.rectCrossesCircle(minLon, minLat, maxLon, maxLat, centerLon, query.centerLat, query.radiusMeters);
|
||||
return GeoRelationUtils.rectCrossesCircle(minLon, minLat, maxLon, maxLat,
|
||||
centerLon, query.centerLat, query.radiusMeters, true);
|
||||
}
|
||||
|
||||
@Override
|
||||
protected boolean cellWithin(final double minLon, final double minLat, final double maxLon, final double maxLat) {
|
||||
return GeoUtils.rectWithinCircle(minLon, minLat, maxLon, maxLat, centerLon, query.centerLat, query.radiusMeters);
|
||||
return GeoRelationUtils.rectWithinCircle(minLon, minLat, maxLon, maxLat,
|
||||
centerLon, query.centerLat, query.radiusMeters, true);
|
||||
}
|
||||
|
||||
@Override
|
||||
|
|
|
|||
|
|
@ -23,6 +23,7 @@ import org.apache.lucene.document.GeoPointField;
|
|||
import org.apache.lucene.index.Terms;
|
||||
import org.apache.lucene.index.TermsEnum;
|
||||
import org.apache.lucene.util.AttributeSource;
|
||||
import org.apache.lucene.util.GeoRelationUtils;
|
||||
import org.apache.lucene.util.GeoUtils;
|
||||
import org.apache.lucene.util.SloppyMath;
|
||||
import org.apache.lucene.util.ToStringUtils;
|
||||
|
|
@ -83,7 +84,7 @@ class GeoPointInBBoxQueryImpl extends GeoPointTermQuery {
|
|||
*/
|
||||
@Override
|
||||
protected boolean cellCrosses(final double minLon, final double minLat, final double maxLon, final double maxLat) {
|
||||
return GeoUtils.rectCrosses(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
|
||||
return GeoRelationUtils.rectCrosses(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
@ -91,7 +92,7 @@ class GeoPointInBBoxQueryImpl extends GeoPointTermQuery {
|
|||
*/
|
||||
@Override
|
||||
protected boolean cellWithin(final double minLon, final double minLat, final double maxLon, final double maxLat) {
|
||||
return GeoUtils.rectWithin(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
|
||||
return GeoRelationUtils.rectWithin(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
|
||||
}
|
||||
|
||||
@Override
|
||||
|
|
@ -101,7 +102,7 @@ class GeoPointInBBoxQueryImpl extends GeoPointTermQuery {
|
|||
|
||||
@Override
|
||||
protected boolean postFilter(final double lon, final double lat) {
|
||||
return GeoUtils.bboxContains(lon, lat, minLon, minLat, maxLon, maxLat);
|
||||
return GeoRelationUtils.pointInRect(lon, lat, minLon, minLat, maxLon, maxLat);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -24,6 +24,7 @@ import org.apache.lucene.index.Terms;
|
|||
import org.apache.lucene.index.TermsEnum;
|
||||
import org.apache.lucene.util.AttributeSource;
|
||||
import org.apache.lucene.util.GeoRect;
|
||||
import org.apache.lucene.util.GeoRelationUtils;
|
||||
import org.apache.lucene.util.GeoUtils;
|
||||
|
||||
/** Implements a simple point in polygon query on a GeoPoint field. This is based on
|
||||
|
|
@ -147,13 +148,13 @@ public final class GeoPointInPolygonQuery extends GeoPointInBBoxQueryImpl {
|
|||
|
||||
@Override
|
||||
protected boolean cellCrosses(final double minLon, final double minLat, final double maxLon, final double maxLat) {
|
||||
return GeoUtils.rectCrossesPoly(minLon, minLat, maxLon, maxLat, x, y, GeoPointInPolygonQuery.this.minLon,
|
||||
return GeoRelationUtils.rectCrossesPoly(minLon, minLat, maxLon, maxLat, x, y, GeoPointInPolygonQuery.this.minLon,
|
||||
GeoPointInPolygonQuery.this.minLat, GeoPointInPolygonQuery.this.maxLon, GeoPointInPolygonQuery.this.maxLat);
|
||||
}
|
||||
|
||||
@Override
|
||||
protected boolean cellWithin(final double minLon, final double minLat, final double maxLon, final double maxLat) {
|
||||
return GeoUtils.rectWithinPoly(minLon, minLat, maxLon, maxLat, x, y, GeoPointInPolygonQuery.this.minLon,
|
||||
return GeoRelationUtils.rectWithinPoly(minLon, minLat, maxLon, maxLat, x, y, GeoPointInPolygonQuery.this.minLon,
|
||||
GeoPointInPolygonQuery.this.minLat, GeoPointInPolygonQuery.this.maxLon, GeoPointInPolygonQuery.this.maxLat);
|
||||
}
|
||||
|
||||
|
|
@ -168,11 +169,11 @@ public final class GeoPointInPolygonQuery extends GeoPointInBBoxQueryImpl {
|
|||
* {@link org.apache.lucene.search.GeoPointTermsEnum#accept} method is called to match
|
||||
* encoded terms that fall within the bounding box of the polygon. Those documents that pass the initial
|
||||
* bounding box filter are then compared to the provided polygon using the
|
||||
* {@link org.apache.lucene.util.GeoUtils#pointInPolygon} method.
|
||||
* {@link org.apache.lucene.util.GeoRelationUtils#pointInPolygon} method.
|
||||
*/
|
||||
@Override
|
||||
protected boolean postFilter(final double lon, final double lat) {
|
||||
return GeoUtils.pointInPolygon(x, y, lat, lon);
|
||||
return GeoRelationUtils.pointInPolygon(x, y, lat, lon);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -26,6 +26,7 @@ import org.apache.lucene.index.FilteredTermsEnum;
|
|||
import org.apache.lucene.index.TermsEnum;
|
||||
import org.apache.lucene.util.BytesRef;
|
||||
import org.apache.lucene.util.BytesRefBuilder;
|
||||
import org.apache.lucene.util.GeoRelationUtils;
|
||||
import org.apache.lucene.util.GeoUtils;
|
||||
import org.apache.lucene.util.NumericUtils;
|
||||
|
||||
|
|
@ -139,14 +140,14 @@ abstract class GeoPointTermsEnum extends FilteredTermsEnum {
|
|||
* Primary driver for cells intersecting shape boundaries
|
||||
*/
|
||||
protected boolean cellIntersectsMBR(final double minLon, final double minLat, final double maxLon, final double maxLat) {
|
||||
return GeoUtils.rectIntersects(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
|
||||
return GeoRelationUtils.rectIntersects(minLon, minLat, maxLon, maxLat, this.minLon, this.minLat, this.maxLon, this.maxLat);
|
||||
}
|
||||
|
||||
/**
|
||||
* Return whether quad-cell contains the bounding box of this shape
|
||||
*/
|
||||
protected boolean cellContains(final double minLon, final double minLat, final double maxLon, final double maxLat) {
|
||||
return GeoUtils.rectWithin(this.minLon, this.minLat, this.maxLon, this.maxLat, minLon, minLat, maxLon, maxLat);
|
||||
return GeoRelationUtils.rectWithin(this.minLon, this.minLat, this.maxLon, this.maxLat, minLon, minLat, maxLon, maxLat);
|
||||
}
|
||||
|
||||
public boolean boundaryTerm() {
|
||||
|
|
|
|||
|
|
@ -17,6 +17,8 @@ package org.apache.lucene.util;
|
|||
* limitations under the License.
|
||||
*/
|
||||
|
||||
import static org.apache.lucene.util.SloppyMath.TO_RADIANS;
|
||||
|
||||
/**
|
||||
* Reusable geo-spatial distance utility methods.
|
||||
*
|
||||
|
|
@ -24,6 +26,28 @@ package org.apache.lucene.util;
|
|||
*/
|
||||
public class GeoDistanceUtils {
|
||||
|
||||
/**
|
||||
* Compute the great-circle distance using original haversine implementation published by Sinnot in:
|
||||
* R.W. Sinnott, "Virtues of the Haversine", Sky and Telescope, vol. 68, no. 2, 1984, p. 159
|
||||
*
|
||||
* NOTE: this differs from {@link org.apache.lucene.util.SloppyMath#haversin} in that it uses the semi-major axis
|
||||
* of the earth instead of an approximation based on the average latitude of the two points (which can introduce an
|
||||
* additional error up to .337%, or ~67.6 km at the equator)
|
||||
*/
|
||||
public static double haversin(double lat1, double lon1, double lat2, double lon2) {
|
||||
double dLat = TO_RADIANS * (lat2 - lat1);
|
||||
double dLon = TO_RADIANS * (lon2 - lon1);
|
||||
lat1 = TO_RADIANS * (lat1);
|
||||
lat2 = TO_RADIANS * (lat2);
|
||||
|
||||
final double sinDLatO2 = SloppyMath.sin(dLat / 2);
|
||||
final double sinDLonO2 = SloppyMath.sin(dLon / 2);
|
||||
|
||||
double a = sinDLatO2*sinDLatO2 + sinDLonO2 * sinDLonO2 * SloppyMath.cos(lat1) * SloppyMath.cos(lat2);
|
||||
double c = 2 * SloppyMath.asin(Math.sqrt(a));
|
||||
return (GeoProjectionUtils.SEMIMAJOR_AXIS * c);
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute the distance between two geo-points using vincenty distance formula
|
||||
* Vincenty uses the oblate spheroid whereas haversine uses unit sphere, this will give roughly
|
||||
|
|
@ -86,6 +110,20 @@ public class GeoDistanceUtils {
|
|||
return (GeoProjectionUtils.SEMIMINOR_AXIS * A * (sigma - deltaSigma));
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes distance between two points in a cartesian (x, y, {z - optional}) coordinate system
|
||||
*/
|
||||
public static double linearDistance(double[] pt1, double[] pt2) {
|
||||
assert pt1 != null && pt2 != null && pt1.length == pt2.length && pt1.length > 1;
|
||||
final double d0 = pt1[0] - pt2[0];
|
||||
final double d1 = pt1[1] - pt2[1];
|
||||
if (pt1.length == 3) {
|
||||
final double d2 = pt1[2] - pt2[2];
|
||||
return Math.sqrt(d0*d0 + d1*d1 + d2*d2);
|
||||
}
|
||||
return Math.sqrt(d0*d0 + d1*d1);
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute the inverse haversine to determine distance in degrees longitude for provided distance in meters
|
||||
* @param lat latitude to compute delta degrees lon
|
||||
|
|
@ -152,9 +190,9 @@ public class GeoDistanceUtils {
|
|||
/** Returns the maximum distance/radius (in meters) from the point 'center' before overlapping */
|
||||
public static double maxRadialDistanceMeters(final double centerLon, final double centerLat) {
|
||||
if (Math.abs(centerLat) == GeoUtils.MAX_LAT_INCL) {
|
||||
return SloppyMath.haversin(centerLat, centerLon, 0, centerLon)*1000.0;
|
||||
return GeoDistanceUtils.haversin(centerLat, centerLon, 0, centerLon);
|
||||
}
|
||||
return SloppyMath.haversin(centerLat, centerLon, centerLat, (GeoUtils.MAX_LON_INCL + centerLon) % 360)*1000.0;
|
||||
return GeoDistanceUtils.haversin(centerLat, centerLon, centerLat, (GeoUtils.MAX_LON_INCL + centerLon) % 360);
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
|
|||
|
|
@ -17,6 +17,10 @@ package org.apache.lucene.util;
|
|||
* limitations under the License.
|
||||
*/
|
||||
|
||||
import static org.apache.lucene.util.SloppyMath.PIO2;
|
||||
import static org.apache.lucene.util.SloppyMath.TO_DEGREES;
|
||||
import static org.apache.lucene.util.SloppyMath.TO_RADIANS;
|
||||
|
||||
/**
|
||||
* Reusable geo-spatial projection utility methods.
|
||||
*
|
||||
|
|
@ -28,13 +32,14 @@ public class GeoProjectionUtils {
|
|||
public static final double FLATTENING = 1.0/298.257223563;
|
||||
public static final double SEMIMINOR_AXIS = SEMIMAJOR_AXIS * (1.0 - FLATTENING); //6_356_752.31420; // [m]
|
||||
public static final double ECCENTRICITY = StrictMath.sqrt((2.0 - FLATTENING) * FLATTENING);
|
||||
static final double PI_OVER_2 = StrictMath.PI / 2.0D;
|
||||
static final double SEMIMAJOR_AXIS2 = SEMIMAJOR_AXIS * SEMIMAJOR_AXIS;
|
||||
static final double SEMIMINOR_AXIS2 = SEMIMINOR_AXIS * SEMIMINOR_AXIS;
|
||||
public static final double MIN_LON_RADIANS = StrictMath.toRadians(GeoUtils.MIN_LON_INCL);
|
||||
public static final double MIN_LAT_RADIANS = StrictMath.toRadians(GeoUtils.MIN_LAT_INCL);
|
||||
public static final double MAX_LON_RADIANS = StrictMath.toRadians(GeoUtils.MAX_LON_INCL);
|
||||
public static final double MAX_LAT_RADIANS = StrictMath.toRadians(GeoUtils.MAX_LAT_INCL);
|
||||
public static final double MIN_LON_RADIANS = TO_RADIANS * GeoUtils.MIN_LON_INCL;
|
||||
public static final double MIN_LAT_RADIANS = TO_RADIANS * GeoUtils.MIN_LAT_INCL;
|
||||
public static final double MAX_LON_RADIANS = TO_RADIANS * GeoUtils.MAX_LON_INCL;
|
||||
public static final double MAX_LAT_RADIANS = TO_RADIANS * GeoUtils.MAX_LAT_INCL;
|
||||
private static final double E2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMAJOR_AXIS2);
|
||||
private static final double EP2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMINOR_AXIS2);
|
||||
|
||||
/**
|
||||
* Converts from geocentric earth-centered earth-fixed to geodesic lat/lon/alt
|
||||
|
|
@ -47,8 +52,6 @@ public class GeoProjectionUtils {
|
|||
public static final double[] ecfToLLA(final double x, final double y, final double z, double[] lla) {
|
||||
boolean atPole = false;
|
||||
final double ad_c = 1.0026000D;
|
||||
final double e2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMAJOR_AXIS2);
|
||||
final double ep2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMINOR_AXIS2);
|
||||
final double cos67P5 = 0.38268343236508977D;
|
||||
|
||||
if (lla == null) {
|
||||
|
|
@ -59,18 +62,18 @@ public class GeoProjectionUtils {
|
|||
lla[0] = StrictMath.atan2(y,x);
|
||||
} else {
|
||||
if (y > 0) {
|
||||
lla[0] = PI_OVER_2;
|
||||
lla[0] = PIO2;
|
||||
} else if (y < 0) {
|
||||
lla[0] = -PI_OVER_2;
|
||||
lla[0] = -PIO2;
|
||||
} else {
|
||||
atPole = true;
|
||||
lla[0] = 0.0D;
|
||||
if (z > 0.0) {
|
||||
lla[1] = PI_OVER_2;
|
||||
lla[1] = PIO2;
|
||||
} else if (z < 0.0) {
|
||||
lla[1] = -PI_OVER_2;
|
||||
lla[1] = -PIO2;
|
||||
} else {
|
||||
lla[1] = PI_OVER_2;
|
||||
lla[1] = PIO2;
|
||||
lla[2] = -SEMIMINOR_AXIS;
|
||||
return lla;
|
||||
}
|
||||
|
|
@ -84,25 +87,25 @@ public class GeoProjectionUtils {
|
|||
final double sinB0 = t0 / s0;
|
||||
final double cosB0 = w / s0;
|
||||
final double sin3B0 = sinB0 * sinB0 * sinB0;
|
||||
final double t1 = z + SEMIMINOR_AXIS * ep2 * sin3B0;
|
||||
final double sum = w - SEMIMAJOR_AXIS * e2 * cosB0 * cosB0 * cosB0;
|
||||
final double t1 = z + SEMIMINOR_AXIS * EP2 * sin3B0;
|
||||
final double sum = w - SEMIMAJOR_AXIS * E2 * cosB0 * cosB0 * cosB0;
|
||||
final double s1 = StrictMath.sqrt(t1 * t1 + sum * sum);
|
||||
final double sinP1 = t1 / s1;
|
||||
final double cosP1 = sum / s1;
|
||||
final double rn = SEMIMAJOR_AXIS / StrictMath.sqrt(1.0D - e2 * sinP1 * sinP1);
|
||||
final double rn = SEMIMAJOR_AXIS / StrictMath.sqrt(1.0D - E2 * sinP1 * sinP1);
|
||||
|
||||
if (cosP1 >= cos67P5) {
|
||||
lla[2] = w / cosP1 - rn;
|
||||
} else if (cosP1 <= -cos67P5) {
|
||||
lla[2] = w / -cosP1 - rn;
|
||||
} else {
|
||||
lla[2] = z / sinP1 + rn * (e2 - 1.0);
|
||||
lla[2] = z / sinP1 + rn * (E2 - 1.0);
|
||||
}
|
||||
if (!atPole) {
|
||||
lla[1] = StrictMath.atan(sinP1/cosP1);
|
||||
}
|
||||
lla[0] = StrictMath.toDegrees(lla[0]);
|
||||
lla[1] = StrictMath.toDegrees(lla[1]);
|
||||
lla[0] = TO_DEGREES * lla[0];
|
||||
lla[1] = TO_DEGREES * lla[1];
|
||||
|
||||
return lla;
|
||||
}
|
||||
|
|
@ -116,33 +119,32 @@ public class GeoProjectionUtils {
|
|||
* @return either a new ecef array or the reusable ecf parameter
|
||||
*/
|
||||
public static final double[] llaToECF(double lon, double lat, double alt, double[] ecf) {
|
||||
lon = StrictMath.toRadians(lon);
|
||||
lat = StrictMath.toRadians(lat);
|
||||
lon = TO_RADIANS * lon;
|
||||
lat = TO_RADIANS * lat;
|
||||
|
||||
final double sl = StrictMath.sin(lat);
|
||||
final double sl = SloppyMath.sin(lat);
|
||||
final double s2 = sl*sl;
|
||||
final double cl = StrictMath.cos(lat);
|
||||
final double ge2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMAJOR_AXIS2);
|
||||
final double cl = SloppyMath.cos(lat);
|
||||
|
||||
if (ecf == null) {
|
||||
ecf = new double[3];
|
||||
}
|
||||
|
||||
if (lat < -PI_OVER_2 && lat > -1.001D * PI_OVER_2) {
|
||||
lat = -PI_OVER_2;
|
||||
} else if (lat > PI_OVER_2 && lat < 1.001D * PI_OVER_2) {
|
||||
lat = PI_OVER_2;
|
||||
if (lat < -PIO2 && lat > -1.001D * PIO2) {
|
||||
lat = -PIO2;
|
||||
} else if (lat > PIO2 && lat < 1.001D * PIO2) {
|
||||
lat = PIO2;
|
||||
}
|
||||
assert (lat >= -PI_OVER_2) || (lat <= PI_OVER_2);
|
||||
assert (lat >= -PIO2) || (lat <= PIO2);
|
||||
|
||||
if (lon > StrictMath.PI) {
|
||||
lon -= (2*StrictMath.PI);
|
||||
}
|
||||
|
||||
final double rn = SEMIMAJOR_AXIS / StrictMath.sqrt(1.0D - ge2 * s2);
|
||||
ecf[0] = (rn+alt) * cl * StrictMath.cos(lon);
|
||||
ecf[1] = (rn+alt) * cl * StrictMath.sin(lon);
|
||||
ecf[2] = ((rn*(1.0-ge2))+alt)*sl;
|
||||
final double rn = SEMIMAJOR_AXIS / StrictMath.sqrt(1.0D - E2 * s2);
|
||||
ecf[0] = (rn+alt) * cl * SloppyMath.cos(lon);
|
||||
ecf[1] = (rn+alt) * cl * SloppyMath.sin(lon);
|
||||
ecf[2] = ((rn*(1.0-E2))+alt)*sl;
|
||||
|
||||
return ecf;
|
||||
}
|
||||
|
|
@ -272,13 +274,13 @@ public class GeoProjectionUtils {
|
|||
phiMatrix = new double[3][3];
|
||||
}
|
||||
|
||||
originLon = StrictMath.toRadians(originLon);
|
||||
originLat = StrictMath.toRadians(originLat);
|
||||
originLon = TO_RADIANS * originLon;
|
||||
originLat = TO_RADIANS * originLat;
|
||||
|
||||
final double sLon = StrictMath.sin(originLon);
|
||||
final double cLon = StrictMath.cos(originLon);
|
||||
final double sLat = StrictMath.sin(originLat);
|
||||
final double cLat = StrictMath.cos(originLat);
|
||||
final double sLon = SloppyMath.sin(originLon);
|
||||
final double cLon = SloppyMath.cos(originLon);
|
||||
final double sLat = SloppyMath.sin(originLat);
|
||||
final double cLat = SloppyMath.cos(originLat);
|
||||
|
||||
phiMatrix[0][0] = -sLon;
|
||||
phiMatrix[0][1] = cLon;
|
||||
|
|
@ -306,13 +308,13 @@ public class GeoProjectionUtils {
|
|||
phiMatrix = new double[3][3];
|
||||
}
|
||||
|
||||
originLon = StrictMath.toRadians(originLon);
|
||||
originLat = StrictMath.toRadians(originLat);
|
||||
originLon = TO_RADIANS * originLon;
|
||||
originLat = TO_RADIANS * originLat;
|
||||
|
||||
final double sLat = StrictMath.sin(originLat);
|
||||
final double cLat = StrictMath.cos(originLat);
|
||||
final double sLon = StrictMath.sin(originLon);
|
||||
final double cLon = StrictMath.cos(originLon);
|
||||
final double sLat = SloppyMath.sin(originLat);
|
||||
final double cLat = SloppyMath.cos(originLat);
|
||||
final double sLon = SloppyMath.sin(originLon);
|
||||
final double cLon = SloppyMath.cos(originLon);
|
||||
|
||||
phiMatrix[0][0] = -sLon;
|
||||
phiMatrix[1][0] = cLon;
|
||||
|
|
@ -337,22 +339,22 @@ public class GeoProjectionUtils {
|
|||
* @param pt resulting point
|
||||
* @return the point along a bearing at a given distance in meters
|
||||
*/
|
||||
public static final double[] pointFromLonLatBearing(double lon, double lat, double bearing, double dist, double[] pt) {
|
||||
public static final double[] pointFromLonLatBearingVincenty(double lon, double lat, double bearing, double dist, double[] pt) {
|
||||
|
||||
if (pt == null) {
|
||||
pt = new double[2];
|
||||
}
|
||||
|
||||
final double alpha1 = StrictMath.toRadians(bearing);
|
||||
final double cosA1 = StrictMath.cos(alpha1);
|
||||
final double sinA1 = StrictMath.sin(alpha1);
|
||||
final double tanU1 = (1-FLATTENING) * StrictMath.tan(StrictMath.toRadians(lat));
|
||||
final double alpha1 = TO_RADIANS * bearing;
|
||||
final double cosA1 = SloppyMath.cos(alpha1);
|
||||
final double sinA1 = SloppyMath.sin(alpha1);
|
||||
final double tanU1 = (1-FLATTENING) * SloppyMath.tan(TO_RADIANS * lat);
|
||||
final double cosU1 = 1 / StrictMath.sqrt((1+tanU1*tanU1));
|
||||
final double sinU1 = tanU1*cosU1;
|
||||
final double sig1 = StrictMath.atan2(tanU1, cosA1);
|
||||
final double sinAlpha = cosU1 * sinA1;
|
||||
final double cosSqAlpha = 1 - sinAlpha*sinAlpha;
|
||||
final double uSq = cosSqAlpha * (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2) / SEMIMINOR_AXIS2;
|
||||
final double uSq = cosSqAlpha * EP2;
|
||||
final double A = 1 + uSq/16384D*(4096D + uSq * (-768D + uSq * (320D - 175D*uSq)));
|
||||
final double B = uSq/1024D * (256D + uSq * (-128D + uSq * (74D - 47D * uSq)));
|
||||
|
||||
|
|
@ -361,9 +363,9 @@ public class GeoProjectionUtils {
|
|||
double sinSigma, cosSigma, cos2SigmaM, deltaSigma;
|
||||
|
||||
do {
|
||||
cos2SigmaM = StrictMath.cos(2*sig1 + sigma);
|
||||
sinSigma = StrictMath.sin(sigma);
|
||||
cosSigma = StrictMath.cos(sigma);
|
||||
cos2SigmaM = SloppyMath.cos(2*sig1 + sigma);
|
||||
sinSigma = SloppyMath.sin(sigma);
|
||||
cosSigma = SloppyMath.cos(sigma);
|
||||
|
||||
deltaSigma = B * sinSigma * (cos2SigmaM + (B/4D) * (cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
|
||||
(B/6) * cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
|
||||
|
|
@ -379,9 +381,42 @@ public class GeoProjectionUtils {
|
|||
|
||||
final double lam = lambda - (1-c) * FLATTENING * sinAlpha *
|
||||
(sigma + c * sinSigma * (cos2SigmaM + c * cosSigma * (-1 + 2* cos2SigmaM*cos2SigmaM)));
|
||||
pt[0] = GeoUtils.normalizeLon(lon + StrictMath.toDegrees(lam));
|
||||
pt[1] = GeoUtils.normalizeLat(StrictMath.toDegrees(lat2));
|
||||
pt[0] = GeoUtils.normalizeLon(lon + TO_DEGREES * lam);
|
||||
pt[1] = GeoUtils.normalizeLat(TO_DEGREES * lat2);
|
||||
|
||||
return pt;
|
||||
}
|
||||
|
||||
/**
|
||||
* Finds a point along a bearing from a given lon,lat geolocation using great circle arc
|
||||
*
|
||||
* @param lon origin longitude in degrees
|
||||
* @param lat origin latitude in degrees
|
||||
* @param bearing azimuthal bearing in degrees
|
||||
* @param dist distance in meters
|
||||
* @param pt resulting point
|
||||
* @return the point along a bearing at a given distance in meters
|
||||
*/
|
||||
public static final double[] pointFromLonLatBearingGreatCircle(double lon, double lat, double bearing, double dist, double[] pt) {
|
||||
|
||||
if (pt == null) {
|
||||
pt = new double[2];
|
||||
}
|
||||
|
||||
lon *= TO_RADIANS;
|
||||
lat *= TO_RADIANS;
|
||||
bearing *= TO_RADIANS;
|
||||
|
||||
final double cLat = SloppyMath.cos(lat);
|
||||
final double sLat = SloppyMath.sin(lat);
|
||||
final double sinDoR = SloppyMath.sin(dist / GeoProjectionUtils.SEMIMAJOR_AXIS);
|
||||
final double cosDoR = SloppyMath.cos(dist / GeoProjectionUtils.SEMIMAJOR_AXIS);
|
||||
|
||||
pt[1] = SloppyMath.asin(sLat*cosDoR + cLat * sinDoR * SloppyMath.cos(bearing));
|
||||
pt[0] = TO_DEGREES * (lon + Math.atan2(SloppyMath.sin(bearing) * sinDoR * cLat, cosDoR - sLat * SloppyMath.sin(pt[1])));
|
||||
pt[1] *= TO_DEGREES;
|
||||
|
||||
return pt;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -0,0 +1,356 @@
|
|||
package org.apache.lucene.util;
|
||||
|
||||
/*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
/**
|
||||
* Reusable geo-relation utility methods
|
||||
*/
|
||||
public class GeoRelationUtils {
|
||||
|
||||
/**
|
||||
* Determine if a bbox (defined by minLon, minLat, maxLon, maxLat) contains the provided point (defined by lon, lat)
|
||||
* NOTE: this is a basic method that does not handle dateline or pole crossing. Unwrapping must be done before
|
||||
* calling this method.
|
||||
*/
|
||||
public static boolean pointInRect(final double lon, final double lat, final double minLon,
|
||||
final double minLat, final double maxLon, final double maxLat) {
|
||||
return (GeoUtils.compare(lon, minLon) >= 0 && GeoUtils.compare(lon, maxLon) <= 0
|
||||
&& GeoUtils.compare(lat, minLat) >= 0 && GeoUtils.compare(lat, maxLat) <= 0);
|
||||
}
|
||||
|
||||
/**
|
||||
* simple even-odd point in polygon computation
|
||||
* 1. Determine if point is contained in the longitudinal range
|
||||
* 2. Determine whether point crosses the edge by computing the latitudinal delta
|
||||
* between the end-point of a parallel vector (originating at the point) and the
|
||||
* y-component of the edge sink
|
||||
*
|
||||
* NOTE: Requires polygon point (x,y) order either clockwise or counter-clockwise
|
||||
*/
|
||||
public static boolean pointInPolygon(double[] x, double[] y, double lat, double lon) {
|
||||
assert x.length == y.length;
|
||||
boolean inPoly = false;
|
||||
/**
|
||||
* Note: This is using a euclidean coordinate system which could result in
|
||||
* upwards of 110KM error at the equator.
|
||||
* TODO convert coordinates to cylindrical projection (e.g. mercator)
|
||||
*/
|
||||
for (int i = 1; i < x.length; i++) {
|
||||
if (x[i] < lon && x[i-1] >= lon || x[i-1] < lon && x[i] >= lon) {
|
||||
if (y[i] + (lon - x[i]) / (x[i-1] - x[i]) * (y[i-1] - y[i]) < lat) {
|
||||
inPoly = !inPoly;
|
||||
}
|
||||
}
|
||||
}
|
||||
return inPoly;
|
||||
}
|
||||
|
||||
public static boolean rectDisjoint(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
|
||||
final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
|
||||
return (aMaxX < bMinX || aMinX > bMaxX || aMaxY < bMinY || aMinY > bMaxY);
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether the first (a) rectangle is wholly within another (b) rectangle (shared boundaries allowed)
|
||||
*/
|
||||
public static boolean rectWithin(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
|
||||
final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
|
||||
return !(aMinX < bMinX || aMinY < bMinY || aMaxX > bMaxX || aMaxY > bMaxY);
|
||||
}
|
||||
|
||||
public static boolean rectCrosses(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
|
||||
final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
|
||||
return !(rectDisjoint(aMinX, aMinY, aMaxX, aMaxY, bMinX, bMinY, bMaxX, bMaxY) ||
|
||||
rectWithin(aMinX, aMinY, aMaxX, aMaxY, bMinX, bMinY, bMaxX, bMaxY));
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether rectangle a contains rectangle b (touching allowed)
|
||||
*/
|
||||
public static boolean rectContains(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
|
||||
final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
|
||||
return !(bMinX < aMinX || bMinY < aMinY || bMaxX > aMaxX || bMaxY > aMaxY);
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether a rectangle intersects another rectangle (crosses, within, touching, etc)
|
||||
*/
|
||||
public static boolean rectIntersects(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
|
||||
final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
|
||||
return !((aMaxX < bMinX || aMinX > bMaxX || aMaxY < bMinY || aMinY > bMaxY) );
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether a rectangle crosses a shape. (touching not allowed)
|
||||
*/
|
||||
public static boolean rectCrossesPoly(final double rMinX, final double rMinY, final double rMaxX,
|
||||
final double rMaxY, final double[] shapeX, final double[] shapeY,
|
||||
final double sMinX, final double sMinY, final double sMaxX,
|
||||
final double sMaxY) {
|
||||
// short-circuit: if the bounding boxes are disjoint then the shape does not cross
|
||||
if (rectDisjoint(rMinX, rMinY, rMaxX, rMaxY, sMinX, sMinY, sMaxX, sMaxY)) {
|
||||
return false;
|
||||
}
|
||||
|
||||
final double[][] bbox = new double[][] { {rMinX, rMinY}, {rMaxX, rMinY}, {rMaxX, rMaxY}, {rMinX, rMaxY}, {rMinX, rMinY} };
|
||||
final int polyLength = shapeX.length-1;
|
||||
double d, s, t, a1, b1, c1, a2, b2, c2;
|
||||
double x00, y00, x01, y01, x10, y10, x11, y11;
|
||||
|
||||
// computes the intersection point between each bbox edge and the polygon edge
|
||||
for (short b=0; b<4; ++b) {
|
||||
a1 = bbox[b+1][1]-bbox[b][1];
|
||||
b1 = bbox[b][0]-bbox[b+1][0];
|
||||
c1 = a1*bbox[b+1][0] + b1*bbox[b+1][1];
|
||||
for (int p=0; p<polyLength; ++p) {
|
||||
a2 = shapeY[p+1]-shapeY[p];
|
||||
b2 = shapeX[p]-shapeX[p+1];
|
||||
// compute determinant
|
||||
d = a1*b2 - a2*b1;
|
||||
if (d != 0) {
|
||||
// lines are not parallel, check intersecting points
|
||||
c2 = a2*shapeX[p+1] + b2*shapeY[p+1];
|
||||
s = (1/d)*(b2*c1 - b1*c2);
|
||||
t = (1/d)*(a1*c2 - a2*c1);
|
||||
x00 = StrictMath.min(bbox[b][0], bbox[b+1][0]) - GeoUtils.TOLERANCE;
|
||||
x01 = StrictMath.max(bbox[b][0], bbox[b+1][0]) + GeoUtils.TOLERANCE;
|
||||
y00 = StrictMath.min(bbox[b][1], bbox[b+1][1]) - GeoUtils.TOLERANCE;
|
||||
y01 = StrictMath.max(bbox[b][1], bbox[b+1][1]) + GeoUtils.TOLERANCE;
|
||||
x10 = StrictMath.min(shapeX[p], shapeX[p+1]) - GeoUtils.TOLERANCE;
|
||||
x11 = StrictMath.max(shapeX[p], shapeX[p+1]) + GeoUtils.TOLERANCE;
|
||||
y10 = StrictMath.min(shapeY[p], shapeY[p+1]) - GeoUtils.TOLERANCE;
|
||||
y11 = StrictMath.max(shapeY[p], shapeY[p+1]) + GeoUtils.TOLERANCE;
|
||||
// check whether the intersection point is touching one of the line segments
|
||||
boolean touching = ((x00 == s && y00 == t) || (x01 == s && y01 == t))
|
||||
|| ((x10 == s && y10 == t) || (x11 == s && y11 == t));
|
||||
// if line segments are not touching and the intersection point is within the range of either segment
|
||||
if (!(touching || x00 > s || x01 < s || y00 > t || y01 < t || x10 > s || x11 < s || y10 > t || y11 < t)) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
} // for each poly edge
|
||||
} // for each bbox edge
|
||||
return false;
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether a rectangle is within a given polygon (shared boundaries allowed)
|
||||
*/
|
||||
public static boolean rectWithinPoly(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double[] shapeX, final double[] shapeY, final double sMinX,
|
||||
final double sMinY, final double sMaxX, final double sMaxY) {
|
||||
// check if rectangle crosses poly (to handle concave/pacman polys), then check that all 4 corners
|
||||
// are contained
|
||||
return !(rectCrossesPoly(rMinX, rMinY, rMaxX, rMaxY, shapeX, shapeY, sMinX, sMinY, sMaxX, sMaxY) ||
|
||||
!pointInPolygon(shapeX, shapeY, rMinY, rMinX) || !pointInPolygon(shapeX, shapeY, rMinY, rMaxX) ||
|
||||
!pointInPolygon(shapeX, shapeY, rMaxY, rMaxX) || !pointInPolygon(shapeX, shapeY, rMaxY, rMinX));
|
||||
}
|
||||
|
||||
private static boolean rectAnyCornersInCircle(final double rMinX, final double rMinY, final double rMaxX,
|
||||
final double rMaxY, final double centerLon, final double centerLat,
|
||||
final double radiusMeters, final boolean approx) {
|
||||
if (approx == true) {
|
||||
return rectAnyCornersInCircleSloppy(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters);
|
||||
}
|
||||
double w = Math.abs(rMaxX - rMinX);
|
||||
if (w <= 90.0) {
|
||||
return GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMinX) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMinX) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMaxX) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMaxX) <= radiusMeters;
|
||||
}
|
||||
// partition
|
||||
w /= 4;
|
||||
final double p1 = rMinX + w;
|
||||
final double p2 = p1 + w;
|
||||
final double p3 = p2 + w;
|
||||
|
||||
return GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMinX) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMinX) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p1) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p1) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p2) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p2) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p3) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p3) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMaxX) <= radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMaxX) <= radiusMeters;
|
||||
}
|
||||
|
||||
private static boolean rectAnyCornersInCircleSloppy(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters) {
|
||||
return SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 <= radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 <= radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 <= radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 <= radiusMeters;
|
||||
}
|
||||
|
||||
private static boolean rectAnyCornersOutsideCircle(final double rMinX, final double rMinY, final double rMaxX,
|
||||
final double rMaxY, final double centerLon, final double centerLat,
|
||||
final double radiusMeters, final boolean approx) {
|
||||
if (approx == true) {
|
||||
return rectAnyCornersOutsideCircleSloppy(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters);
|
||||
}
|
||||
double w = Math.abs(rMaxX - rMinX);
|
||||
if (w <= 90.0) {
|
||||
return GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMinX) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMinX) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMaxX) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMaxX) > radiusMeters;
|
||||
}
|
||||
// partition
|
||||
w /= 4;
|
||||
final double p1 = rMinX + w;
|
||||
final double p2 = p1 + w;
|
||||
final double p3 = p2 + w;
|
||||
|
||||
return GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMinX) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMinX) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p1) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p1) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p2) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p2) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, p3) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, p3) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMaxY, rMaxX) > radiusMeters
|
||||
|| GeoDistanceUtils.haversin(centerLat, centerLon, rMinY, rMaxX) > radiusMeters;
|
||||
}
|
||||
|
||||
private static boolean rectAnyCornersOutsideCircleSloppy(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters) {
|
||||
return SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 > radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 > radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 > radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 > radiusMeters;
|
||||
}
|
||||
|
||||
public static boolean rectWithinCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters) {
|
||||
return rectWithinCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, false);
|
||||
}
|
||||
|
||||
public static boolean rectWithinCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters,
|
||||
final boolean approx) {
|
||||
return rectAnyCornersOutsideCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx) == false;
|
||||
}
|
||||
|
||||
/**
|
||||
* Determine if a bbox (defined by minLon, minLat, maxLon, maxLat) contains the provided point (defined by lon, lat)
|
||||
* NOTE: this is basic method that does not handle dateline or pole crossing. Unwrapping must be done before
|
||||
* calling this method.
|
||||
*/
|
||||
public static boolean rectCrossesCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters) {
|
||||
return rectCrossesCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, false);
|
||||
}
|
||||
|
||||
public static boolean rectCrossesCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters,
|
||||
final boolean approx) {
|
||||
if (approx == true) {
|
||||
return rectAnyCornersInCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx)
|
||||
|| isClosestPointOnRectWithinRange(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx);
|
||||
}
|
||||
|
||||
return (rectAnyCornersInCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx) &&
|
||||
rectAnyCornersOutsideCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx))
|
||||
|| isClosestPointOnRectWithinRange(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters, approx);
|
||||
}
|
||||
|
||||
private static boolean isClosestPointOnRectWithinRange(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters,
|
||||
final boolean approx) {
|
||||
double[] closestPt = {0, 0};
|
||||
GeoDistanceUtils.closestPointOnBBox(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, closestPt);
|
||||
boolean haverShortCut = GeoDistanceUtils.haversin(centerLat, centerLon, closestPt[1], closestPt[0]) <= radiusMeters;
|
||||
if (approx == true || haverShortCut == true) {
|
||||
return haverShortCut;
|
||||
}
|
||||
double lon1 = rMinX;
|
||||
double lon2 = rMaxX;
|
||||
double lat1 = rMinY;
|
||||
double lat2 = rMaxY;
|
||||
if (closestPt[0] == rMinX || closestPt[0] == rMaxX) {
|
||||
lon1 = closestPt[0];
|
||||
lon2 = lon1;
|
||||
} else if (closestPt[1] == rMinY || closestPt[1] == rMaxY) {
|
||||
lat1 = closestPt[1];
|
||||
lat2 = lat1;
|
||||
}
|
||||
|
||||
return lineCrossesSphere(lon1, lat1, 0, lon2, lat2, 0, centerLon, centerLat, 0, radiusMeters);
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether or a 3dimensional line segment intersects or crosses a sphere
|
||||
*
|
||||
* @param lon1 longitudinal location of the line segment start point (in degrees)
|
||||
* @param lat1 latitudinal location of the line segment start point (in degrees)
|
||||
* @param alt1 altitude of the line segment start point (in degrees)
|
||||
* @param lon2 longitudinal location of the line segment end point (in degrees)
|
||||
* @param lat2 latitudinal location of the line segment end point (in degrees)
|
||||
* @param alt2 altitude of the line segment end point (in degrees)
|
||||
* @param centerLon longitudinal location of center search point (in degrees)
|
||||
* @param centerLat latitudinal location of center search point (in degrees)
|
||||
* @param centerAlt altitude of the center point (in meters)
|
||||
* @param radiusMeters search sphere radius (in meters)
|
||||
* @return whether the provided line segment is a secant of the
|
||||
*/
|
||||
private static boolean lineCrossesSphere(double lon1, double lat1, double alt1, double lon2,
|
||||
double lat2, double alt2, double centerLon, double centerLat,
|
||||
double centerAlt, double radiusMeters) {
|
||||
// convert to cartesian 3d (in meters)
|
||||
double[] ecf1 = GeoProjectionUtils.llaToECF(lon1, lat1, alt1, null);
|
||||
double[] ecf2 = GeoProjectionUtils.llaToECF(lon2, lat2, alt2, null);
|
||||
double[] cntr = GeoProjectionUtils.llaToECF(centerLon, centerLat, centerAlt, null);
|
||||
|
||||
// convert radius from arc radius to cartesian radius
|
||||
double[] oneEighty = GeoProjectionUtils.pointFromLonLatBearingGreatCircle(centerLon, centerLat, 180.0d, radiusMeters, new double[3]);
|
||||
GeoProjectionUtils.llaToECF(oneEighty[0], oneEighty[1], 0, oneEighty);
|
||||
|
||||
radiusMeters = GeoDistanceUtils.linearDistance(oneEighty, cntr);// Math.sqrt(oneEighty[0]*cntr[0] + oneEighty[1]*cntr[1] + oneEighty[2]*cntr[2]);
|
||||
|
||||
final double dX = ecf2[0] - ecf1[0];
|
||||
final double dY = ecf2[1] - ecf1[1];
|
||||
final double dZ = ecf2[2] - ecf1[2];
|
||||
final double fX = ecf1[0] - cntr[0];
|
||||
final double fY = ecf1[1] - cntr[1];
|
||||
final double fZ = ecf1[2] - cntr[2];
|
||||
|
||||
final double a = dX*dX + dY*dY + dZ*dZ;
|
||||
final double b = 2 * (fX*dX + fY*dY + fZ*dZ);
|
||||
final double c = (fX*fX + fY*fY + fZ*fZ) - (radiusMeters*radiusMeters);
|
||||
|
||||
double discrim = (b*b)-(4*a*c);
|
||||
if (discrim < 0) {
|
||||
return false;
|
||||
}
|
||||
|
||||
discrim = StrictMath.sqrt(discrim);
|
||||
final double a2 = 2*a;
|
||||
final double t1 = (-b - discrim)/a2;
|
||||
final double t2 = (-b + discrim)/a2;
|
||||
|
||||
if ( (t1 < 0 || t1 > 1) ) {
|
||||
return !(t2 < 0 || t2 > 1);
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
|
@ -19,6 +19,9 @@ package org.apache.lucene.util;
|
|||
|
||||
import java.util.ArrayList;
|
||||
|
||||
import static org.apache.lucene.util.SloppyMath.TO_DEGREES;
|
||||
import static org.apache.lucene.util.SloppyMath.TO_RADIANS;
|
||||
|
||||
/**
|
||||
* Basic reusable geo-spatial utility methods
|
||||
*
|
||||
|
|
@ -74,6 +77,9 @@ public final class GeoUtils {
|
|||
return (val / LAT_SCALE) + MIN_LAT_INCL;
|
||||
}
|
||||
|
||||
/**
|
||||
* Compare two position values within a {@link org.apache.lucene.util.GeoUtils#TOLERANCE} factor
|
||||
*/
|
||||
public static double compare(final double v1, final double v2) {
|
||||
final double delta = v1-v2;
|
||||
return Math.abs(delta) <= TOLERANCE ? 0 : delta;
|
||||
|
|
@ -107,44 +113,6 @@ public final class GeoUtils {
|
|||
return (off <= 180 ? off : 360-off) - 90;
|
||||
}
|
||||
|
||||
/**
|
||||
* Determine if a bbox (defined by minLon, minLat, maxLon, maxLat) contains the provided point (defined by lon, lat)
|
||||
* NOTE: this is a basic method that does not handle dateline or pole crossing. Unwrapping must be done before
|
||||
* calling this method.
|
||||
*/
|
||||
public static boolean bboxContains(final double lon, final double lat, final double minLon,
|
||||
final double minLat, final double maxLon, final double maxLat) {
|
||||
return (compare(lon, minLon) >= 0 && compare(lon, maxLon) <= 0
|
||||
&& compare(lat, minLat) >= 0 && compare(lat, maxLat) <= 0);
|
||||
}
|
||||
|
||||
/**
|
||||
* simple even-odd point in polygon computation
|
||||
* 1. Determine if point is contained in the longitudinal range
|
||||
* 2. Determine whether point crosses the edge by computing the latitudinal delta
|
||||
* between the end-point of a parallel vector (originating at the point) and the
|
||||
* y-component of the edge sink
|
||||
*
|
||||
* NOTE: Requires polygon point (x,y) order either clockwise or counter-clockwise
|
||||
*/
|
||||
public static boolean pointInPolygon(double[] x, double[] y, double lat, double lon) {
|
||||
assert x.length == y.length;
|
||||
boolean inPoly = false;
|
||||
/**
|
||||
* Note: This is using a euclidean coordinate system which could result in
|
||||
* upwards of 110KM error at the equator.
|
||||
* TODO convert coordinates to cylindrical projection (e.g. mercator)
|
||||
*/
|
||||
for (int i = 1; i < x.length; i++) {
|
||||
if (x[i] < lon && x[i-1] >= lon || x[i-1] < lon && x[i] >= lon) {
|
||||
if (y[i] + (lon - x[i]) / (x[i-1] - x[i]) * (y[i-1] - y[i]) < lat) {
|
||||
inPoly = !inPoly;
|
||||
}
|
||||
}
|
||||
}
|
||||
return inPoly;
|
||||
}
|
||||
|
||||
public static String geoTermToString(long term) {
|
||||
StringBuilder s = new StringBuilder(64);
|
||||
final int numberOfLeadingZeros = Long.numberOfLeadingZeros(term);
|
||||
|
|
@ -157,95 +125,6 @@ public final class GeoUtils {
|
|||
return s.toString();
|
||||
}
|
||||
|
||||
|
||||
public static boolean rectDisjoint(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
|
||||
final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
|
||||
return (aMaxX < bMinX || aMinX > bMaxX || aMaxY < bMinY || aMinY > bMaxY);
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether the first (a) rectangle is wholly within another (b) rectangle (shared boundaries allowed)
|
||||
*/
|
||||
public static boolean rectWithin(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
|
||||
final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
|
||||
return !(aMinX < bMinX || aMinY < bMinY || aMaxX > bMaxX || aMaxY > bMaxY);
|
||||
}
|
||||
|
||||
public static boolean rectCrosses(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
|
||||
final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
|
||||
return !(rectDisjoint(aMinX, aMinY, aMaxX, aMaxY, bMinX, bMinY, bMaxX, bMaxY) ||
|
||||
rectWithin(aMinX, aMinY, aMaxX, aMaxY, bMinX, bMinY, bMaxX, bMaxY));
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether rectangle a contains rectangle b (touching allowed)
|
||||
*/
|
||||
public static boolean rectContains(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
|
||||
final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
|
||||
return !(bMinX < aMinX || bMinY < aMinY || bMaxX > aMaxX || bMaxY > aMaxY);
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether a rectangle intersects another rectangle (crosses, within, touching, etc)
|
||||
*/
|
||||
public static boolean rectIntersects(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
|
||||
final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
|
||||
return !((aMaxX < bMinX || aMinX > bMaxX || aMaxY < bMinY || aMinY > bMaxY) );
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether a rectangle crosses a shape. (touching not allowed)
|
||||
*/
|
||||
public static boolean rectCrossesPoly(final double rMinX, final double rMinY, final double rMaxX,
|
||||
final double rMaxY, final double[] shapeX, final double[] shapeY,
|
||||
final double sMinX, final double sMinY, final double sMaxX,
|
||||
final double sMaxY) {
|
||||
// short-circuit: if the bounding boxes are disjoint then the shape does not cross
|
||||
if (rectDisjoint(rMinX, rMinY, rMaxX, rMaxY, sMinX, sMinY, sMaxX, sMaxY)) {
|
||||
return false;
|
||||
}
|
||||
|
||||
final double[][] bbox = new double[][] { {rMinX, rMinY}, {rMaxX, rMinY}, {rMaxX, rMaxY}, {rMinX, rMaxY}, {rMinX, rMinY} };
|
||||
final int polyLength = shapeX.length-1;
|
||||
double d, s, t, a1, b1, c1, a2, b2, c2;
|
||||
double x00, y00, x01, y01, x10, y10, x11, y11;
|
||||
|
||||
// computes the intersection point between each bbox edge and the polygon edge
|
||||
for (short b=0; b<4; ++b) {
|
||||
a1 = bbox[b+1][1]-bbox[b][1];
|
||||
b1 = bbox[b][0]-bbox[b+1][0];
|
||||
c1 = a1*bbox[b+1][0] + b1*bbox[b+1][1];
|
||||
for (int p=0; p<polyLength; ++p) {
|
||||
a2 = shapeY[p+1]-shapeY[p];
|
||||
b2 = shapeX[p]-shapeX[p+1];
|
||||
// compute determinant
|
||||
d = a1*b2 - a2*b1;
|
||||
if (d != 0) {
|
||||
// lines are not parallel, check intersecting points
|
||||
c2 = a2*shapeX[p+1] + b2*shapeY[p+1];
|
||||
s = (1/d)*(b2*c1 - b1*c2);
|
||||
t = (1/d)*(a1*c2 - a2*c1);
|
||||
x00 = StrictMath.min(bbox[b][0], bbox[b+1][0]) - TOLERANCE;
|
||||
x01 = StrictMath.max(bbox[b][0], bbox[b+1][0]) + TOLERANCE;
|
||||
y00 = StrictMath.min(bbox[b][1], bbox[b+1][1]) - TOLERANCE;
|
||||
y01 = StrictMath.max(bbox[b][1], bbox[b+1][1]) + TOLERANCE;
|
||||
x10 = StrictMath.min(shapeX[p], shapeX[p+1]) - TOLERANCE;
|
||||
x11 = StrictMath.max(shapeX[p], shapeX[p+1]) + TOLERANCE;
|
||||
y10 = StrictMath.min(shapeY[p], shapeY[p+1]) - TOLERANCE;
|
||||
y11 = StrictMath.max(shapeY[p], shapeY[p+1]) + TOLERANCE;
|
||||
// check whether the intersection point is touching one of the line segments
|
||||
boolean touching = ((x00 == s && y00 == t) || (x01 == s && y01 == t))
|
||||
|| ((x10 == s && y10 == t) || (x11 == s && y11 == t));
|
||||
// if line segments are not touching and the intersection point is within the range of either segment
|
||||
if (!(touching || x00 > s || x01 < s || y00 > t || y01 < t || x10 > s || x11 < s || y10 > t || y11 < t)) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
} // for each poly edge
|
||||
} // for each bbox edge
|
||||
return false;
|
||||
}
|
||||
|
||||
/**
|
||||
* Converts a given circle (defined as a point/radius) to an approximated line-segment polygon
|
||||
*
|
||||
|
|
@ -267,7 +146,7 @@ public final class GeoUtils {
|
|||
final int sidesLen = sides-1;
|
||||
for (int i=0; i<sidesLen; ++i) {
|
||||
angle = (i*360/sides);
|
||||
pt = GeoProjectionUtils.pointFromLonLatBearing(lon, lat, angle, radiusMeters, pt);
|
||||
pt = GeoProjectionUtils.pointFromLonLatBearingGreatCircle(lon, lat, angle, radiusMeters, pt);
|
||||
lons[i] = pt[0];
|
||||
lats[i] = pt[1];
|
||||
}
|
||||
|
|
@ -280,64 +159,12 @@ public final class GeoUtils {
|
|||
return geometry;
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes whether a rectangle is within a given polygon (shared boundaries allowed)
|
||||
*/
|
||||
public static boolean rectWithinPoly(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double[] shapeX, final double[] shapeY, final double sMinX,
|
||||
final double sMinY, final double sMaxX, final double sMaxY) {
|
||||
// check if rectangle crosses poly (to handle concave/pacman polys), then check that all 4 corners
|
||||
// are contained
|
||||
return !(rectCrossesPoly(rMinX, rMinY, rMaxX, rMaxY, shapeX, shapeY, sMinX, sMinY, sMaxX, sMaxY) ||
|
||||
!pointInPolygon(shapeX, shapeY, rMinY, rMinX) || !pointInPolygon(shapeX, shapeY, rMinY, rMaxX) ||
|
||||
!pointInPolygon(shapeX, shapeY, rMaxY, rMaxX) || !pointInPolygon(shapeX, shapeY, rMaxY, rMinX));
|
||||
}
|
||||
|
||||
private static boolean rectAnyCornersOutsideCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters) {
|
||||
return SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 > radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 > radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 > radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 > radiusMeters;
|
||||
}
|
||||
|
||||
private static boolean rectAnyCornersInCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters) {
|
||||
return SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 <= radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 <= radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 <= radiusMeters
|
||||
|| SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 <= radiusMeters;
|
||||
}
|
||||
|
||||
public static boolean rectWithinCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters) {
|
||||
return rectAnyCornersOutsideCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters) == false;
|
||||
}
|
||||
|
||||
/**
|
||||
* Determine if a bbox (defined by minLon, minLat, maxLon, maxLat) contains the provided point (defined by lon, lat)
|
||||
* NOTE: this is basic method that does not handle dateline or pole crossing. Unwrapping must be done before
|
||||
* calling this method.
|
||||
*/
|
||||
public static boolean rectCrossesCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters) {
|
||||
return rectAnyCornersInCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters)
|
||||
|| isClosestPointOnRectWithinRange(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters);
|
||||
}
|
||||
|
||||
private static boolean isClosestPointOnRectWithinRange(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
|
||||
final double centerLon, final double centerLat, final double radiusMeters) {
|
||||
double[] closestPt = {0, 0};
|
||||
GeoDistanceUtils.closestPointOnBBox(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, closestPt);
|
||||
return SloppyMath.haversin(centerLat, centerLon, closestPt[1], closestPt[0])*1000.0 <= radiusMeters;
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute Bounding Box for a circle using WGS-84 parameters
|
||||
*/
|
||||
public static GeoRect circleToBBox(final double centerLon, final double centerLat, final double radiusMeters) {
|
||||
final double radLat = StrictMath.toRadians(centerLat);
|
||||
final double radLon = StrictMath.toRadians(centerLon);
|
||||
final double radLat = TO_RADIANS * centerLat;
|
||||
final double radLon = TO_RADIANS * centerLon;
|
||||
double radDistance = radiusMeters / GeoProjectionUtils.SEMIMAJOR_AXIS;
|
||||
double minLat = radLat - radDistance;
|
||||
double maxLat = radLat + radDistance;
|
||||
|
|
@ -345,7 +172,7 @@ public final class GeoUtils {
|
|||
double maxLon;
|
||||
|
||||
if (minLat > GeoProjectionUtils.MIN_LAT_RADIANS && maxLat < GeoProjectionUtils.MAX_LAT_RADIANS) {
|
||||
double deltaLon = StrictMath.asin(StrictMath.sin(radDistance) / StrictMath.cos(radLat));
|
||||
double deltaLon = SloppyMath.asin(SloppyMath.sin(radDistance) / SloppyMath.cos(radLat));
|
||||
minLon = radLon - deltaLon;
|
||||
if (minLon < GeoProjectionUtils.MIN_LON_RADIANS) {
|
||||
minLon += 2d * StrictMath.PI;
|
||||
|
|
@ -362,10 +189,12 @@ public final class GeoUtils {
|
|||
maxLon = GeoProjectionUtils.MAX_LON_RADIANS;
|
||||
}
|
||||
|
||||
return new GeoRect(StrictMath.toDegrees(minLon), StrictMath.toDegrees(maxLon),
|
||||
StrictMath.toDegrees(minLat), StrictMath.toDegrees(maxLat));
|
||||
return new GeoRect(TO_DEGREES * minLon, TO_DEGREES * maxLon, TO_DEGREES * minLat, TO_DEGREES * maxLat);
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute Bounding Box for a polygon using WGS-84 parameters
|
||||
*/
|
||||
public static GeoRect polyToBBox(double[] polyLons, double[] polyLats) {
|
||||
if (polyLons.length != polyLats.length) {
|
||||
throw new IllegalArgumentException("polyLons and polyLats must be equal length");
|
||||
|
|
@ -393,64 +222,16 @@ public final class GeoUtils {
|
|||
GeoUtils.unscaleLat(GeoUtils.scaleLat(minLat)), GeoUtils.unscaleLat(GeoUtils.scaleLat(maxLat)));
|
||||
}
|
||||
|
||||
/*
|
||||
/**
|
||||
* Computes whether or a 3dimensional line segment intersects or crosses a sphere
|
||||
*
|
||||
* @param lon1 longitudinal location of the line segment start point (in degrees)
|
||||
* @param lat1 latitudinal location of the line segment start point (in degrees)
|
||||
* @param alt1 altitude of the line segment start point (in degrees)
|
||||
* @param lon2 longitudinal location of the line segment end point (in degrees)
|
||||
* @param lat2 latitudinal location of the line segment end point (in degrees)
|
||||
* @param alt2 altitude of the line segment end point (in degrees)
|
||||
* @param centerLon longitudinal location of center search point (in degrees)
|
||||
* @param centerLat latitudinal location of center search point (in degrees)
|
||||
* @param centerAlt altitude of the center point (in meters)
|
||||
* @param radiusMeters search sphere radius (in meters)
|
||||
* @return whether the provided line segment is a secant of the
|
||||
* /
|
||||
// NOTE: not used for 2d at the moment. used for 3d w/ altitude (we can keep or add back)
|
||||
private static boolean lineCrossesSphere(double lon1, double lat1, double alt1, double lon2,
|
||||
double lat2, double alt2, double centerLon, double centerLat,
|
||||
double centerAlt, double radiusMeters) {
|
||||
// convert to cartesian 3d (in meters)
|
||||
double[] ecf1 = GeoProjectionUtils.llaToECF(lon1, lat1, alt1, null);
|
||||
double[] ecf2 = GeoProjectionUtils.llaToECF(lon2, lat2, alt2, null);
|
||||
double[] cntr = GeoProjectionUtils.llaToECF(centerLon, centerLat, centerAlt, null);
|
||||
|
||||
final double dX = ecf2[0] - ecf1[0];
|
||||
final double dY = ecf2[1] - ecf1[1];
|
||||
final double dZ = ecf2[2] - ecf1[2];
|
||||
final double fX = ecf1[0] - cntr[0];
|
||||
final double fY = ecf1[1] - cntr[1];
|
||||
final double fZ = ecf1[2] - cntr[2];
|
||||
|
||||
final double a = dX*dX + dY*dY + dZ*dZ;
|
||||
final double b = 2 * (fX*dX + fY*dY + fZ*dZ);
|
||||
final double c = (fX*fX + fY*fY + fZ*fZ) - (radiusMeters*radiusMeters);
|
||||
|
||||
double discrim = (b*b)-(4*a*c);
|
||||
if (discrim < 0) {
|
||||
return false;
|
||||
}
|
||||
|
||||
discrim = StrictMath.sqrt(discrim);
|
||||
final double a2 = 2*a;
|
||||
final double t1 = (-b - discrim)/a2;
|
||||
final double t2 = (-b + discrim)/a2;
|
||||
|
||||
if ( (t1 < 0 || t1 > 1) ) {
|
||||
return !(t2 < 0 || t2 > 1);
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
*/
|
||||
|
||||
* validates latitude value is within standard +/-90 coordinate bounds
|
||||
*/
|
||||
public static boolean isValidLat(double lat) {
|
||||
return Double.isNaN(lat) == false && lat >= MIN_LAT_INCL && lat <= MAX_LAT_INCL;
|
||||
}
|
||||
|
||||
/**
|
||||
* validates longitude value is within standard +/-180 coordinate bounds
|
||||
*/
|
||||
public static boolean isValidLon(double lon) {
|
||||
return Double.isNaN(lon) == false && lon >= MIN_LON_INCL && lon <= MAX_LON_INCL;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -21,6 +21,7 @@ import org.apache.lucene.document.DimensionalLatLonField;
|
|||
import org.apache.lucene.document.Document;
|
||||
import org.apache.lucene.store.Directory;
|
||||
import org.apache.lucene.util.BaseGeoPointTestCase;
|
||||
import org.apache.lucene.util.GeoDistanceUtils;
|
||||
import org.apache.lucene.util.GeoRect;
|
||||
import org.apache.lucene.util.SloppyMath;
|
||||
|
||||
|
|
@ -96,15 +97,15 @@ public class TestDimensionalQueries extends BaseGeoPointTestCase {
|
|||
|
||||
@Override
|
||||
protected Boolean circleContainsPoint(double centerLat, double centerLon, double radiusMeters, double pointLat, double pointLon) {
|
||||
double distanceKM = SloppyMath.haversin(centerLat, centerLon, pointLat, pointLon);
|
||||
boolean result = distanceKM*1000.0 <= radiusMeters;
|
||||
double distanceMeters = GeoDistanceUtils.haversin(centerLat, centerLon, pointLat, pointLon);
|
||||
boolean result = distanceMeters <= radiusMeters;
|
||||
//System.out.println(" shouldMatch? centerLon=" + centerLon + " centerLat=" + centerLat + " pointLon=" + pointLon + " pointLat=" + pointLat + " result=" + result + " distanceMeters=" + (distanceKM * 1000));
|
||||
return result;
|
||||
}
|
||||
|
||||
@Override
|
||||
protected Boolean distanceRangeContainsPoint(double centerLat, double centerLon, double minRadiusMeters, double radiusMeters, double pointLat, double pointLon) {
|
||||
final double d = SloppyMath.haversin(centerLat, centerLon, pointLat, pointLon)*1000.0;
|
||||
final double d = GeoDistanceUtils.haversin(centerLat, centerLon, pointLat, pointLon);
|
||||
return d >= minRadiusMeters && d <= radiusMeters;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -28,6 +28,7 @@ import org.apache.lucene.index.RandomIndexWriter;
|
|||
import org.apache.lucene.store.Directory;
|
||||
import org.apache.lucene.util.BaseGeoPointTestCase;
|
||||
import org.apache.lucene.util.GeoRect;
|
||||
import org.apache.lucene.util.GeoRelationUtils;
|
||||
import org.apache.lucene.util.GeoUtils;
|
||||
import org.apache.lucene.util.SloppyMath;
|
||||
import org.apache.lucene.util.TestUtil;
|
||||
|
|
@ -171,11 +172,11 @@ public class TestGeoPointQuery extends BaseGeoPointTestCase {
|
|||
}
|
||||
|
||||
if (rect.minLon < rect.maxLon) {
|
||||
return GeoUtils.bboxContains(pointLon, pointLat, rect.minLon, rect.minLat, rect.maxLon, rect.maxLat);
|
||||
return GeoRelationUtils.pointInRect(pointLon, pointLat, rect.minLon, rect.minLat, rect.maxLon, rect.maxLat);
|
||||
} else {
|
||||
// Rect crosses dateline:
|
||||
return GeoUtils.bboxContains(pointLon, pointLat, -180.0, rect.minLat, rect.maxLon, rect.maxLat)
|
||||
|| GeoUtils.bboxContains(pointLon, pointLat, rect.minLon, rect.minLat, 180.0, rect.maxLat);
|
||||
return GeoRelationUtils.pointInRect(pointLon, pointLat, -180.0, rect.minLat, rect.maxLon, rect.maxLat)
|
||||
|| GeoRelationUtils.pointInRect(pointLon, pointLat, rect.minLon, rect.minLat, 180.0, rect.maxLat);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -221,7 +222,7 @@ public class TestGeoPointQuery extends BaseGeoPointTestCase {
|
|||
}
|
||||
|
||||
public void testRectCrossesCircle() throws Exception {
|
||||
assertTrue(GeoUtils.rectCrossesCircle(-180, -90, 180, 0.0, 0.667, 0.0, 88000.0));
|
||||
assertTrue(GeoRelationUtils.rectCrossesCircle(-180, -90, 180, 0.0, 0.667, 0.0, 88000.0));
|
||||
}
|
||||
|
||||
private TopDocs geoDistanceRangeQuery(double lon, double lat, double minRadius, double maxRadius, int limit)
|
||||
|
|
@ -261,9 +262,9 @@ public class TestGeoPointQuery extends BaseGeoPointTestCase {
|
|||
double yMax = 1;//5;
|
||||
|
||||
// test cell crossing poly
|
||||
assertTrue(GeoUtils.rectCrossesPoly(xMin, yMin, xMax, yMax, px, py, xMinA, yMinA, xMaxA, yMaxA));
|
||||
assertFalse(GeoUtils.rectCrossesPoly(-5, 0, 0.000001, 5, px, py, xMin, yMin, xMax, yMax));
|
||||
assertTrue(GeoUtils.rectWithinPoly(-5, 0, -2, 5, px, py, xMin, yMin, xMax, yMax));
|
||||
assertTrue(GeoRelationUtils.rectCrossesPoly(xMin, yMin, xMax, yMax, px, py, xMinA, yMinA, xMaxA, yMaxA));
|
||||
assertFalse(GeoRelationUtils.rectCrossesPoly(-5, 0, 0.000001, 5, px, py, xMin, yMin, xMax, yMax));
|
||||
assertTrue(GeoRelationUtils.rectWithinPoly(-5, 0, -2, 5, px, py, xMin, yMin, xMax, yMax));
|
||||
}
|
||||
|
||||
public void testBBoxCrossDateline() throws Exception {
|
||||
|
|
|
|||
|
|
@ -312,7 +312,7 @@ public class TestGeoUtils extends LuceneTestCase {
|
|||
// Leaf cell: brute force check all docs that fall within this cell:
|
||||
for(int docID=0;docID<docLons.length;docID++) {
|
||||
if (cell.contains(docLons[docID], docLats[docID])) {
|
||||
double distanceMeters = SloppyMath.haversin(centerLat, centerLon, docLats[docID], docLons[docID]) * 1000.0;
|
||||
double distanceMeters = GeoDistanceUtils.haversin(centerLat, centerLon, docLats[docID], docLons[docID]);
|
||||
if (distanceMeters <= radiusMeters) {
|
||||
if (VERBOSE) {
|
||||
log.println(" check doc=" + docID + ": match!");
|
||||
|
|
@ -327,7 +327,7 @@ public class TestGeoUtils extends LuceneTestCase {
|
|||
}
|
||||
} else {
|
||||
|
||||
if (GeoUtils.rectWithinCircle(cell.minLon, cell.minLat, cell.maxLon, cell.maxLat, centerLon, centerLat, radiusMeters)) {
|
||||
if (GeoRelationUtils.rectWithinCircle(cell.minLon, cell.minLat, cell.maxLon, cell.maxLat, centerLon, centerLat, radiusMeters)) {
|
||||
// Query circle fully contains this cell, just addAll:
|
||||
if (VERBOSE) {
|
||||
log.println(" circle fully contains cell: now addAll");
|
||||
|
|
@ -341,13 +341,13 @@ public class TestGeoUtils extends LuceneTestCase {
|
|||
}
|
||||
}
|
||||
continue;
|
||||
} else if (GeoUtils.rectWithin(root.minLon, root.minLat, root.maxLon, root.maxLat,
|
||||
} else if (GeoRelationUtils.rectWithin(root.minLon, root.minLat, root.maxLon, root.maxLat,
|
||||
cell.minLon, cell.minLat, cell.maxLon, cell.maxLat)) {
|
||||
// Fall through below to "recurse"
|
||||
if (VERBOSE) {
|
||||
log.println(" cell fully contains circle: keep splitting");
|
||||
}
|
||||
} else if (GeoUtils.rectCrossesCircle(cell.minLon, cell.minLat, cell.maxLon, cell.maxLat,
|
||||
} else if (GeoRelationUtils.rectCrossesCircle(cell.minLon, cell.minLat, cell.maxLon, cell.maxLat,
|
||||
centerLon, centerLat, radiusMeters)) {
|
||||
// Fall through below to "recurse"
|
||||
if (VERBOSE) {
|
||||
|
|
@ -415,15 +415,11 @@ public class TestGeoUtils extends LuceneTestCase {
|
|||
}
|
||||
|
||||
/** Tests consistency of GeoUtils.rectWithinCircle, .rectCrossesCircle, .rectWithin and SloppyMath.haversine distance check */
|
||||
@AwaitsFix(bugUrl = "https://issues.apache.org/jira/browse/LUCENE-6846")
|
||||
public void testGeoRelations() throws Exception {
|
||||
|
||||
int numDocs = atLeast(1000);
|
||||
|
||||
// boolean useSmallRanges = random().nextBoolean();
|
||||
|
||||
// TODO: the GeoUtils APIs have bugs if you use large distances:
|
||||
boolean useSmallRanges = true;
|
||||
boolean useSmallRanges = random().nextBoolean();
|
||||
|
||||
if (VERBOSE) {
|
||||
System.out.println("TEST: " + numDocs + " docs useSmallRanges=" + useSmallRanges);
|
||||
|
|
@ -475,22 +471,29 @@ public class TestGeoUtils extends LuceneTestCase {
|
|||
|
||||
if (bbox.maxLon < bbox.minLon) {
|
||||
// Crosses dateline
|
||||
log.println(" circle crosses dateline; first right query");
|
||||
log.println(" circle crosses dateline; first left query");
|
||||
double unwrappedLon = centerLon;
|
||||
if (unwrappedLon > bbox.maxLon) {
|
||||
// unwrap left
|
||||
unwrappedLon += -360.0D;
|
||||
}
|
||||
findMatches(hits, log,
|
||||
new Cell(null,
|
||||
-180, bbox.minLat,
|
||||
bbox.maxLon, bbox.maxLat,
|
||||
0),
|
||||
centerLon, centerLat, radiusMeters,
|
||||
docLons, docLats);
|
||||
log.println(" circle crosses dateline; now left query");
|
||||
unwrappedLon, centerLat, radiusMeters, docLons, docLats);
|
||||
log.println(" circle crosses dateline; now right query");
|
||||
if (unwrappedLon < bbox.maxLon) {
|
||||
// unwrap right
|
||||
unwrappedLon += 360.0D;
|
||||
}
|
||||
findMatches(hits, log,
|
||||
new Cell(null,
|
||||
bbox.minLon, bbox.minLat,
|
||||
180, bbox.maxLat,
|
||||
0),
|
||||
centerLon, centerLat, radiusMeters,
|
||||
docLons, docLats);
|
||||
unwrappedLon, centerLat, radiusMeters, docLons, docLats);
|
||||
} else {
|
||||
// Start with the root cell that fully contains the shape:
|
||||
findMatches(hits, log,
|
||||
|
|
@ -510,15 +513,15 @@ public class TestGeoUtils extends LuceneTestCase {
|
|||
|
||||
// Done matching, now verify:
|
||||
for(int docID=0;docID<numDocs;docID++) {
|
||||
double distanceMeters = SloppyMath.haversin(centerLat, centerLon, docLats[docID], docLons[docID]) * 1000.0;
|
||||
double distanceMeters = GeoDistanceUtils.haversin(centerLat, centerLon, docLats[docID], docLons[docID]);
|
||||
boolean expected = distanceMeters <= radiusMeters;
|
||||
|
||||
boolean actual = hits.contains(docID);
|
||||
if (actual != expected) {
|
||||
if (actual) {
|
||||
log.println("doc=" + docID + " matched but should not");
|
||||
log.println("doc=" + docID + " matched but should not on iteration " + iter);
|
||||
} else {
|
||||
log.println("doc=" + docID + " did not match but should");
|
||||
log.println("doc=" + docID + " did not match but should on iteration " + iter);
|
||||
}
|
||||
log.println(" lon=" + docLons[docID] + " lat=" + docLats[docID] + " distanceMeters=" + distanceMeters + " vs radiusMeters=" + radiusMeters);
|
||||
failCount++;
|
||||
|
|
|
|||
Loading…
Reference in New Issue