Last active
July 24, 2019 18:38
-
-
Save rynop/1f695feb343d7113ed16c82893fac1d9 to your computer and use it in GitHub Desktop.
S2: get bounding box of center point (in degrees) and radius
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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("\nBounding sq %+v\n", boundingSq) | |
coveringCells := getCoveringCells(boundingSq) | |
fmt.Printf("Covering Cells (%d):\n", 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\n", cell, fullHash, hashPrefix, uint64(cell.RangeMin()), uint64(cell.RangeMax())) | |
} | |
fmt.Printf("\tPrefix Range from loop: %-10d - %-10d\n", 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\n", 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 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment