Merge pull request #585 from mschoch/fix-geo-point-dist
fix geo point distance search
This commit is contained in:
commit
11a45d6f9c
139
geo/geo.go
139
geo/geo.go
@ -15,6 +15,7 @@
|
|||||||
package geo
|
package geo
|
||||||
|
|
||||||
import (
|
import (
|
||||||
|
"fmt"
|
||||||
"math"
|
"math"
|
||||||
|
|
||||||
"github.com/blevesearch/bleve/numeric"
|
"github.com/blevesearch/bleve/numeric"
|
||||||
@ -26,6 +27,12 @@ var GeoBits uint = 32
|
|||||||
|
|
||||||
var minLon = -180.0
|
var minLon = -180.0
|
||||||
var minLat = -90.0
|
var minLat = -90.0
|
||||||
|
var maxLon = 180.0
|
||||||
|
var maxLat = 90.0
|
||||||
|
var minLonRad = minLon * degreesToRadian
|
||||||
|
var minLatRad = minLat * degreesToRadian
|
||||||
|
var maxLonRad = maxLon * degreesToRadian
|
||||||
|
var maxLatRad = maxLat * degreesToRadian
|
||||||
var geoTolerance = 1E-6
|
var geoTolerance = 1E-6
|
||||||
var lonScale = float64((uint64(0x1)<<GeoBits)-1) / 360.0
|
var lonScale = float64((uint64(0x1)<<GeoBits)-1) / 360.0
|
||||||
var latScale = float64((uint64(0x1)<<GeoBits)-1) / 180.0
|
var latScale = float64((uint64(0x1)<<GeoBits)-1) / 180.0
|
||||||
@ -91,26 +98,8 @@ func BoundingBoxContains(lon, lat, minLon, minLat, maxLon, maxLat float64) bool
|
|||||||
compareGeo(lat, minLat) >= 0 && compareGeo(lat, maxLat) <= 0
|
compareGeo(lat, minLat) >= 0 && compareGeo(lat, maxLat) <= 0
|
||||||
}
|
}
|
||||||
|
|
||||||
// ComputeBoundingBox will compute a bounding box around the provided point
|
|
||||||
// which surrounds a circle of the provided radius (in meters).
|
|
||||||
func ComputeBoundingBox(centerLon, centerLat,
|
|
||||||
radius float64) (upperLeftLon float64, upperLeftLat float64,
|
|
||||||
lowerRightLon float64, lowerRightLat float64) {
|
|
||||||
_, tlat := pointFromLonLatBearing(centerLon, centerLat, 0, radius)
|
|
||||||
rlon, _ := pointFromLonLatBearing(centerLon, centerLat, 90, radius)
|
|
||||||
_, blat := pointFromLonLatBearing(centerLon, centerLat, 180, radius)
|
|
||||||
llon, _ := pointFromLonLatBearing(centerLon, centerLat, 270, radius)
|
|
||||||
return normalizeLon(llon), normalizeLat(tlat),
|
|
||||||
normalizeLon(rlon), normalizeLat(blat)
|
|
||||||
}
|
|
||||||
|
|
||||||
const degreesToRadian = math.Pi / 180
|
const degreesToRadian = math.Pi / 180
|
||||||
const radiansToDegrees = 180 / math.Pi
|
const radiansToDegrees = 180 / math.Pi
|
||||||
const flattening = 1.0 / 298.257223563
|
|
||||||
const semiMajorAxis = 6378137
|
|
||||||
const semiMinorAxis = semiMajorAxis * (1.0 - flattening)
|
|
||||||
const semiMajorAxis2 = semiMajorAxis * semiMajorAxis
|
|
||||||
const semiMinorAxis2 = semiMinorAxis * semiMinorAxis
|
|
||||||
|
|
||||||
// DegreesToRadians converts an angle in degrees to radians
|
// DegreesToRadians converts an angle in degrees to radians
|
||||||
func DegreesToRadians(d float64) float64 {
|
func DegreesToRadians(d float64) float64 {
|
||||||
@ -122,84 +111,60 @@ func RadiansToDegrees(r float64) float64 {
|
|||||||
return r * radiansToDegrees
|
return r * radiansToDegrees
|
||||||
}
|
}
|
||||||
|
|
||||||
// pointFromLonLatBearing starts that the provide lon,lat
|
var earthMeanRadiusMeters = 6371008.7714
|
||||||
// then moves in the bearing direction (in degrees)
|
|
||||||
// this move continues for the provided distance (in meters)
|
|
||||||
// The lon, lat of this destination location is returned.
|
|
||||||
func pointFromLonLatBearing(lon, lat, bearing,
|
|
||||||
dist float64) (float64, float64) {
|
|
||||||
|
|
||||||
alpha1 := DegreesToRadians(bearing)
|
func RectFromPointDistance(lon, lat, dist float64) (float64, float64, float64, float64, error) {
|
||||||
cosA1 := math.Cos(alpha1)
|
err := checkLongitude(lon)
|
||||||
sinA1 := math.Sin(alpha1)
|
if err != nil {
|
||||||
tanU1 := (1 - flattening) * math.Tan(DegreesToRadians(lat))
|
return 0, 0, 0, 0, err
|
||||||
cosU1 := 1 / math.Sqrt(1+tanU1*tanU1)
|
}
|
||||||
sinU1 := tanU1 * cosU1
|
err = checkLatitude(lat)
|
||||||
sig1 := math.Atan2(tanU1, cosA1)
|
if err != nil {
|
||||||
sinAlpha := cosU1 * sinA1
|
return 0, 0, 0, 0, err
|
||||||
cosSqAlpha := 1 - sinAlpha*sinAlpha
|
}
|
||||||
uSq := cosSqAlpha * (semiMajorAxis2 - semiMinorAxis2) / semiMinorAxis2
|
radLon := DegreesToRadians(lon)
|
||||||
A := 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
|
radLat := DegreesToRadians(lat)
|
||||||
B := uSq / 1024 * (256 + uSq*(-128+uSq*(74-47*uSq)))
|
radDistance := (dist + 7e-2) / earthMeanRadiusMeters
|
||||||
|
|
||||||
sigma := dist / (semiMinorAxis * A)
|
minLatL := radLat - radDistance
|
||||||
|
maxLatL := radLat + radDistance
|
||||||
|
|
||||||
cos25SigmaM := math.Cos(2*sig1 + sigma)
|
var minLonL, maxLonL float64
|
||||||
sinSigma := math.Sin(sigma)
|
if minLatL > minLatRad && maxLatL < maxLatRad {
|
||||||
cosSigma := math.Cos(sigma)
|
deltaLon := asin(sin(radDistance) / cos(radLat))
|
||||||
deltaSigma := B * sinSigma * (cos25SigmaM + (B/4)*
|
minLonL = radLon - deltaLon
|
||||||
(cosSigma*(-1+2*cos25SigmaM*cos25SigmaM)-(B/6)*cos25SigmaM*
|
if minLonL < minLonRad {
|
||||||
(-1+4*sinSigma*sinSigma)*(-3+4*cos25SigmaM*cos25SigmaM)))
|
minLonL += 2 * math.Pi
|
||||||
sigmaP := sigma
|
}
|
||||||
sigma = dist/(semiMinorAxis*A) + deltaSigma
|
maxLonL = radLon + deltaLon
|
||||||
for math.Abs(sigma-sigmaP) > 1E-12 {
|
if maxLonL > maxLonRad {
|
||||||
cos25SigmaM = math.Cos(2*sig1 + sigma)
|
maxLonL -= 2 * math.Pi
|
||||||
sinSigma = math.Sin(sigma)
|
}
|
||||||
cosSigma = math.Cos(sigma)
|
} else {
|
||||||
deltaSigma = B * sinSigma * (cos25SigmaM + (B/4)*
|
// pole is inside distance
|
||||||
(cosSigma*(-1+2*cos25SigmaM*cos25SigmaM)-(B/6)*cos25SigmaM*
|
minLatL = math.Max(minLatL, minLatRad)
|
||||||
(-1+4*sinSigma*sinSigma)*(-3+4*cos25SigmaM*cos25SigmaM)))
|
maxLatL = math.Min(maxLatL, maxLatRad)
|
||||||
sigmaP = sigma
|
minLonL = minLonRad
|
||||||
sigma = dist/(semiMinorAxis*A) + deltaSigma
|
maxLonL = maxLonRad
|
||||||
}
|
}
|
||||||
|
|
||||||
tmp := sinU1*sinSigma - cosU1*cosSigma*cosA1
|
return RadiansToDegrees(minLonL),
|
||||||
lat2 := math.Atan2(sinU1*cosSigma+cosU1*sinSigma*cosA1,
|
RadiansToDegrees(maxLatL),
|
||||||
(1-flattening)*math.Sqrt(sinAlpha*sinAlpha+tmp*tmp))
|
RadiansToDegrees(maxLonL),
|
||||||
lamda := math.Atan2(sinSigma*sinA1, cosU1*cosSigma-sinU1*sinSigma*cosA1)
|
RadiansToDegrees(minLatL),
|
||||||
c := flattening / 16 * cosSqAlpha * (4 + flattening*(4-3*cosSqAlpha))
|
nil
|
||||||
lam := lamda - (1-c)*flattening*sinAlpha*
|
|
||||||
(sigma+c*sinSigma*(cos25SigmaM+c*cosSigma*(-1+2*cos25SigmaM*cos25SigmaM)))
|
|
||||||
|
|
||||||
rvlon := lon + RadiansToDegrees(lam)
|
|
||||||
rvlat := RadiansToDegrees(lat2)
|
|
||||||
|
|
||||||
return rvlon, rvlat
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// normalizeLon normalizes a longitude value within the -180 to 180 range
|
func checkLatitude(latitude float64) error {
|
||||||
func normalizeLon(lonDeg float64) float64 {
|
if math.IsNaN(latitude) || latitude < minLat || latitude > maxLat {
|
||||||
if lonDeg >= -180 && lonDeg <= 180 {
|
return fmt.Errorf("invalid latitude %f; must be between %f and %f", latitude, minLat, maxLat)
|
||||||
return lonDeg
|
|
||||||
}
|
}
|
||||||
|
return nil
|
||||||
off := math.Mod(lonDeg+180, 360)
|
|
||||||
if off < 0 {
|
|
||||||
return 180 + off
|
|
||||||
} else if off == 0 && lonDeg > 0 {
|
|
||||||
return 180
|
|
||||||
}
|
|
||||||
return -180 + off
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// normalizeLat normalizes a latitude value within the -90 to 90 range
|
func checkLongitude(longitude float64) error {
|
||||||
func normalizeLat(latDeg float64) float64 {
|
if math.IsNaN(longitude) || longitude < minLon || longitude > maxLon {
|
||||||
if latDeg >= -90 && latDeg <= 90 {
|
return fmt.Errorf("invalid longitude %f; must be between %f and %f", longitude, minLon, maxLon)
|
||||||
return latDeg
|
|
||||||
}
|
}
|
||||||
off := math.Abs(math.Mod(latDeg+90, 360))
|
return nil
|
||||||
if off <= 180 {
|
|
||||||
return off - 90
|
|
||||||
}
|
|
||||||
return (360 - off) - 90
|
|
||||||
}
|
}
|
||||||
|
@ -83,13 +83,16 @@ func TestScaleLatUnscaleLat(t *testing.T) {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
func TestComputeBoundingBoxCheckLatitudeAtEquator(t *testing.T) {
|
func TestRectFromPointDistance(t *testing.T) {
|
||||||
// at the equator 1 degree of latitude is about 110567 meters
|
// at the equator 1 degree of latitude is about 110567 meters
|
||||||
_, upperLeftLat, _, lowerRightLat := ComputeBoundingBox(0, 0, 110567)
|
_, upperLeftLat, _, lowerRightLat, err := RectFromPointDistance(0, 0, 110567)
|
||||||
if math.Abs(upperLeftLat-1) > 1E-4 {
|
if err != nil {
|
||||||
|
t.Fatal(err)
|
||||||
|
}
|
||||||
|
if math.Abs(upperLeftLat-1) > 1E-2 {
|
||||||
t.Errorf("expected bounding box upper left lat to be almost 1, got %f", upperLeftLat)
|
t.Errorf("expected bounding box upper left lat to be almost 1, got %f", upperLeftLat)
|
||||||
}
|
}
|
||||||
if math.Abs(lowerRightLat+1) > 1E-4 {
|
if math.Abs(lowerRightLat+1) > 1E-2 {
|
||||||
t.Errorf("expected bounding box lower right lat to be almost -1, got %f", lowerRightLat)
|
t.Errorf("expected bounding box lower right lat to be almost -1, got %f", lowerRightLat)
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@ -178,45 +181,3 @@ func TestBoundingBoxContains(t *testing.T) {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
func TestNormalizeLon(t *testing.T) {
|
|
||||||
tests := []struct {
|
|
||||||
lon float64
|
|
||||||
want float64
|
|
||||||
}{
|
|
||||||
{-180, -180},
|
|
||||||
{0, 0},
|
|
||||||
{180, 180},
|
|
||||||
{181, -179},
|
|
||||||
{-181, 179},
|
|
||||||
{540, 180},
|
|
||||||
}
|
|
||||||
|
|
||||||
for _, test := range tests {
|
|
||||||
got := normalizeLon(test.lon)
|
|
||||||
if test.want != got {
|
|
||||||
t.Errorf("expected normalizedLon %f, got %f for %f", test.want, got, test.lon)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
func TestNormalizeLat(t *testing.T) {
|
|
||||||
tests := []struct {
|
|
||||||
lat float64
|
|
||||||
want float64
|
|
||||||
}{
|
|
||||||
{-90, -90},
|
|
||||||
{0, 0},
|
|
||||||
{90, 90},
|
|
||||||
// somewhat unexpected, but double-checked against lucene
|
|
||||||
{91, 89},
|
|
||||||
{-91, -89},
|
|
||||||
}
|
|
||||||
|
|
||||||
for _, test := range tests {
|
|
||||||
got := normalizeLat(test.lat)
|
|
||||||
if test.want != got {
|
|
||||||
t.Errorf("expected normalizedLat %f, got %f for %f", test.want, got, test.lat)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
@ -146,6 +146,12 @@ func earthDiameter(lat float64) float64 {
|
|||||||
return earthDiameterPerLatitude[int(index)]
|
return earthDiameterPerLatitude[int(index)]
|
||||||
}
|
}
|
||||||
|
|
||||||
|
var pio2 = math.Pi / 2
|
||||||
|
|
||||||
|
func sin(a float64) float64 {
|
||||||
|
return cos(a - pio2)
|
||||||
|
}
|
||||||
|
|
||||||
// cos is a sloppy math (faster) implementation of math.Cos
|
// cos is a sloppy math (faster) implementation of math.Cos
|
||||||
func cos(a float64) float64 {
|
func cos(a float64) float64 {
|
||||||
if a < 0.0 {
|
if a < 0.0 {
|
||||||
|
@ -24,10 +24,12 @@ import (
|
|||||||
func NewGeoPointDistanceSearcher(indexReader index.IndexReader, centerLon,
|
func NewGeoPointDistanceSearcher(indexReader index.IndexReader, centerLon,
|
||||||
centerLat, dist float64, field string, boost float64,
|
centerLat, dist float64, field string, boost float64,
|
||||||
options search.SearcherOptions) (search.Searcher, error) {
|
options search.SearcherOptions) (search.Searcher, error) {
|
||||||
|
|
||||||
// compute bounding box containing the circle
|
// compute bounding box containing the circle
|
||||||
topLeftLon, topLeftLat, bottomRightLon, bottomRightLat :=
|
topLeftLon, topLeftLat, bottomRightLon, bottomRightLat, err :=
|
||||||
geo.ComputeBoundingBox(centerLon, centerLat, dist)
|
geo.RectFromPointDistance(centerLon, centerLat, dist)
|
||||||
|
if err != nil {
|
||||||
|
return nil, err
|
||||||
|
}
|
||||||
|
|
||||||
// build a searcher for the box
|
// build a searcher for the box
|
||||||
boxSearcher, err := boxSearcher(indexReader,
|
boxSearcher, err := boxSearcher(indexReader,
|
||||||
|
Loading…
Reference in New Issue
Block a user