将Golang S2 Geometry库与dynamodb一起使用

Question

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
}

Answer 1

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.