LUCENE-6547: add GeoPointDistanceQuery and fix GeoPointInBBoxQuery to handle dateline crossing

git-svn-id: https://svn.apache.org/repos/asf/lucene/dev/trunk@1691659 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Michael McCandless 2015-07-18 01:00:21 +00:00
parent f917ac900f
commit ca60372d47
13 changed files with 1524 additions and 432 deletions

View File

@ -137,6 +137,10 @@ New Features
* LUCENE-6659: Remove IndexWriter's unnecessary hard limit on max concurrency
(Robert Muir, Mike McCandless)
* LUCENE-6547: Add GeoPointDistanceQuery, matching all points within
the specified distance from the center point. Fix
GeoPointInBBoxQuery to handle dateline crossing.
API Changes
* LUCENE-6508: Simplify Lock api, there is now just

View File

@ -32,7 +32,8 @@ import org.apache.lucene.util.GeoUtils;
* </pre>
*
* <p>To perform simple geospatial queries against a <code>GeoPointField</code>,
* see {@link org.apache.lucene.search.GeoPointInBBoxQuery} or {@link org.apache.lucene.search.GeoPointInPolygonQuery}
* see {@link org.apache.lucene.search.GeoPointInBBoxQuery}, {@link org.apache.lucene.search.GeoPointInPolygonQuery},
* or {@link org.apache.lucene.search.GeoPointDistanceQuery}
*
* NOTE: This indexes only high precision encoded terms which may result in visiting a high number
* of terms for large queries. See LUCENE-6481 for a future improvement.

View File

@ -0,0 +1,169 @@
package org.apache.lucene.search;
/*
* 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.
*/
import org.apache.lucene.index.IndexReader;
import org.apache.lucene.util.GeoDistanceUtils;
import org.apache.lucene.util.GeoUtils;
import org.apache.lucene.util.ToStringUtils;
/** Implements a simple point distance query on a GeoPoint field. This is based on
* {@link org.apache.lucene.search.GeoPointInBBoxQuery} and is implemented using a two phase approach. First,
* like {@code GeoPointInBBoxQueryImpl} candidate terms are queried using the numeric ranges based on
* the morton codes of the min and max lat/lon pairs that intersect the boundary of the point-radius
* circle (see {@link org.apache.lucene.util.GeoUtils#lineCrossesSphere}. Terms
* passing this initial filter are then passed to a secondary filter that verifies whether the
* decoded lat/lon point fall within the specified query distance (see {@link org.apache.lucene.util.SloppyMath#haversin}.
* All morton value comparisons are subject to the same precision tolerance defined in
* {@value org.apache.lucene.util.GeoUtils#TOLERANCE} and distance comparisons are subject to the accuracy of the
* haversine formula (from R.W. Sinnott, "Virtues of the Haversine", Sky and Telescope, vol. 68, no. 2, 1984, p. 159)
*
* NOTE: This query works best for point-radius queries that do not cross the dateline, there is a penalty for crossing
* the dateline as the bounding box is effectively split into two separate queries, and the point-radius is converted
* to a euclidean spherical search to handle a wrapping coordinate system (TODO split the point radius at the dateline?)
*
* This query also currently uses haversine which is a sloppy distance calculation. For large queries one can expect
* upwards of 400m error. Vincenty shrinks this to ~40m error but pays a penalty for computing using the spheroid
*
* @lucene.experimental
*/
public final class GeoPointDistanceQuery extends GeoPointInBBoxQuery {
protected final double centerLon;
protected final double centerLat;
protected final double radius;
/** NOTE: radius is in meters. */
public GeoPointDistanceQuery(final String field, final double centerLon, final double centerLat, final double radius) {
this(field, computeBBox(centerLon, centerLat, radius), centerLon, centerLat, radius);
}
private GeoPointDistanceQuery(final String field, GeoBoundingBox bbox, final double centerLon,
final double centerLat, final double radius) {
super(field, bbox.minLon, bbox.minLat, bbox.maxLon, bbox.maxLat);
if (GeoUtils.isValidLon(centerLon) == false) {
throw new IllegalArgumentException("invalid centerLon " + centerLon);
}
if (GeoUtils.isValidLat(centerLat) == false) {
throw new IllegalArgumentException("invalid centerLat " + centerLat);
}
this.centerLon = centerLon;
this.centerLat = centerLat;
this.radius = radius;
}
@Override
public Query rewrite(IndexReader reader) {
if (maxLon < minLon) {
BooleanQuery.Builder bqb = new BooleanQuery.Builder();
GeoPointDistanceQueryImpl left = new GeoPointDistanceQueryImpl(field, this, new GeoBoundingBox(-180.0D, maxLon,
minLat, maxLat));
left.setBoost(getBoost());
bqb.add(new BooleanClause(left, BooleanClause.Occur.SHOULD));
GeoPointDistanceQueryImpl right = new GeoPointDistanceQueryImpl(field, this, new GeoBoundingBox(minLon, 180.0D,
minLat, maxLat));
right.setBoost(getBoost());
bqb.add(new BooleanClause(right, BooleanClause.Occur.SHOULD));
return bqb.build();
}
return new GeoPointDistanceQueryImpl(field, this, new GeoBoundingBox(this.minLon, this.maxLon, this.minLat, this.maxLat));
}
static GeoBoundingBox computeBBox(final double centerLon, final double centerLat, final double radius) {
final double lonDistDeg = GeoDistanceUtils.distanceToDegreesLon(centerLat, radius);
final double latDistDeg = GeoDistanceUtils.distanceToDegreesLat(centerLat, radius);
return new GeoBoundingBox(GeoUtils.normalizeLon(centerLon - lonDistDeg), GeoUtils.normalizeLon(centerLon + lonDistDeg),
GeoUtils.normalizeLat(centerLat - latDistDeg), GeoUtils.normalizeLat(centerLat + latDistDeg));
}
@Override
public boolean equals(Object o) {
if (this == o) return true;
if (!(o instanceof GeoPointDistanceQuery)) return false;
if (!super.equals(o)) return false;
GeoPointDistanceQuery that = (GeoPointDistanceQuery) o;
if (Double.compare(that.centerLat, centerLat) != 0) return false;
if (Double.compare(that.centerLon, centerLon) != 0) return false;
if (Double.compare(that.radius, radius) != 0) return false;
return true;
}
@Override
public int hashCode() {
int result = super.hashCode();
long temp;
temp = Double.doubleToLongBits(centerLon);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(centerLat);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(radius);
result = 31 * result + (int) (temp ^ (temp >>> 32));
return result;
}
@Override
public String toString(String field) {
final StringBuilder sb = new StringBuilder();
sb.append(getClass().getSimpleName());
sb.append(':');
if (!this.field.equals(field)) {
sb.append(" field=");
sb.append(this.field);
sb.append(':');
}
return sb.append( " Center: [")
.append(centerLon)
.append(',')
.append(centerLat)
.append(']')
.append(" Distance: ")
.append(radius)
.append(" m")
.append(" Lower Left: [")
.append(minLon)
.append(',')
.append(minLat)
.append(']')
.append(" Upper Right: [")
.append(maxLon)
.append(',')
.append(maxLat)
.append("]")
.append(ToStringUtils.boost(getBoost()))
.toString();
}
public double getCenterLon() {
return this.centerLon;
}
public double getCenterLat() {
return this.centerLat;
}
public double getRadius() {
return this.radius;
}
}

View File

@ -0,0 +1,119 @@
package org.apache.lucene.search;
/*
* 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.
*/
import java.io.IOException;
import org.apache.lucene.index.Terms;
import org.apache.lucene.index.TermsEnum;
import org.apache.lucene.util.AttributeSource;
import org.apache.lucene.util.BytesRef;
import org.apache.lucene.util.GeoUtils;
import org.apache.lucene.util.NumericUtils;
import org.apache.lucene.util.SloppyMath;
/** Package private implementation for the public facing GeoPointDistanceQuery delegate class.
*
* @lucene.experimental
*/
final class GeoPointDistanceQueryImpl extends GeoPointInBBoxQueryImpl {
private final GeoPointDistanceQuery query;
GeoPointDistanceQueryImpl(final String field, final GeoPointDistanceQuery q, final GeoBoundingBox bbox) {
super(field, bbox.minLon, bbox.minLat, bbox.maxLon, bbox.maxLat);
query = q;
}
@Override @SuppressWarnings("unchecked")
protected TermsEnum getTermsEnum(final Terms terms, AttributeSource atts) throws IOException {
return new GeoPointRadiusTermsEnum(terms.iterator(), this.minLon, this.minLat, this.maxLon, this.maxLat);
}
@Override
public void setRewriteMethod(RewriteMethod method) {
throw new UnsupportedOperationException("cannot change rewrite method");
}
private final class GeoPointRadiusTermsEnum extends GeoPointTermsEnum {
GeoPointRadiusTermsEnum(final TermsEnum tenum, final double minLon, final double minLat,
final double maxLon, final double maxLat) {
super(tenum, minLon, minLat, maxLon, maxLat);
}
@Override
protected boolean cellCrosses(final double minLon, final double minLat, final double maxLon, final double maxLat) {
return GeoUtils.rectCrossesCircle(minLon, minLat, maxLon, maxLat, query.centerLon, query.centerLat, query.radius);
}
@Override
protected boolean cellWithin(final double minLon, final double minLat, final double maxLon, final double maxLat) {
return GeoUtils.rectWithinCircle(minLon, minLat, maxLon, maxLat, query.centerLon, query.centerLat, query.radius);
}
/**
* The two-phase query approach. The parent
* {@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.
*
* @param term term for candidate document
* @return match status
*/
@Override
protected AcceptStatus accept(BytesRef term) {
// first filter by bounding box
AcceptStatus status = super.accept(term);
assert status != AcceptStatus.YES_AND_SEEK;
if (status != AcceptStatus.YES) {
return status;
}
final long val = NumericUtils.prefixCodedToLong(term);
final double lon = GeoUtils.mortonUnhashLon(val);
final double lat = GeoUtils.mortonUnhashLat(val);
// post-filter by distance
if (!(SloppyMath.haversin(query.centerLat, query.centerLon, lat, lon)*1000.0 <= query.radius)) {
return AcceptStatus.NO;
}
return AcceptStatus.YES;
}
}
@Override
public boolean equals(Object o) {
if (this == o) return true;
if (!(o instanceof GeoPointDistanceQueryImpl)) return false;
if (!super.equals(o)) return false;
GeoPointDistanceQueryImpl that = (GeoPointDistanceQueryImpl) o;
if (!query.equals(that.query)) return false;
return true;
}
@Override
public int hashCode() {
int result = super.hashCode();
result = 31 * result + query.hashCode();
return result;
}
}

View File

@ -17,12 +17,7 @@ package org.apache.lucene.search;
* limitations under the License.
*/
import java.io.IOException;
import org.apache.lucene.index.Terms;
import org.apache.lucene.index.TermsEnum;
import org.apache.lucene.util.AttributeSource;
import org.apache.lucene.util.GeoUtils;
import org.apache.lucene.index.IndexReader;
import org.apache.lucene.util.ToStringUtils;
/** Implements a simple bounding box query on a GeoPoint field. This is inspired by
@ -44,82 +39,35 @@ import org.apache.lucene.util.ToStringUtils;
*
* @lucene.experimental
*/
public class GeoPointInBBoxQuery extends MultiTermQuery {
// simple bounding box optimization - no objects used to avoid dependencies
public class GeoPointInBBoxQuery extends Query {
protected final String field;
protected final double minLon;
protected final double minLat;
protected final double maxLon;
protected final double maxLat;
/**
* Constructs a new GeoBBoxQuery that will match encoded GeoPoint terms that fall within or on the boundary
* of the bounding box defined by the input parameters
* @param field the field name
* @param minLon lower longitude (x) value of the bounding box
* @param minLat lower latitude (y) value of the bounding box
* @param maxLon upper longitude (x) value of the bounding box
* @param maxLat upper latitude (y) value of the bounding box
*/
public GeoPointInBBoxQuery(final String field, final double minLon, final double minLat, final double maxLon, final double maxLat) {
super(field);
if (GeoUtils.isValidLon(minLon) == false) {
throw new IllegalArgumentException("invalid minLon " + minLon);
}
if (GeoUtils.isValidLon(maxLon) == false) {
throw new IllegalArgumentException("invalid maxLon " + maxLon);
}
if (GeoUtils.isValidLat(minLat) == false) {
throw new IllegalArgumentException("invalid minLat " + minLat);
}
if (GeoUtils.isValidLat(maxLat) == false) {
throw new IllegalArgumentException("invalid maxLat " + maxLat);
}
this.field = field;
this.minLon = minLon;
this.minLat = minLat;
this.maxLon = maxLon;
this.maxLat = maxLat;
}
@Override @SuppressWarnings("unchecked")
protected TermsEnum getTermsEnum(final Terms terms, AttributeSource atts) throws IOException {
final Long min = GeoUtils.mortonHash(minLon, minLat);
final Long max = Math.abs(GeoUtils.mortonHash(maxLon, maxLat));
if (min != null && max != null && min.compareTo(max) > 0) {
return TermsEnum.EMPTY;
@Override
public Query rewrite(IndexReader reader) {
if (maxLon < minLon) {
BooleanQuery.Builder bqb = new BooleanQuery.Builder();
GeoPointInBBoxQueryImpl left = new GeoPointInBBoxQueryImpl(field, -180.0D, minLat, maxLon, maxLat);
left.setBoost(getBoost());
bqb.add(new BooleanClause(left, BooleanClause.Occur.SHOULD));
GeoPointInBBoxQueryImpl right = new GeoPointInBBoxQueryImpl(field, minLon, minLat, 180.0D, maxLat);
right.setBoost(getBoost());
bqb.add(new BooleanClause(right, BooleanClause.Occur.SHOULD));
return bqb.build();
}
return new GeoPointTermsEnum(terms.iterator(), atts, minLon, minLat, maxLon, maxLat);
}
@Override
@SuppressWarnings({"unchecked","rawtypes"})
public boolean equals(Object o) {
if (this == o) return true;
if (o == null || getClass() != o.getClass()) return false;
if (!super.equals(o)) return false;
GeoPointInBBoxQuery that = (GeoPointInBBoxQuery) o;
if (Double.compare(that.maxLat, maxLat) != 0) return false;
if (Double.compare(that.maxLon, maxLon) != 0) return false;
if (Double.compare(that.minLat, minLat) != 0) return false;
if (Double.compare(that.minLon, minLon) != 0) return false;
return true;
}
@Override
public int hashCode() {
int result = super.hashCode();
long temp;
temp = Double.doubleToLongBits(minLon);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(minLat);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(maxLon);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(maxLat);
result = 31 * result + (int) (temp ^ (temp >>> 32));
return result;
return new GeoPointInBBoxQueryImpl(field, minLon, minLat, maxLon, maxLat);
}
@Override
@ -127,9 +75,9 @@ public class GeoPointInBBoxQuery extends MultiTermQuery {
final StringBuilder sb = new StringBuilder();
sb.append(getClass().getSimpleName());
sb.append(':');
if (!getField().equals(field)) {
if (!this.field.equals(field)) {
sb.append(" field=");
sb.append(getField());
sb.append(this.field);
sb.append(':');
}
return sb.append(" Lower Left: [")
@ -145,4 +93,57 @@ public class GeoPointInBBoxQuery extends MultiTermQuery {
.append(ToStringUtils.boost(getBoost()))
.toString();
}
@Override
public boolean equals(Object o) {
if (this == o) return true;
if (!(o instanceof GeoPointInBBoxQuery)) return false;
if (!super.equals(o)) return false;
GeoPointInBBoxQuery that = (GeoPointInBBoxQuery) o;
if (Double.compare(that.maxLat, maxLat) != 0) return false;
if (Double.compare(that.maxLon, maxLon) != 0) return false;
if (Double.compare(that.minLat, minLat) != 0) return false;
if (Double.compare(that.minLon, minLon) != 0) return false;
if (!field.equals(that.field)) return false;
return true;
}
@Override
public int hashCode() {
int result = super.hashCode();
long temp;
result = 31 * result + field.hashCode();
temp = Double.doubleToLongBits(minLon);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(minLat);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(maxLon);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(maxLat);
result = 31 * result + (int) (temp ^ (temp >>> 32));
return result;
}
public final String getField() {
return this.field;
}
public final double getMinLon() {
return this.minLon;
}
public final double getMinLat() {
return this.minLat;
}
public final double getMaxLon() {
return this.maxLon;
}
public final double getMaxLat() {
return this.maxLat;
}
}

View File

@ -0,0 +1,110 @@
package org.apache.lucene.search;
/*
* 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.
*/
import java.io.IOException;
import org.apache.lucene.index.Terms;
import org.apache.lucene.index.TermsEnum;
import org.apache.lucene.util.AttributeSource;
import org.apache.lucene.util.ToStringUtils;
/** Package private implementation for the public facing GeoPointInBBoxQuery delegate class.
*
* @lucene.experimental
*/
class GeoPointInBBoxQueryImpl extends GeoPointTermQuery {
/**
* Constructs a new GeoBBoxQuery that will match encoded GeoPoint terms that fall within or on the boundary
* of the bounding box defined by the input parameters
* @param field the field name
* @param minLon lower longitude (x) value of the bounding box
* @param minLat lower latitude (y) value of the bounding box
* @param maxLon upper longitude (x) value of the bounding box
* @param maxLat upper latitude (y) value of the bounding box
*/
GeoPointInBBoxQueryImpl(final String field, final double minLon, final double minLat, final double maxLon, final double maxLat) {
super(field, minLon, minLat, maxLon, maxLat);
}
@Override @SuppressWarnings("unchecked")
protected TermsEnum getTermsEnum(final Terms terms, AttributeSource atts) throws IOException {
return new GeoPointTermsEnum(terms.iterator(), minLon, minLat, maxLon, maxLat);
}
@Override
public void setRewriteMethod(RewriteMethod method) {
throw new UnsupportedOperationException("cannot change rewrite method");
}
@Override
@SuppressWarnings({"unchecked","rawtypes"})
public boolean equals(Object o) {
if (this == o) return true;
if (o == null || getClass() != o.getClass()) return false;
if (!super.equals(o)) return false;
GeoPointInBBoxQueryImpl that = (GeoPointInBBoxQueryImpl) o;
if (Double.compare(that.maxLat, maxLat) != 0) return false;
if (Double.compare(that.maxLon, maxLon) != 0) return false;
if (Double.compare(that.minLat, minLat) != 0) return false;
if (Double.compare(that.minLon, minLon) != 0) return false;
return true;
}
@Override
public int hashCode() {
int result = super.hashCode();
long temp;
temp = Double.doubleToLongBits(minLon);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(minLat);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(maxLon);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(maxLat);
result = 31 * result + (int) (temp ^ (temp >>> 32));
return result;
}
@Override
public String toString(String field) {
final StringBuilder sb = new StringBuilder();
sb.append(getClass().getSimpleName());
sb.append(':');
if (!getField().equals(field)) {
sb.append(" field=");
sb.append(getField());
sb.append(':');
}
return sb.append(" Lower Left: [")
.append(minLon)
.append(',')
.append(minLat)
.append(']')
.append(" Upper Right: [")
.append(maxLon)
.append(',')
.append(maxLat)
.append("]")
.append(ToStringUtils.boost(getBoost()))
.toString();
}
}

View File

@ -29,8 +29,8 @@ import java.io.IOException;
import java.util.Arrays;
/** Implements a simple point in polygon query on a GeoPoint field. This is based on
* {@link GeoPointInBBoxQuery} and is implemented using a
* three phase approach. First, like {@link GeoPointInBBoxQuery}
* {@code GeoPointInBBoxQueryImpl} and is implemented using a
* three phase approach. First, like {@code GeoPointInBBoxQueryImpl}
* candidate terms are queried using a numeric range based on the morton codes
* of the min and max lat/lon pairs. Terms passing this initial filter are passed
* to a secondary filter that verifies whether the decoded lat/lon point falls within
@ -47,7 +47,7 @@ import java.util.Arrays;
*
* @lucene.experimental
*/
public final class GeoPointInPolygonQuery extends GeoPointInBBoxQuery {
public final class GeoPointInPolygonQuery extends GeoPointInBBoxQueryImpl {
// polygon position arrays - this avoids the use of any objects or
// or geo library dependencies
private final double[] x;
@ -61,27 +61,6 @@ public final class GeoPointInPolygonQuery extends GeoPointInBBoxQuery {
this(field, computeBBox(polyLons, polyLats), polyLons, polyLats);
}
/**
* Expert: constructs a new GeoPolygonQuery that will match encoded {@link org.apache.lucene.document.GeoPointField} terms
* that fall within or on the boundary of the polygon defined by the input parameters. This constructor requires a
* precomputed bounding box. As an alternative, {@link #GeoPointInPolygonQuery(String,double[],double[])} can
* be used to compute the bounding box during construction
*
* @param field the field name
* @param minLon lower longitude (x) value of the bounding box optimizer
* @param minLat lower latitude (y) value of the bounding box optimizer
* @param maxLon upper longitude (x) value of the bounding box optimizer
* @param maxLat upper latitude (y) value of the bounding box optimizer
* @param polyLons array containing all longitude values for the polygon
* @param polyLats array containing all latitude values for the polygon
*/
public GeoPointInPolygonQuery(final String field, final double minLon, final double minLat, final double maxLon,
final double maxLat, final double[] polyLons, final double[] polyLats) {
// TODO: should we remove this? It's dangerous .. app could accidentally provide too-small bbox?
// we should at least verify that bbox does in fact fully contain the poly?
this(field, new GeoBoundingBox(minLon, maxLon, minLat, maxLat), polyLons, polyLats);
}
/** Common constructor, used only internally. */
private GeoPointInPolygonQuery(final String field, GeoBoundingBox bbox, final double[] polyLons, final double[] polyLats) {
super(field, bbox.minLon, bbox.minLat, bbox.maxLon, bbox.maxLat);
@ -112,16 +91,16 @@ public final class GeoPointInPolygonQuery extends GeoPointInBBoxQuery {
@Override @SuppressWarnings("unchecked")
protected TermsEnum getTermsEnum(final Terms terms, AttributeSource atts) throws IOException {
final Long min = GeoUtils.mortonHash(minLon, minLat);
final Long max = Math.abs(GeoUtils.mortonHash(maxLon, maxLat));
if (min != null && max != null && min.compareTo(max) > 0) {
return TermsEnum.EMPTY;
}
return new GeoPolygonTermsEnum(terms.iterator(), atts, minLon, minLat, maxLon, maxLat);
return new GeoPolygonTermsEnum(terms.iterator(), this.minLon, this.minLat, this.maxLon, this.maxLat);
}
@Override
public final boolean equals(Object o) {
public void setRewriteMethod(RewriteMethod method) {
throw new UnsupportedOperationException("cannot change rewrite method");
}
@Override
public boolean equals(Object o) {
if (this == o) return true;
if (o == null || getClass() != o.getClass()) return false;
if (!super.equals(o)) return false;
@ -135,7 +114,7 @@ public final class GeoPointInPolygonQuery extends GeoPointInBBoxQuery {
}
@Override
public final int hashCode() {
public int hashCode() {
int result = super.hashCode();
result = 31 * result + (x != null ? Arrays.hashCode(x) : 0);
result = 31 * result + (y != null ? Arrays.hashCode(y) : 0);
@ -167,10 +146,14 @@ public final class GeoPointInPolygonQuery extends GeoPointInBBoxQuery {
return sb.toString();
}
/**
* Custom {@link org.apache.lucene.index.TermsEnum} that computes morton hash ranges based on the defined edges of
* the provided polygon.
*/
private final class GeoPolygonTermsEnum extends GeoPointTermsEnum {
GeoPolygonTermsEnum(final TermsEnum tenum, AttributeSource atts, final double minLon, final double minLat,
GeoPolygonTermsEnum(final TermsEnum tenum, final double minLon, final double minLat,
final double maxLon, final double maxLat) {
super(tenum, atts, minLon, minLat, maxLon, maxLat);
super(tenum, minLon, minLat, maxLon, maxLat);
}
@Override
@ -202,7 +185,7 @@ public final class GeoPointInPolygonQuery extends GeoPointInBBoxQuery {
* @return match status
*/
@Override
protected final AcceptStatus accept(BytesRef term) {
protected AcceptStatus accept(BytesRef term) {
// first filter by bounding box
AcceptStatus status = super.accept(term);
assert status != AcceptStatus.YES_AND_SEEK;
@ -247,4 +230,20 @@ public final class GeoPointInPolygonQuery extends GeoPointInBBoxQuery {
return new GeoBoundingBox(minLon, maxLon, minLat, maxLat);
}
/**
* API utility method for returning the array of longitudinal values for this GeoPolygon
* The returned array is not a copy so do not change it!
*/
public double[] getLons() {
return this.x;
}
/**
* API utility method for returning the array of latitudinal values for this GeoPolygon
* The returned array is not a copy so do not change it!
*/
public double[] getLats() {
return this.y;
}
}

View File

@ -0,0 +1,59 @@
package org.apache.lucene.search;
/*
* 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.
*/
import org.apache.lucene.util.GeoUtils;
/**
* Abstract base class used by {@code GeoPointInBBoxQueryImpl}
*
* @lucene.experimental
*/
// TODO: remove this? Just absorb into its base class
abstract class GeoPointTermQuery extends MultiTermQuery {
// simple bounding box optimization - no objects used to avoid dependencies
protected final double minLon;
protected final double minLat;
protected final double maxLon;
protected final double maxLat;
/**
* Constructs a query matching terms that cannot be represented with a single
* Term.
*/
public GeoPointTermQuery(String field, final double minLon, final double minLat, final double maxLon, final double maxLat) {
super(field);
if (GeoUtils.isValidLon(minLon) == false) {
throw new IllegalArgumentException("invalid minLon " + minLon);
}
if (GeoUtils.isValidLon(maxLon) == false) {
throw new IllegalArgumentException("invalid maxLon " + maxLon);
}
if (GeoUtils.isValidLat(minLat) == false) {
throw new IllegalArgumentException("invalid minLat " + minLat);
}
if (GeoUtils.isValidLat(maxLat) == false) {
throw new IllegalArgumentException("invalid maxLat " + maxLat);
}
this.minLon = minLon;
this.minLat = minLat;
this.maxLon = maxLon;
this.maxLat = maxLat;
}
}

View File

@ -21,14 +21,9 @@ import java.util.Collections;
import java.util.LinkedList;
import java.util.List;
import org.apache.lucene.analysis.tokenattributes.TermToBytesRefAttribute;
import org.apache.lucene.document.GeoPointField;
import org.apache.lucene.index.FilteredTermsEnum;
import org.apache.lucene.index.TermsEnum;
import org.apache.lucene.util.Attribute;
import org.apache.lucene.util.AttributeImpl;
import org.apache.lucene.util.AttributeReflector;
import org.apache.lucene.util.AttributeSource;
import org.apache.lucene.util.BytesRef;
import org.apache.lucene.util.BytesRefBuilder;
import org.apache.lucene.util.GeoUtils;
@ -47,13 +42,11 @@ class GeoPointTermsEnum extends FilteredTermsEnum {
private Range currentRange;
private BytesRef currentLowerBound, currentUpperBound;
private final ComputedRangesAttribute rangesAtt;
private final List<Range> rangeBounds = new LinkedList<>();
private final LinkedList<Range> rangeBounds;
protected static final short DETAIL_LEVEL = 16;
private static final short DETAIL_LEVEL = 16;
GeoPointTermsEnum(final TermsEnum tenum, AttributeSource atts, final double minLon, final double minLat,
GeoPointTermsEnum(final TermsEnum tenum, final double minLon, final double minLat,
final double maxLon, final double maxLat) {
super(tenum);
final long rectMinHash = GeoUtils.mortonHash(minLon, minLat);
@ -63,13 +56,8 @@ class GeoPointTermsEnum extends FilteredTermsEnum {
this.maxLon = GeoUtils.mortonUnhashLon(rectMaxHash);
this.maxLat = GeoUtils.mortonUnhashLat(rectMaxHash);
this.rangesAtt = atts.addAttribute(ComputedRangesAttribute.class);
this.rangeBounds = rangesAtt.ranges();
if (rangeBounds.isEmpty()) {
computeRange(0L, (short) (((GeoUtils.BITS) << 1) - 1));
Collections.sort(rangeBounds);
}
computeRange(0L, (short) (((GeoUtils.BITS) << 1) - 1));
Collections.sort(rangeBounds);
}
/**
@ -100,9 +88,7 @@ class GeoPointTermsEnum extends FilteredTermsEnum {
final short level = (short)(62-res>>>1);
// if cell is within and a factor of the precision step, add the range
// if cell cellCrosses
// if cell is within and a factor of the precision step, or it crosses the edge of the shape add the range
final boolean within = res% GeoPointField.PRECISION_STEP == 0 && cellWithin(minLon, minLat, maxLon, maxLat);
if (within || (level == DETAIL_LEVEL && cellCrosses(minLon, minLat, maxLon, maxLat))) {
rangeBounds.add(new Range(start, end, res, level, !within));
@ -124,7 +110,7 @@ class GeoPointTermsEnum extends FilteredTermsEnum {
}
private void nextRange() {
currentRange = rangeBounds.removeFirst();
currentRange = rangeBounds.remove(0);
currentLowerBound = currentRange.lower;
assert currentUpperBound == null || currentUpperBound.compareTo(currentRange.lower) <= 0 :
"The current upper bound must be <= the new lower bound";
@ -169,11 +155,13 @@ class GeoPointTermsEnum extends FilteredTermsEnum {
protected AcceptStatus accept(BytesRef term) {
// validate value is in range
while (currentUpperBound == null || term.compareTo(currentUpperBound) > 0) {
if (rangeBounds.isEmpty())
if (rangeBounds.isEmpty()) {
return AcceptStatus.END;
}
// peek next sub-range, only seek if the current term is smaller than next lower bound
if (term.compareTo(rangeBounds.getFirst().lower) < 0)
if (term.compareTo(rangeBounds.get(0).lower) < 0) {
return AcceptStatus.NO_AND_SEEK;
}
// step forward to next range without seeking, as next lower range bound is less or equal current term
nextRange();
}
@ -190,62 +178,10 @@ class GeoPointTermsEnum extends FilteredTermsEnum {
return AcceptStatus.YES;
}
public static interface ComputedRangesAttribute extends Attribute {
public LinkedList<Range> ranges();
}
@SuppressWarnings({"unchecked","rawtypes"})
public static final class ComputedRangesAttributeImpl extends AttributeImpl implements ComputedRangesAttribute {
public final LinkedList<Range> rangeBounds = new LinkedList();
@Override
public LinkedList<Range> ranges() {
return rangeBounds;
}
@Override
public void clear() {
rangeBounds.clear();
}
@Override
public int hashCode() {
return rangeBounds.hashCode();
}
@Override
public boolean equals(Object other) {
if (this == other)
return true;
if (!(other instanceof ComputedRangesAttributeImpl))
return false;
return rangeBounds.equals(((ComputedRangesAttributeImpl)other).rangeBounds);
}
@Override
public void copyTo(AttributeImpl target) {
final List<Range> targetRanges = ((ComputedRangesAttribute)target).ranges();
targetRanges.clear();
targetRanges.addAll(rangeBounds);
}
@Override
public AttributeImpl clone() {
ComputedRangesAttributeImpl c = (ComputedRangesAttributeImpl) super.clone();;
copyTo(c);
return c;
}
@Override
public void reflectWith(AttributeReflector reflector) {
reflector.reflect(ComputedRangesAttribute.class, "rangeBounds", rangeBounds);
}
}
/**
* Internal class to represent a range along the space filling curve
*/
private final class Range implements Comparable<Range> {
protected final class Range implements Comparable<Range> {
final BytesRef lower;
final BytesRef upper;
final short level;
@ -263,7 +199,7 @@ class GeoPointTermsEnum extends FilteredTermsEnum {
}
@Override
public final int compareTo(Range other) {
public int compareTo(Range other) {
return this.lower.compareTo(other.lower);
}
}

View File

@ -0,0 +1,130 @@
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-spatial distance utility methods.
*
* @lucene.experimental
*/
public class GeoDistanceUtils {
/**
* 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
* 22m better accuracy (in worst case) than haversine
*
* @param lonA longitudinal coordinate of point A (in degrees)
* @param latA latitudinal coordinate of point A (in degrees)
* @param lonB longitudinal coordinate of point B (in degrees)
* @param latB latitudinal coordinate of point B (in degrees)
* @return distance (in meters) between point A and point B
*/
public static final double vincentyDistance(final double lonA, final double latA, final double lonB, final double latB) {
final double L = StrictMath.toRadians(lonB - lonA);
final double oF = 1 - GeoProjectionUtils.FLATTENING;
final double U1 = StrictMath.atan(oF * StrictMath.tan(StrictMath.toRadians(latA)));
final double U2 = StrictMath.atan(oF * StrictMath.tan(StrictMath.toRadians(latB)));
final double sU1 = StrictMath.sin(U1);
final double cU1 = StrictMath.cos(U1);
final double sU2 = StrictMath.sin(U2);
final double cU2 = StrictMath.cos(U2);
double sigma, sinSigma, cosSigma;
double sinAlpha, cos2Alpha, cos2SigmaM;
double lambda = L;
double lambdaP;
double iters = 100;
double sinLambda, cosLambda, c;
do {
sinLambda = StrictMath.sin(lambda);
cosLambda = Math.cos(lambda);
sinSigma = Math.sqrt((cU2 * sinLambda) * (cU2 * sinLambda) + (cU1 * sU2 - sU1 * cU2 * cosLambda)
* (cU1 * sU2 - sU1 * cU2 * cosLambda));
if (sinSigma == 0) {
return 0;
}
cosSigma = sU1 * sU2 + cU1 * cU2 * cosLambda;
sigma = Math.atan2(sinSigma, cosSigma);
sinAlpha = cU1 * cU2 * sinLambda / sinSigma;
cos2Alpha = 1 - sinAlpha * sinAlpha;
cos2SigmaM = cosSigma - 2 * sU1 * sU2 / cos2Alpha;
c = GeoProjectionUtils.FLATTENING/16 * cos2Alpha * (4 + GeoProjectionUtils.FLATTENING * (4 - 3 * cos2Alpha));
lambdaP = lambda;
lambda = L + (1 - c) * GeoProjectionUtils.FLATTENING * sinAlpha * (sigma + c * sinSigma * (cos2SigmaM + c * cosSigma *
(-1 + 2 * cos2SigmaM * cos2SigmaM)));
} while (StrictMath.abs(lambda - lambdaP) > 1E-12 && --iters > 0);
if (iters == 0) {
return 0;
}
final double uSq = cos2Alpha * (GeoProjectionUtils.SEMIMAJOR_AXIS2 - GeoProjectionUtils.SEMIMINOR_AXIS2) / (GeoProjectionUtils.SEMIMINOR_AXIS2);
final double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
final double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
final double deltaSigma = B * sinSigma * (cos2SigmaM + B/4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B/6 * cos2SigmaM
* (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
return (GeoProjectionUtils.SEMIMINOR_AXIS * A * (sigma - deltaSigma));
}
/**
* Compute the inverse haversine to determine distance in degrees longitude for provided distance in meters
* @param lat latitude to compute delta degrees lon
* @param distance distance in meters to convert to degrees lon
* @return Sloppy distance in degrees longitude for provided distance in meters
*/
public static double distanceToDegreesLon(double lat, double distance) {
distance /= 1000.0;
// convert latitude to radians
lat = StrictMath.toRadians(lat);
// get the diameter at the latitude
final double diameter = SloppyMath.earthDiameter(StrictMath.toRadians(lat));
// compute inverse haversine
double a = StrictMath.sin(distance/diameter);
double h = StrictMath.min(1, a);
h *= h;
double cLat = StrictMath.cos(lat);
return StrictMath.toDegrees(StrictMath.acos(1-((2d*h)/(cLat*cLat))));
}
/**
* Compute the inverse haversine to determine distance in degrees longitude for provided distance in meters
* @param lat latitude to compute delta degrees lon
* @param distance distance in meters to convert to degrees lon
* @return Sloppy distance in degrees longitude for provided distance in meters
*/
public static double distanceToDegreesLat(double lat, double distance) {
// get the diameter at the latitude
final double diameter = SloppyMath.earthDiameter(StrictMath.toRadians(lat));
distance /= 1000.0;
// compute inverse haversine
double a = StrictMath.sin(distance/diameter);
double h = StrictMath.min(1, a);
h *= h;
return StrictMath.toDegrees(StrictMath.acos(1-(2d*h)));
}
}

View File

@ -0,0 +1,383 @@
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-spatial projection utility methods.
*
* @lucene.experimental
*/
public class GeoProjectionUtils {
// WGS84 earth-ellipsoid major (a) minor (b) radius, (f) flattening and eccentricity (e)
static final double SEMIMAJOR_AXIS = 6_378_137; // [m]
static final double FLATTENING = 1.0/298.257223563;
static final double SEMIMINOR_AXIS = SEMIMAJOR_AXIS * (1.0 - FLATTENING); //6_356_752.31420; // [m]
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;
/**
* Converts from geocentric earth-centered earth-fixed to geodesic lat/lon/alt
* @param x Cartesian x coordinate
* @param y Cartesian y coordinate
* @param z Cartesian z coordinate
* @param lla 0: longitude 1: latitude: 2: altitude
* @return double array as 0: longitude 1: latitude 2: altitude
*/
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) {
lla = new double[3];
}
if (x != 0.0) {
lla[0] = StrictMath.atan2(y,x);
} else {
if (y > 0) {
lla[0] = PI_OVER_2;
} else if (y < 0) {
lla[0] = -PI_OVER_2;
} else {
atPole = true;
lla[0] = 0.0D;
if (z > 0.0) {
lla[1] = PI_OVER_2;
} else if (z < 0.0) {
lla[1] = -PI_OVER_2;
} else {
lla[1] = PI_OVER_2;
lla[2] = -SEMIMINOR_AXIS;
return lla;
}
}
}
final double w2 = x*x + y*y;
final double w = StrictMath.sqrt(w2);
final double t0 = z * ad_c;
final double s0 = StrictMath.sqrt(t0 * t0 + w2);
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 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);
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);
}
if (!atPole) {
lla[1] = StrictMath.atan(sinP1/cosP1);
}
lla[0] = StrictMath.toDegrees(lla[0]);
lla[1] = StrictMath.toDegrees(lla[1]);
return lla;
}
/**
* Converts from geodesic lon lat alt to geocentric earth-centered earth-fixed
* @param lon geodesic longitude
* @param lat geodesic latitude
* @param alt geodesic altitude
* @param ecf reusable earth-centered earth-fixed result
* @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);
final double sl = StrictMath.sin(lat);
final double s2 = sl*sl;
final double cl = StrictMath.cos(lat);
final double ge2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMAJOR_AXIS2);
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;
}
assert (lat >= -PI_OVER_2) || (lat <= PI_OVER_2);
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;
return ecf;
}
/**
* Converts from lat lon alt (in degrees) to East North Up right-hand coordinate system
* @param lon longitude in degrees
* @param lat latitude in degrees
* @param alt altitude in meters
* @param centerLon reference point longitude in degrees
* @param centerLat reference point latitude in degrees
* @param centerAlt reference point altitude in meters
* @param enu result east, north, up coordinate
* @return east, north, up coordinate
*/
public static double[] llaToENU(final double lon, final double lat, final double alt, double centerLon,
double centerLat, final double centerAlt, double[] enu) {
if (enu == null) {
enu = new double[3];
}
// convert point to ecf coordinates
final double[] ecf = llaToECF(lon, lat, alt, null);
// convert from ecf to enu
return ecfToENU(ecf[0], ecf[1], ecf[2], centerLon, centerLat, centerAlt, enu);
}
/**
* Converts from East North Up right-hand rule to lat lon alt in degrees
* @param x easting (in meters)
* @param y northing (in meters)
* @param z up (in meters)
* @param centerLon reference point longitude (in degrees)
* @param centerLat reference point latitude (in degrees)
* @param centerAlt reference point altitude (in meters)
* @param lla resulting lat, lon, alt point (in degrees)
* @return lat, lon, alt point (in degrees)
*/
public static double[] enuToLLA(final double x, final double y, final double z, final double centerLon,
final double centerLat, final double centerAlt, double[] lla) {
// convert enuToECF
if (lla == null) {
lla = new double[3];
}
// convert enuToECF, storing intermediate result in lla
lla = enuToECF(x, y, z, centerLon, centerLat, centerAlt, lla);
// convert ecf to LLA
return ecfToLLA(lla[0], lla[1], lla[2], lla);
}
/**
* Convert from Earth-Centered-Fixed to Easting, Northing, Up Right Hand System
* @param x ECF X coordinate (in meters)
* @param y ECF Y coordinate (in meters)
* @param z ECF Z coordinate (in meters)
* @param centerLon ENU origin longitude (in degrees)
* @param centerLat ENU origin latitude (in degrees)
* @param centerAlt ENU altitude (in meters)
* @param enu reusable enu result
* @return Easting, Northing, Up coordinate
*/
public static double[] ecfToENU(double x, double y, double z, final double centerLon,
final double centerLat, final double centerAlt, double[] enu) {
if (enu == null) {
enu = new double[3];
}
// create rotation matrix and rotate to enu orientation
final double[][] phi = createPhiTransform(centerLon, centerLat, null);
// convert origin to ENU
final double[] originECF = llaToECF(centerLon, centerLat, centerAlt, null);
final double[] originENU = new double[3];
originENU[0] = ((phi[0][0] * originECF[0]) + (phi[0][1] * originECF[1]) + (phi[0][2] * originECF[2]));
originENU[1] = ((phi[1][0] * originECF[0]) + (phi[1][1] * originECF[1]) + (phi[1][2] * originECF[2]));
originENU[2] = ((phi[2][0] * originECF[0]) + (phi[2][1] * originECF[1]) + (phi[2][2] * originECF[2]));
// rotate then translate
enu[0] = ((phi[0][0] * x) + (phi[0][1] * y) + (phi[0][2] * z)) - originENU[0];
enu[1] = ((phi[1][0] * x) + (phi[1][1] * y) + (phi[1][2] * z)) - originENU[1];
enu[2] = ((phi[2][0] * x) + (phi[2][1] * y) + (phi[2][2] * z)) - originENU[2];
return enu;
}
/**
* Convert from Easting, Northing, Up Right-Handed system to Earth Centered Fixed system
* @param x ENU x coordinate (in meters)
* @param y ENU y coordinate (in meters)
* @param z ENU z coordinate (in meters)
* @param centerLon ENU origin longitude (in degrees)
* @param centerLat ENU origin latitude (in degrees)
* @param centerAlt ENU origin altitude (in meters)
* @param ecf reusable ecf result
* @return ecf result coordinate
*/
public static double[] enuToECF(final double x, final double y, final double z, double centerLon,
double centerLat, final double centerAlt, double[] ecf) {
if (ecf == null) {
ecf = new double[3];
}
double[][] phi = createTransposedPhiTransform(centerLon, centerLat, null);
double[] ecfOrigin = llaToECF(centerLon, centerLat, centerAlt, null);
// rotate and translate
ecf[0] = (phi[0][0]*x + phi[0][1]*y + phi[0][2]*z) + ecfOrigin[0];
ecf[1] = (phi[1][0]*x + phi[1][1]*y + phi[1][2]*z) + ecfOrigin[1];
ecf[2] = (phi[2][0]*x + phi[2][1]*y + phi[2][2]*z) + ecfOrigin[2];
return ecf;
}
/**
* Create the rotation matrix for converting Earth Centered Fixed to Easting Northing Up
* @param originLon ENU origin longitude (in degrees)
* @param originLat ENU origin latitude (in degrees)
* @param phiMatrix reusable phi matrix result
* @return phi rotation matrix
*/
private static double[][] createPhiTransform(double originLon, double originLat, double[][] phiMatrix) {
if (phiMatrix == null) {
phiMatrix = new double[3][3];
}
originLon = StrictMath.toRadians(originLon);
originLat = StrictMath.toRadians(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);
phiMatrix[0][0] = -sLon;
phiMatrix[0][1] = cLon;
phiMatrix[0][2] = 0.0D;
phiMatrix[1][0] = -sLat * cLon;
phiMatrix[1][1] = -sLat * sLon;
phiMatrix[1][2] = cLat;
phiMatrix[2][0] = cLat * cLon;
phiMatrix[2][1] = cLat * sLon;
phiMatrix[2][2] = sLat;
return phiMatrix;
}
/**
* Create the transposed rotation matrix for converting Easting Northing Up coordinates to Earth Centered Fixed
* @param originLon ENU origin longitude (in degrees)
* @param originLat ENU origin latitude (in degrees)
* @param phiMatrix reusable phi rotation matrix result
* @return transposed phi rotation matrix
*/
private static double[][] createTransposedPhiTransform(double originLon, double originLat, double[][] phiMatrix) {
if (phiMatrix == null) {
phiMatrix = new double[3][3];
}
originLon = StrictMath.toRadians(originLon);
originLat = StrictMath.toRadians(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);
phiMatrix[0][0] = -sLon;
phiMatrix[1][0] = cLon;
phiMatrix[2][0] = 0.0D;
phiMatrix[0][1] = -sLat * cLon;
phiMatrix[1][1] = -sLat * sLon;
phiMatrix[2][1] = cLat;
phiMatrix[0][2] = cLat * cLon;
phiMatrix[1][2] = cLat * sLon;
phiMatrix[2][2] = sLat;
return phiMatrix;
}
/**
* Finds a point along a bearing from a given lon,lat geolocation using vincenty's distance formula
*
* @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[] pointFromLonLatBearing(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 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 A = 1 + uSq/16384D*(4096D + uSq * (-768D + uSq * (320D - 175D*uSq)));
final double B = uSq/1024D * (256D + uSq * (-128D + uSq * (74D - 47D * uSq)));
double sigma = dist / (SEMIMINOR_AXIS*A);
double sigmaP;
double sinSigma, cosSigma, cos2SigmaM, deltaSigma;
do {
cos2SigmaM = StrictMath.cos(2*sig1 + sigma);
sinSigma = StrictMath.sin(sigma);
cosSigma = StrictMath.cos(sigma);
deltaSigma = B * sinSigma * (cos2SigmaM + (B/4D) * (cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
(B/6) * cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
sigmaP = sigma;
sigma = dist / (SEMIMINOR_AXIS*A) + deltaSigma;
} while (StrictMath.abs(sigma-sigmaP) > 1E-12);
final double tmp = sinU1*sinSigma - cosU1*cosSigma*cosA1;
final double lat2 = StrictMath.atan2(sinU1*cosSigma + cosU1*sinSigma*cosA1,
(1-FLATTENING) * StrictMath.sqrt(sinAlpha*sinAlpha + tmp*tmp));
final double lambda = StrictMath.atan2(sinSigma*sinA1, cosU1*cosSigma - sinU1*sinSigma*cosA1);
final double c = FLATTENING/16 * cosSqAlpha * (4 + FLATTENING * (4 - 3 * cosSqAlpha));
final double lam = lambda - (1-c) * FLATTENING * sinAlpha *
(sigma + c * sinSigma * (cos2SigmaM + c * cosSigma * (-1 + 2* cos2SigmaM*cos2SigmaM)));
pt[0] = lon + StrictMath.toDegrees(lam);
pt[1] = StrictMath.toDegrees(lat2);
return pt;
}
}

View File

@ -17,21 +17,14 @@ package org.apache.lucene.util;
* limitations under the License.
*/
import java.util.ArrayList;
/**
* Basic reusable geo-spatial utility methods
*
* @lucene.experimental
*/
public final class GeoUtils {
// WGS84 earth-ellipsoid major (a) minor (b) radius, (f) flattening and eccentricity (e)
private static final double SEMIMAJOR_AXIS = 6_378_137; // [m]
private static final double FLATTENING = 1.0/298.257223563;
private static final double SEMIMINOR_AXIS = SEMIMAJOR_AXIS * (1.0 - FLATTENING); //6_356_752.31420; // [m]
private static final double ECCENTRICITY = StrictMath.sqrt((2.0 - FLATTENING) * FLATTENING);
private static final double PI_OVER_2 = StrictMath.PI / 2.0D;
private static final double SEMIMAJOR_AXIS2 = SEMIMAJOR_AXIS * SEMIMINOR_AXIS;
private static final double SEMIMINOR_AXIS2 = SEMIMINOR_AXIS * SEMIMINOR_AXIS;
private static final short MIN_LON = -180;
private static final short MIN_LAT = -90;
public static final short BITS = 31;
@ -55,15 +48,15 @@ public final class GeoUtils {
private GeoUtils() {
}
public static final Long mortonHash(final double lon, final double lat) {
public static Long mortonHash(final double lon, final double lat) {
return BitUtil.interleave(scaleLon(lon), scaleLat(lat));
}
public static final double mortonUnhashLon(final long hash) {
public static double mortonUnhashLon(final long hash) {
return unscaleLon(BitUtil.deinterleave(hash));
}
public static final double mortonUnhashLat(final long hash) {
public static double mortonUnhashLat(final long hash) {
return unscaleLat(BitUtil.deinterleave(hash >>> 1));
}
@ -83,126 +76,45 @@ public final class GeoUtils {
return (val / LAT_SCALE) + MIN_LAT;
}
public static final double compare(final double v1, final double v2) {
public static double compare(final double v1, final double v2) {
final double compare = v1-v2;
return Math.abs(compare) <= TOLERANCE ? 0 : compare;
}
/**
* Puts longitude in range of -180 to +180.
*/
public static double normalizeLon(double lon_deg) {
if (lon_deg >= -180 && lon_deg <= 180) {
return lon_deg; //common case, and avoids slight double precision shifting
}
double off = (lon_deg + 180) % 360;
if (off < 0) {
return 180 + off;
} else if (off == 0 && lon_deg > 0) {
return 180;
} else {
return -180 + off;
}
}
/**
* Puts latitude in range of -90 to 90.
*/
public static double normalizeLat(double lat_deg) {
if (lat_deg >= -90 && lat_deg <= 90) {
return lat_deg; //common case, and avoids slight double precision shifting
}
double off = Math.abs((lat_deg + 90) % 360);
return (off <= 180 ? off : 360-off) - 90;
}
public static final 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);
}
/**
* Converts from geodesic lon lat alt to geocentric earth-centered earth-fixed
* @param lon geodesic longitude
* @param lat geodesic latitude
* @param alt geodesic altitude
* @param ecf reusable earth-centered earth-fixed result
* @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);
final double sl = StrictMath.sin(lat);
final double s2 = sl*sl;
final double cl = StrictMath.cos(lat);
final double ge2 = (SEMIMAJOR_AXIS2 - SEMIMINOR_AXIS2)/(SEMIMAJOR_AXIS2);
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;
}
assert ((lat >= -PI_OVER_2) || (lat <= PI_OVER_2));
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;
return ecf;
}
/**
* Converts from geocentric earth-centered earth-fixed to geodesic lat/lon/alt
* @param x Cartesian x coordinate
* @param y Cartesian y coordinate
* @param z Cartesian z coordinate
* @param lla 0: longitude 1: latitude: 2: altitude
* @return double array as 0: longitude 1: latitude 2: altitude
*/
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)
lla = new double[3];
if (x != 0.0) {
lla[0] = StrictMath.atan2(y,x);
} else {
if (y > 0) {
lla[0] = PI_OVER_2;
} else if (y < 0) {
lla[0] = -PI_OVER_2;
} else {
atPole = true;
lla[0] = 0.0D;
if (z > 0.0) {
lla[1] = PI_OVER_2;
} else if (z < 0.0) {
lla[1] = -PI_OVER_2;
} else {
lla[1] = PI_OVER_2;
lla[2] = -SEMIMINOR_AXIS;
return lla;
}
}
}
final double w2 = x*x + y*y;
final double w = StrictMath.sqrt(w2);
final double t0 = z * ad_c;
final double s0 = StrictMath.sqrt(t0 * t0 + w2);
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 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);
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);
}
if (!atPole) {
lla[1] = StrictMath.atan(sinP1/cosP1);
}
lla[0] = StrictMath.toDegrees(lla[0]);
lla[1] = StrictMath.toDegrees(lla[1]);
return lla;
}
/**
* simple even-odd point in polygon computation
* 1. Determine if point is contained in the longitudinal range
@ -236,8 +148,9 @@ public final class GeoUtils {
for (int i = 0; i < numberOfLeadingZeros; i++) {
s.append('0');
}
if (term != 0)
if (term != 0) {
s.append(Long.toBinaryString(term));
}
return s.toString();
}
@ -280,13 +193,14 @@ public final class GeoUtils {
/**
* Computes whether a rectangle crosses a shape. (touching not allowed)
*/
public static final 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) {
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))
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;
@ -320,14 +234,49 @@ public final class GeoUtils {
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))
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
*
* @param lon longitudinal center of circle (in degrees)
* @param lat latitudinal center of circle (in degrees)
* @param radius distance radius of circle (in meters)
* @return a list of lon/lat points representing the circle
*/
@SuppressWarnings({"unchecked","rawtypes"})
public static ArrayList<double[]> circleToPoly(final double lon, final double lat, final double radius) {
double angle;
// a little under-sampling (to limit the number of polygonal points): using archimedes estimation of pi
final int sides = 25;
ArrayList<double[]> geometry = new ArrayList();
double[] lons = new double[sides];
double[] lats = new double[sides];
double[] pt = new double[2];
final int sidesLen = sides-1;
for (int i=0; i<sidesLen; ++i) {
angle = (i*360/sides);
pt = GeoProjectionUtils.pointFromLonLatBearing(lon, lat, angle, radius, pt);
lons[i] = pt[0];
lats[i] = pt[1];
}
// close the poly
lons[sidesLen] = lons[0];
lats[sidesLen] = lats[0];
geometry.add(lons);
geometry.add(lats);
return geometry;
}
/**
* Computes whether a rectangle is within a given polygon (shared boundaries allowed)
*/
@ -341,6 +290,83 @@ public final class GeoUtils {
!pointInPolygon(shapeX, shapeY, rMaxY, rMaxX) || !pointInPolygon(shapeX, shapeY, rMaxY, rMinX));
}
private static boolean rectAnyCornersInCirlce( final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
final double centerLon, final double centerLat, final double radius) {
return (SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 > radius
|| SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 > radius
|| SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 > radius
|| SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 > radius);
}
public static boolean rectWithinCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
final double centerLon, final double centerLat, final double radius) {
return !(rectAnyCornersInCirlce(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radius));
}
/**
* Computes whether a rectangle crosses a circle
*/
public static boolean rectCrossesCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
final double centerLon, final double centerLat, final double radius) {
return rectAnyCornersInCirlce(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radius)
|| lineCrossesSphere(rMinX, rMinY, 0, rMaxX, rMinY, 0, centerLon, centerLat, 0, radius)
|| lineCrossesSphere(rMaxX, rMinY, 0, rMaxX, rMaxY, 0, centerLon, centerLat, 0, radius)
|| lineCrossesSphere(rMaxX, rMaxY, 0, rMinX, rMaxY, 0, centerLon, centerLat, 0, radius)
|| lineCrossesSphere(rMinX, rMaxY, 0, rMinX, rMinY, 0, centerLon, centerLat, 0, radius);
}
/**
* 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 radius 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 radius) {
// 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) - (radius*radius);
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;
}
public static boolean isValidLat(double lat) {
return Double.isNaN(lat) == false && lat >= MIN_LAT_INCL && lat <= MAX_LAT_INCL;
}

View File

@ -44,6 +44,7 @@ import org.apache.lucene.util.FixedBitSet;
import org.apache.lucene.util.GeoUtils;
import org.apache.lucene.util.IOUtils;
import org.apache.lucene.util.LuceneTestCase;
import org.apache.lucene.util.SloppyMath;
import org.apache.lucene.util.TestUtil;
import org.junit.AfterClass;
import org.junit.BeforeClass;
@ -61,10 +62,15 @@ public class TestGeoPointQuery extends LuceneTestCase {
private static final String FIELD_NAME = "geoField";
// error threshold for point-distance queries (in meters)
private static final int DISTANCE_ERR = 700;
// Global bounding box we will "cover" in the random test; we have to make this "smallish" else the queries take very long:
private static double originLat;
private static double originLon;
private static double range;
// private static double range;
private static double lonRange;
private static double latRange;
@BeforeClass
public static void beforeClass() throws Exception {
@ -74,11 +80,21 @@ public class TestGeoPointQuery extends LuceneTestCase {
// number of ranges that can be created in degenerate cases.
// Between 1.0 and 3.0:
range = 2*(random().nextDouble() + 0.5);
originLon = GeoUtils.MIN_LON_INCL + range + (GeoUtils.MAX_LON_INCL - GeoUtils.MIN_LON_INCL - 2*range) * random().nextDouble();
originLat = GeoUtils.MIN_LAT_INCL + range + (GeoUtils.MAX_LAT_INCL - GeoUtils.MIN_LAT_INCL - 2*range) * random().nextDouble();
// range = 2*(random().nextDouble() + 0.5);
// Between 1.0 and 90.0
//lonRange = 1.0 + (90.0 - 1.0) * random().nextDouble();
//latRange = 1.0 + (45.0 - 1.0) * random().nextDouble();
// Between 1.0 and 3.0:
lonRange = 2*(random().nextDouble() + 0.5);
latRange = 2*(random().nextDouble() + 0.5);
originLon = GeoUtils.MIN_LON_INCL + lonRange + (GeoUtils.MAX_LON_INCL - GeoUtils.MIN_LON_INCL - 2*lonRange) * random().nextDouble();
originLon = GeoUtils.normalizeLon(originLon);
originLat = GeoUtils.MIN_LAT_INCL + latRange + (GeoUtils.MAX_LAT_INCL - GeoUtils.MIN_LAT_INCL - 2*latRange) * random().nextDouble();
originLat = GeoUtils.normalizeLat(originLat);
if (VERBOSE) {
System.out.println("TEST: originLon=" + originLon + " originLat=" + originLat + " range=" + range);
System.out.println("TEST: originLon=" + originLon + " lonRange= " + lonRange + " originLat=" + originLat + " latRange=" + latRange);
}
RandomIndexWriter writer = new RandomIndexWriter(random(), directory,
newIndexWriterConfig(new MockAnalyzer(random()))
@ -99,7 +115,10 @@ public class TestGeoPointQuery extends LuceneTestCase {
new GeoPointField(FIELD_NAME, -83.99724648980559, 58.29438379542874, storedPoint),
new GeoPointField(FIELD_NAME, -26.779373834241003, 33.541429799076354, storedPoint),
new GeoPointField(FIELD_NAME, -77.35379276106497, 26.774024500421728, storedPoint),
new GeoPointField(FIELD_NAME, -14.796283808944777, -62.455081198245665, storedPoint)};
new GeoPointField(FIELD_NAME, -14.796283808944777, -62.455081198245665, storedPoint),
new GeoPointField(FIELD_NAME, -178.8538113027811, 32.94823588839368, storedPoint),
new GeoPointField(FIELD_NAME, 178.8538113027811, 32.94823588839368, storedPoint),
new GeoPointField(FIELD_NAME, -179.5, -44.5, storedPoint)};
for (GeoPointField p : pts) {
Document doc = new Document();
@ -130,6 +149,11 @@ public class TestGeoPointQuery extends LuceneTestCase {
return searcher.search(q, limit);
}
private TopDocs geoDistanceQuery(double lon, double lat, double radius, int limit) throws Exception {
GeoPointDistanceQuery q = new GeoPointDistanceQuery(FIELD_NAME, lon, lat, radius);
return searcher.search(q, limit);
}
@Test
public void testBBoxQuery() throws Exception {
TopDocs td = bboxQuery(-96.7772, 32.778650, -96.77690000, 32.778950, 5);
@ -138,11 +162,11 @@ public class TestGeoPointQuery extends LuceneTestCase {
@Test
public void testPolyQuery() throws Exception {
TopDocs td = polygonQuery( new double[] {-96.7682647, -96.8280029, -96.6288757, -96.4929199,
-96.6041564, -96.7449188, -96.76826477, -96.7682647},
new double[] { 33.073130, 32.9942669, 32.938386, 33.0374494,
33.1369762, 33.1162747, 33.073130, 33.073130}, 5);
assertEquals("GeoPolygonQuery failed", td.totalHits, 1);
TopDocs td = polygonQuery(new double[]{-96.7682647, -96.8280029, -96.6288757, -96.4929199,
-96.6041564, -96.7449188, -96.76826477, -96.7682647},
new double[]{33.073130, 32.9942669, 32.938386, 33.0374494,
33.1369762, 33.1162747, 33.073130, 33.073130}, 5);
assertEquals("GeoPolygonQuery failed", 1, td.totalHits);
}
@Test
@ -169,6 +193,50 @@ public class TestGeoPointQuery extends LuceneTestCase {
assertTrue(GeoUtils.rectWithinPoly(-5, 0, -2, 5, px, py, xMin, yMin, xMax, yMax));
}
@Test
public void testBBoxCrossDateline() throws Exception {
TopDocs td = bboxQuery(179.0, -45.0, -179.0, -44.0, 20);
assertEquals("BBoxCrossDateline query failed", 1, td.totalHits);
}
@Test
public void testWholeMap() throws Exception {
TopDocs td = bboxQuery(-179.9, -89.9, 179.9, 89.9, 20);
assertEquals("testWholeMap failed", 14, td.totalHits);
}
@Test
public void testInvalidBBox() throws Exception {
try {
bboxQuery(179.0, -92.0, 181.0, -91.0, 20);
} catch(Exception e) {
return;
}
throw new Exception("GeoBoundingBox should not accept invalid lat/lon");
}
@Test
public void testGeoDistanceQuery() throws Exception {
TopDocs td = geoDistanceQuery(-96.4538113027811, 32.94823588839368, 600000, 20);
assertEquals("GeoDistanceQuery failed", 6, td.totalHits);
}
@Test
public void testGeoDistanceQueryCrossDateline() throws Exception {
TopDocs td = geoDistanceQuery(-179.9538113027811, 32.94823588839368, 120000, 20);
assertEquals("GeoDistanceQuery failed", 2, td.totalHits);
}
@Test
public void testInvalidGeoDistanceQuery() throws Exception {
try {
geoDistanceQuery(181.0, 92.0, 120000, 20);
} catch (Exception e) {
return;
}
throw new Exception("GeoDistanceQuery should not accept invalid lat/lon as origin");
}
public void testRandomTiny() throws Exception {
// Make sure single-leaf-node case is OK:
doTestRandom(10);
@ -298,7 +366,7 @@ public class TestGeoPointQuery extends LuceneTestCase {
int numThreads = TestUtil.nextInt(random(), 2, 5);
List<Thread> threads = new ArrayList<>();
final int iters = atLeast(100);
final int iters = atLeast(10);
final CountDownLatch startingGun = new CountDownLatch(1);
@ -319,112 +387,102 @@ public class TestGeoPointQuery extends LuceneTestCase {
NumericDocValues docIDToID = MultiDocValues.getNumericValues(r, "id");
for (int iter=0;iter<iters;iter++) {
double lat0 = randomLat();
double lat1 = randomLat();
double lon0 = randomLon();
double lon1 = randomLon();
if (lat1 < lat0) {
double x = lat0;
lat0 = lat1;
lat1 = x;
}
if (lon1 < lon0) {
double x = lon0;
lon0 = lon1;
lon1 = x;
}
if (VERBOSE) {
System.out.println("\nTEST: iter=" + iter + " lat=" + lat0 + " TO " + lat1 + " lon=" + lon0 + " TO " + lon1);
System.out.println("\nTEST: iter=" + iter);
}
Query query;
boolean tooBigBBox = false;
boolean polySearch = false;
double bboxLat0 = lat0;
double bboxLat1 = lat1;
double bboxLon0 = lon0;
double bboxLon1 = lon1;
VerifyHits verifyHits;
if (random().nextBoolean()) {
query = new GeoPointInBBoxQuery(FIELD_NAME, lon0, lat0, lon1, lat1);
} else {
polySearch = true;
if (random().nextBoolean()) {
// Intentionally pass a "too big" bounding box:
double pct = random().nextDouble()*0.5;
double width = lon1-lon0;
bboxLon0 = Math.max(-180.0, lon0-width*pct);
bboxLon1 = Math.min(180.0, lon1+width*pct);
double height = lat1-lat0;
bboxLat0 = Math.max(-90.0, lat0-height*pct);
bboxLat1 = Math.min(90.0, lat1+height*pct);
tooBigBBox = true;
GeoBoundingBox bbox = randomBBox();
query = new GeoPointInBBoxQuery(FIELD_NAME, bbox.minLon, bbox.minLat, bbox.maxLon, bbox.maxLat);
verifyHits = new VerifyHits() {
@Override
protected Boolean shouldMatch(double pointLat, double pointLon) {
// morton encode & decode to compare apples to apples (that is, compare with same hash precision error
// present in the index)
long pointHash = GeoUtils.mortonHash(pointLon, pointLat);
pointLon = GeoUtils.mortonUnhashLon(pointHash);
pointLat = GeoUtils.mortonUnhashLat(pointHash);
if (bboxQueryCanBeWrong(bbox, pointLat, pointLon)) {
return null;
} else {
return rectContainsPointEnc(bbox, pointLat, pointLon);
}
}
};
} else if (random().nextBoolean()) {
// generate a random bounding box
GeoBoundingBox bbox = randomBBox();
double centerLat = bbox.minLat + ((bbox.maxLat - bbox.minLat)/2.0);
double centerLon = bbox.minLon + ((bbox.maxLon - bbox.minLon)/2.0);
// radius (in meters) as a function of the random generated bbox
// TODO: change 100 back to 1000
final double radius = SloppyMath.haversin(centerLat, centerLon, bbox.minLat, centerLon)*100;
if (VERBOSE) {
System.out.println("\t radius = " + radius);
}
// query using the centroid of the bounding box
query = new GeoPointDistanceQuery(FIELD_NAME, centerLon, centerLat, radius);
verifyHits = new VerifyHits() {
@Override
protected Boolean shouldMatch(double pointLat, double pointLon) {
if (Double.isNaN(pointLat) || Double.isNaN(pointLon)) {
return null;
}
if (radiusQueryCanBeWrong(centerLat, centerLon, pointLon, pointLat, radius)) {
return null;
} else {
return distanceContainsPt(centerLon, centerLat, pointLon, pointLat, radius);
}
}
};
} else {
GeoBoundingBox bbox = randomBBox();
double[] pLats = new double[5];
double[] pLons = new double[5];
pLats[0] = bboxLat0;
pLons[0] = bboxLon0;
pLats[1] = bboxLat1;
pLons[1] = bboxLon0;
pLats[2] = bboxLat1;
pLons[2] = bboxLon1;
pLats[3] = bboxLat0;
pLons[3] = bboxLon1;
pLats[4] = bboxLat0;
pLons[4] = bboxLon0;
query = new GeoPointInPolygonQuery(FIELD_NAME, bboxLon0, bboxLat0, bboxLon1, bboxLat1, pLons, pLats);
}
pLats[0] = bbox.minLat;
pLons[0] = bbox.minLon;
pLats[1] = bbox.maxLat;
pLons[1] = bbox.minLon;
pLats[2] = bbox.maxLat;
pLons[2] = bbox.maxLon;
pLats[3] = bbox.minLat;
pLons[3] = bbox.maxLon;
pLats[4] = bbox.minLat;
pLons[4] = bbox.minLon;
query = new GeoPointInPolygonQuery(FIELD_NAME, pLons, pLats);
final FixedBitSet hits = new FixedBitSet(r.maxDoc());
s.search(query, new SimpleCollector() {
verifyHits = new VerifyHits() {
@Override
protected Boolean shouldMatch(double pointLat, double pointLon) {
// morton encode & decode to compare apples to apples (that is, compare with same hash precision error
// present in the index)
long pointHash = GeoUtils.mortonHash(pointLon, pointLat);
pointLon = GeoUtils.mortonUnhashLon(pointHash);
pointLat = GeoUtils.mortonUnhashLat(pointHash);
private int docBase;
@Override
public boolean needsScores() {
return false;
}
@Override
protected void doSetNextReader(LeafReaderContext context) throws IOException {
docBase = context.docBase;
}
@Override
public void collect(int doc) {
hits.set(docBase+doc);
}
});
for(int docID=0;docID<r.maxDoc();docID++) {
int id = (int) docIDToID.get(docID);
if (polySearch) {
lat0 = bboxLat0;
lon0 = bboxLon0;
lat1 = bboxLat1;
lon1 = bboxLon1;
}
// morton encode & decode to compare apples to apples (that is, compare with same hash precision error
// present in the index)
final long pointHash = GeoUtils.mortonHash(lons[id], lats[id]);
final double pointLon = GeoUtils.mortonUnhashLon(pointHash);
final double pointLat = GeoUtils.mortonUnhashLat(pointHash);
if (!tolerateIgnorance(lat0, lat1, lon0, lon1, pointLat, pointLon)) {
boolean expected = (deleted.contains(id) == false) &&
rectContainsPointEnc(lat0, lat1, lon0, lon1, pointLat, pointLon);
if (hits.get(docID) != expected) {
System.out.println(Thread.currentThread().getName() + ": iter=" + iter + " id=" + id + " docID=" + docID + " lat=" + pointLat + " lon=" + pointLon + " (bbox: lat=" + lat0 + " TO " + lat1 + " lon=" + lon0 + " TO " + lon1 + ") expected " + expected + " but got: " + hits.get(docID) + " deleted?=" + deleted.contains(id) + " query=" + query);
if (tooBigBBox) {
System.out.println(" passed too-big bbox: lat=" + bboxLat0 + " TO " + bboxLat1 + " lon=" + bboxLon0 + " TO " + bboxLon1);
if (bboxQueryCanBeWrong(bbox, pointLat, pointLon)) {
return null;
} else {
return rectContainsPointEnc(bbox, pointLat, pointLon);
}
}
fail("wrong result");
}
}
};
}
verifyHits.test(s, docIDToID, deleted, query, lats, lons);
}
}
};
@ -441,35 +499,132 @@ public class TestGeoPointQuery extends LuceneTestCase {
IOUtils.close(r, dir);
}
private static boolean rectContainsPointEnc(double rectLatMin, double rectLatMax,
double rectLonMin, double rectLonMax,
double pointLat, double pointLon) {
if (Double.isNaN(pointLat)) {
return false;
private static abstract class VerifyHits {
public void test(IndexSearcher s, NumericDocValues docIDToID, Set<Integer> deleted, Query query, double[] lats, double[] lons) throws Exception {
int maxDoc = s.getIndexReader().maxDoc();
final FixedBitSet hits = new FixedBitSet(maxDoc);
s.search(query, new SimpleCollector() {
private int docBase;
@Override
public boolean needsScores() {
return false;
}
@Override
protected void doSetNextReader(LeafReaderContext context) throws IOException {
docBase = context.docBase;
}
@Override
public void collect(int doc) {
hits.set(docBase+doc);
}
});
for(int docID=0;docID<maxDoc;docID++) {
int id = (int) docIDToID.get(docID);
Boolean expected;
if (deleted.contains(id)) {
expected = false;
} else {
expected = shouldMatch(lats[id], lons[id]);
}
// null means it's a borderline case which is allowed to be wrong:
if (expected != null) {
if (hits.get(docID) != expected) {
System.out.println(Thread.currentThread().getName() + ": id=" + id +
" docID=" + docID + " lat=" + lats[id] + " lon=" + lons[id] +
" deleted?=" + deleted.contains(id) + " expected=" + expected + " but got " + hits.get(docID) +
" query=" + query);
fail("wrong hit");
}
}
}
}
return GeoUtils.bboxContains(pointLon, pointLat, rectLonMin, rectLatMin, rectLonMax, rectLatMax);
/** Return true if we definitely should match, false if we definitely
* should not match, and null if it's a borderline case which might
* go either way. */
protected abstract Boolean shouldMatch(double lat, double lon);
}
private static boolean tolerateIgnorance(final double minLat, final double maxLat,
final double minLon, final double maxLon,
final double lat, final double lon) {
private static boolean distanceContainsPt(double lonA, double latA, double lonB, double latB, final double radius) {
final long hashedPtA = GeoUtils.mortonHash(lonA, latA);
lonA = GeoUtils.mortonUnhashLon(hashedPtA);
latA = GeoUtils.mortonUnhashLat(hashedPtA);
final long hashedPtB = GeoUtils.mortonHash(lonB, latB);
lonB = GeoUtils.mortonUnhashLon(hashedPtB);
latB = GeoUtils.mortonUnhashLat(hashedPtB);
return (SloppyMath.haversin(latA, lonA, latB, lonB)*1000.0 <= radius);
}
private static boolean rectContainsPointEnc(GeoBoundingBox bbox, double pointLat, double pointLon) {
// We should never see a deleted doc here?
assert Double.isNaN(pointLat) == false;
return GeoUtils.bboxContains(pointLon, pointLat, bbox.minLon, bbox.minLat, bbox.maxLon, bbox.maxLat);
}
private static boolean radiusQueryCanBeWrong(double centerLat, double centerLon, double ptLon, double ptLat,
final double radius) {
final long hashedCntr = GeoUtils.mortonHash(centerLon, centerLat);
centerLon = GeoUtils.mortonUnhashLon(hashedCntr);
centerLat = GeoUtils.mortonUnhashLat(hashedCntr);
final long hashedPt = GeoUtils.mortonHash(ptLon, ptLat);
ptLon = GeoUtils.mortonUnhashLon(hashedPt);
ptLat = GeoUtils.mortonUnhashLat(hashedPt);
double ptDistance = SloppyMath.haversin(centerLat, centerLon, ptLat, ptLon)*1000.0;
double delta = StrictMath.abs(ptDistance - radius);
// if its within the distance error then it can be wrong
return delta < DISTANCE_ERR;
}
private static boolean bboxQueryCanBeWrong(GeoBoundingBox bbox, double lat, double lon) {
// we can tolerate variance at the GeoUtils.TOLERANCE decimal place
final int tLon = (int)(lon/(GeoUtils.TOLERANCE-1));
final int tLat = (int)(lat/(GeoUtils.TOLERANCE-1));
final int tMinLon = (int)(minLon/(GeoUtils.TOLERANCE-1));
final int tMinLat = (int)(minLat/(GeoUtils.TOLERANCE-1));
final int tMaxLon = (int)(maxLon/(GeoUtils.TOLERANCE-1));
final int tMaxLat = (int)(maxLat/(GeoUtils.TOLERANCE-1));
final int tMinLon = (int)(bbox.minLon/(GeoUtils.TOLERANCE-1));
final int tMinLat = (int)(bbox.minLat/(GeoUtils.TOLERANCE-1));
final int tMaxLon = (int)(bbox.maxLon/(GeoUtils.TOLERANCE-1));
final int tMaxLat = (int)(bbox.maxLat/(GeoUtils.TOLERANCE-1));
return ((tMinLon - tLon) == 0 || (tMinLat - tLat) == 0
|| (tMaxLon - tLon) == 0 || (tMaxLat - tLat) == 0);
}
private static double randomLat() {
return originLat + range * (random().nextDouble()-0.5);
return GeoUtils.normalizeLat(originLat + latRange * (random().nextDouble() - 0.5));
}
private static double randomLon() {
return originLon + range * (random().nextDouble()-0.5);
return GeoUtils.normalizeLon(originLon + lonRange * (random().nextDouble() - 0.5));
}
private static GeoBoundingBox randomBBox() {
double lat0 = randomLat();
double lat1 = randomLat();
double lon0 = randomLon();
double lon1 = randomLon();
if (lat1 < lat0) {
double x = lat0;
lat0 = lat1;
lat1 = x;
}
if (lon1 < lon0) {
double x = lon0;
lon0 = lon1;
lon1 = x;
}
return new GeoBoundingBox(lon0, lon1, lat0, lat1);
}
}