I'm looking to port a portion of the dynamo-geo.js library to golang, in order to query the closest points (stored in DyanmoDB) to a given point.
Querying by radius is ideal, but if querying by rectangle is a more straightforward algorithm I'm cool with that too.
Here is the code I came up with to query by radius, but I can't seem to get a non-empty list of covering cells.
What is wrong with my algorithm?
// Query a circular area constructed by a center point and its radius.
// @see https://github.com/rh389/dynamodb-geo.js/blob/6c388b9070014a096885e00fff6c3fc933d9853f/src/GeoDataManager.ts#L229
func queryRadius(lat float64, lng float64, radiusMeters float64) (error) {
earthRadiusMeters := 6367000.0
// Step1: Get the bounding region (rectangle) from the center and the radius
// @see https://github.com/rh389/dynamodb-geo.js/blob/6c388b9070014a096885e00fff6c3fc933d9853f/src/s2/S2Util.ts#L23
centerLatLng := s2.LatLngFromDegrees(lat, lng)
latReferenceUnit := 1.0
if lat > 0.0 {
latReferenceUnit = -1.0
}
latReferenceLatLng := s2.LatLngFromDegrees(lat+latReferenceUnit, lng)
lngReferenceUnit := 1.0
if lng > 0.0 {
lngReferenceUnit = -1.0
}
lngReferenceLatLng := s2.LatLngFromDegrees(lat, lng+lngReferenceUnit)
latForRadius := radiusMeters / centerLatLng.Distance(latReferenceLatLng).Radians() * earthRadiusMeters
lngForRadius := radiusMeters / centerLatLng.Distance(lngReferenceLatLng).Radians() * earthRadiusMeters
minLatLng := s2.LatLngFromDegrees(lat-latForRadius, lng-lngForRadius)
maxLatLng := s2.LatLngFromDegrees(lat+latForRadius, lng+lngForRadius)
boundingRect := s2.RectFromLatLng(minLatLng)
boundingRect = boundingRect.AddPoint(maxLatLng)
// Step2: Compute the CellIDs for the region we want to cover.
// defaults per https://github.com/vekexasia/nodes2-ts/blob/1952d8c1f6cb4a862731ace2d5f74d472ec22e55/src/S2RegionCoverer.ts#L101
rc := &s2.RegionCoverer{MaxLevel: 30, MaxCells: 8, LevelMod: 1}
r := s2.Region(boundingRect.CapBound())
coveringCells := rc.Covering(r)
for _, c := range coveringCells {
log.WithFields(log.Fields{
"Covering Cell": c,
}).Info("=>")
}
return nil
}
Note: This is an answer to my specific original question, however I've been unable to find a complete solution to the problem I'm trying to solve (query the closest points, stored in DyanmoDB, to a given point using S2). I'll update this answer if/when I get a complete solution. I'm currently blocked on this issue. Any help appreciated.
Here is a complete go program that computes the covering cells from a point (in degrees) and a radius (in meters).
FWIW the algorithm to determine the bounding square, from a point and a radius, is not very accurate. SO to Martin F for providing the bounding box algorithm.
package main
import (
"fmt"
"math"
"strconv"
"github.com/golang/geo/s2"
)
const earthRadiusM = 6371000 // per https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
const hashLength = 8 // < 1km per https://github.com/rh389/dynamodb-geo.js/blob/master/test/integration/hashKeyLength.ts
func main() {
lowPrefix := uint64(0)
highPrefix := uint64(0)
ctrLat := 52.225730 // Cambridge UK
ctrLng := 0.149593
boundingSq := squareFromCenterAndRadius(ctrLat, ctrLng, 500)
fmt.Printf("
Bounding sq %+v
", boundingSq)
coveringCells := getCoveringCells(boundingSq)
fmt.Printf("Covering Cells (%d):
", len(coveringCells))
for idx, cell := range coveringCells {
// cell is the UUID of the center of this cell
fullHash, hashPrefix := genCellIntPrefix(cell)
if 0 == idx {
lowPrefix = hashPrefix
highPrefix = hashPrefix
} else if hashPrefix < lowPrefix {
lowPrefix = hashPrefix
} else if hashPrefix > highPrefix {
highPrefix = hashPrefix
}
fmt.Printf("\tID:%19v uint64: %-19d prefix: %-10d Range: %-19d - %-19d
", cell, fullHash, hashPrefix, uint64(cell.RangeMin()), uint64(cell.RangeMax()))
}
fmt.Printf("\tPrefix Range from loop: %-10d - %-10d
", lowPrefix, highPrefix)
// TODO: Assuming covering cells are sorted. Correct assumption?
_, lowPrefix = genCellIntPrefix(coveringCells[0].RangeMin())
_, highPrefix = genCellIntPrefix(coveringCells[len(coveringCells)-1].RangeMax())
fmt.Printf("\tPrefix Range direct: %-10d - %-10d
", lowPrefix, highPrefix)
}
// Get bounding box square from center point and radius
// Boundnig box is not extremely accurate to the radiusMeters passed in
// @see https://gis.stackexchange.com/questions/80809/calculating-bounding-box-coordinates-based-on-center-and-radius
func squareFromCenterAndRadius(centerLatDegrees float64, centerLngDegrees float64, radiusMeters float32) s2.Rect {
latLng := s2.LatLngFromDegrees(centerLatDegrees, centerLngDegrees)
deltaLng := float64(360 * radiusMeters / earthRadiusM) //Search Radius, difference in lat
deltaLat := deltaLng * math.Cos(latLng.Lng.Radians()) //Search Radius, difference in lng
lowerLeftLatDeg := centerLatDegrees - deltaLat
lowerLeftLngDeg := centerLngDegrees - deltaLng
lowerLeft := s2.LatLngFromDegrees(lowerLeftLatDeg, lowerLeftLngDeg) // AKA s2.Rect.Lo
upperRightLatDeg := centerLatDegrees + deltaLat
upperRightLngDeg := centerLngDegrees + deltaLng
upperRight := s2.LatLngFromDegrees(upperRightLatDeg, upperRightLngDeg) // AKA s2.Rect.Hi
boundingSquare := s2.RectFromLatLng(lowerLeft).AddPoint(upperRight)
return boundingSquare
}
func getCoveringCells(boundingRect s2.Rect) s2.CellUnion {
// defaults per https://github.com/vekexasia/nodes2-ts/blob/1952d8c1f6cb4a862731ace2d5f74d472ec22e55/src/S2RegionCoverer.ts#L101
rc := &s2.RegionCoverer{
MinLevel: 12, // 3km^2 per http://s2geometry.io/resources/s2cell_statistics
MaxLevel: 20, // 46m^2 per http://s2geometry.io/resources/s2cell_statistics
MaxCells: 8,
LevelMod: 1,
}
return rc.Covering(boundingRect)
}
func genCellIntPrefix(cell s2.CellID) (hash uint64, prefix uint64) {
hash = uint64(cell)
geohashString := strconv.FormatUint(hash, 10)
denominator := math.Pow10(len(geohashString) - hashLength)
prefix = hash / uint64(denominator)
return
}
I am by no means a GIS expert, so any suggestions/comments for improvement are welcome.