Updated geojson packages

This commit is contained in:
tidwall 2018-10-27 09:23:29 -07:00
parent b17bbbd829
commit 745579b56b
8 changed files with 372 additions and 55 deletions

4
Gopkg.lock generated
View File

@ -227,7 +227,7 @@
[[projects]]
branch = "master"
digest = "1:4be7626fb8f801eb85aa7494ce3a504c9b4121a07f4ec19d7d204185bd6397d5"
digest = "1:c5ac96e72d3ff6694602f3273dd71ef04a67c9591465aac92dc1aa8c821b8f91"
name = "github.com/tidwall/geojson"
packages = [
".",
@ -235,7 +235,7 @@
"geometry",
]
pruneopts = ""
revision = "928ede3da18d831dea0af0bb26adeb025145c23b"
revision = "32782c39ca84f98113436a297f14601e4fee527d"
[[projects]]
digest = "1:3ddca2bd5496c6922a2a9e636530e178a43c2a534ea6634211acdc7d10222794"

1
vendor/github.com/tidwall/geojson/.travis.yml generated vendored Normal file
View File

@ -0,0 +1 @@
language: go

View File

@ -1,6 +1,6 @@
# `GeoJSON`
[![GoDoc](https://godoc.org/github.com/tidwall/geojson?status.svg)](https://godoc.org/github.com/tidwall/geojson)
[![Build Status](https://travis-ci.org/tidwall/geojson.svg?branch=master)](https://travis-ci.org/tidwall/geojson) [![GoDoc](https://godoc.org/github.com/tidwall/geojson?status.svg)](https://godoc.org/github.com/tidwall/geojson)
This package provides GeoJSON utilties for Go. It's designed for [Tile38](https://github.com/tidwall/tile38).

View File

@ -10,11 +10,12 @@ import (
// Circle ...
type Circle struct {
Object
center geometry.Point
meters float64
steps int
km bool
extra *extra
center geometry.Point
meters float64
haversine float64
steps int
km bool
extra *extra
}
// NewCircle returns an circle object
@ -29,6 +30,7 @@ func NewCircle(center geometry.Point, meters float64, steps int) *Circle {
if meters <= 0 {
g.Object = NewPoint(center)
} else {
meters = geo.NormalizeDistance(meters)
var points []geometry.Point
step := 360.0 / float64(steps)
i := 0
@ -44,6 +46,7 @@ func NewCircle(center geometry.Point, meters float64, steps int) *Circle {
g.Object = NewPolygon(
geometry.NewPoly(points, nil, geometry.DefaultIndexOptions),
)
g.haversine = geo.DistanceToHaversine(meters)
}
return g
}
@ -71,10 +74,110 @@ func (g *Circle) String() string {
return string(g.AppendJSON(nil))
}
// Meters returns the circle's radius
func (g *Circle) Meters() float64 {
return g.meters
}
// Center returns the circle's center point
func (g *Circle) Center() geometry.Point {
return g.center
}
// Within returns true if circle is contained inside object
func (g *Circle) Within(obj Object) bool {
return obj.Contains(g)
}
func (g *Circle) contains(p geometry.Point, allowOnEdge bool) bool {
h := geo.Haversine(p.Y, p.X, g.center.Y, g.center.X)
if allowOnEdge {
return h <= g.haversine
}
return h < g.haversine
}
// Contains returns true if the circle contains other object
func (g *Circle) Contains(obj Object) bool {
switch other := obj.(type) {
case *Point:
return g.contains(other.Center(), false)
case *Circle:
return other.Distance(g) < (other.meters + g.meters)
case *LineString:
for i := 0; i < other.base.NumPoints(); i++ {
if geoDistancePoints(other.base.PointAt(i), g.center) > g.meters {
return false
}
}
return true
case Collection:
for _, p := range other.Children() {
if !g.Contains(p) {
return false
}
}
return true
default:
// No simple cases, so using polygon approximation.
return g.Object.Contains(other)
}
}
func (g *Circle) intersectsSegment(seg geometry.Segment) bool {
start, end := seg.A, seg.B
// These are faster checks. If they succeed there's no need do complicate things.
if g.contains(start, true) {
return true
}
if g.contains(end, true) {
return true
}
// Distance between start and end
l := geo.DistanceTo(start.Y, start.X, end.Y, end.X)
// Unit direction vector
dx := (end.X - start.X) / l
dy := (end.Y - start.Y) / l
// Point of the line closest to the center
t := dx*(g.center.X-start.X) + dy*(g.center.Y-start.Y)
px := t*dx + start.X
py := t*dy + start.Y
if px < start.X || px > end.X || py < start.Y || py > end.Y {
// closest point is outside the segment
return false
}
// Distance from the closest point to the center
return g.contains(geometry.Point{X: px, Y: py}, true)
}
// Intersects returns true the circle intersects other object
func (g *Circle) Intersects(obj Object) bool {
switch other := obj.(type) {
case *Point:
return g.contains(other.Center(), true)
case *Circle:
return other.Distance(g) <= (other.meters + g.meters)
case *LineString:
for i := 0; i < other.base.NumSegments(); i++ {
if g.intersectsSegment(other.base.SegmentAt(i)) {
return true
}
}
return false
case Collection:
for _, p := range other.Children() {
if g.Intersects(p) {
return true
}
}
return false
default:
// No simple cases, so using polygon approximation.
return g.Object.Intersects(obj)
}
}

View File

@ -1,22 +1,173 @@
package geojson
import "testing"
import (
"testing"
func TestCircle(t *testing.T) {
"github.com/tidwall/geojson/geometry"
)
func TestCircleNew(t *testing.T) {
expectJSON(t,
`{"type":"Feature","geometry":{"type":"Point","coordinates":[-112,33]},"properties":{"type":"Circle","radius":"5000"}}`,
`{"type":"Feature","geometry":{"type":"Point","coordinates":[-112,33]},"properties":{"type":"Circle","radius":5000,"radius_units":"m"}}`,
)
g, err := Parse(`{
"type":"Feature",
"geometry":{"type":"Point","coordinates":[-112.2693,33.5123]},
"properties": {
"type": "Circle",
"radius": 1000
}
}`, nil)
"type":"Feature",
"geometry":{"type":"Point","coordinates":[-112.2693,33.5123]},
"properties": {
"type": "Circle",
"radius": 1000
}
}`, nil)
if err != nil {
t.Fatal(err)
}
expect(t, g.Contains(PO(-112.26, 33.51)))
circle := NewCircle(P(-112, 33), 123456.654321, 64)
expectJSON(t, circle.JSON(), `{"type":"Feature","geometry":{"type":"Point","coordinates":[-112,33]},"properties":{"type":"Circle","radius":123456.654321,"radius_units":"m"}}`)
}
func TestCircleContains(t *testing.T) {
g := NewCircle(P(-122.4412, 37.7335), 1000, 64)
expect(t, g.Contains(PO(-122.4412, 37.7335)))
expect(t, g.Contains(PO(-122.44121, 37.7335)))
expect(t, g.Contains(
MPO([]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.4408378, 37.733)})))
expect(t, g.Contains(
NewCircle(P(-122.44121, 37.7335), 500, 64)))
expect(t, g.Contains(
LO([]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.4408378, 37.733)})))
expect(t, g.Contains(
MLO([]*geometry.Line{
L([]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.4408378, 37.733),
}),
L([]geometry.Point{
P(-122.44, 37.733),
P(-122.44, 37.7341129),
})})))
expect(t, g.Contains(
PPO(
[]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.4408378, 37.733),
P(-122.44, 37.733),
P(-122.44, 37.7341129),
P(-122.4408378, 37.7341129),
},
[][]geometry.Point{})))
// Does-not-contain
expect(t, !g.Contains(PO(-122.265, 37.826)))
expect(t, !g.Contains(
NewCircle(P(-122.265, 37.826), 100, 64)))
expect(t, !g.Contains(
LO([]geometry.Point{
P(-122.265, 37.826),
P(-122.210, 37.860)})))
expect(t, !g.Contains(
MPO([]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.198181, 37.7490)})))
expect(t, !g.Contains(
MLO([]*geometry.Line{
L([]geometry.Point{
P(-122.265, 37.826),
P(-122.265, 37.860),
}),
L([]geometry.Point{
P(-122.44, 37.733),
P(-122.44, 37.7341129),
})})))
expect(t, !g.Contains(PPO(
[]geometry.Point{
P(-122.265, 37.826),
P(-122.265, 37.860),
P(-122.210, 37.860),
P(-122.210, 37.826),
P(-122.265, 37.826),
},
[][]geometry.Point{})))
}
func TestCircleIntersects(t *testing.T) {
g := NewCircle(P(-122.4412, 37.7335), 1000, 64)
expect(t, g.Intersects(PO(-122.4412, 37.7335)))
expect(t, g.Intersects(PO(-122.44121, 37.7335)))
expect(t, g.Intersects(
MPO([]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.4408378, 37.733)})))
expect(t, g.Intersects(
NewCircle(P(-122.44121, 37.7335), 500, 64)))
expect(t, g.Intersects(
LO([]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.4408378, 37.733)})))
expect(t, g.Intersects(
LO([]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.265, 37.826)})))
expect(t, g.Intersects(
MLO([]*geometry.Line{
L([]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.265, 37.826),
}),
L([]geometry.Point{
P(-122.44, 37.733),
P(-122.44, 37.7341129),
})})))
expect(t, g.Intersects(
PPO(
[]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.265, 37.860),
P(-122.210, 37.826),
P(-122.44, 37.7341129),
P(-122.4408378, 37.7341129),
},
[][]geometry.Point{})))
expect(t, g.Intersects(
MPO([]geometry.Point{
P(-122.4408378, 37.7341129),
P(-122.198181, 37.7490)})))
expect(t, g.Intersects(
MLO([]*geometry.Line{
L([]geometry.Point{
P(-122.265, 37.826),
P(-122.265, 37.860),
}),
L([]geometry.Point{
P(-122.44, 37.733),
P(-122.44, 37.7341129),
})})))
// Does-not-intersect
expect(t, !g.Intersects(PO(-122.265, 37.826)))
expect(t, !g.Intersects(
NewCircle(P(-122.265, 37.826), 100, 64)))
expect(t, !g.Intersects(
LO([]geometry.Point{
P(-122.265, 37.826),
P(-122.210, 37.860)})))
}
// This snippet tests 100M comparisons.
// On my box this takes 24.5s without haversine trick, and 13.7s with the trick.
//
//func TestCircle_Performance(t *testing.T) {
// g := NewCircle(P(-122.4412, 37.7335), 1000, 64)
// r := rand.New(rand.NewSource(42))
// for i:= 0; i < 100000000; i++ {
// g.Contains(PO((r.Float64() - 0.5) * 180, r.Float64() * 90))
// }
// expect(t, true)
//}

View File

@ -12,18 +12,42 @@ const (
earthRadius = 6371e3
radians = math.Pi / 180
degrees = 180 / math.Pi
piR = math.Pi * earthRadius
twoPiR = 2 * piR
)
// DistanceTo return the distance in meteres between two point.
func DistanceTo(latA, lonA, latB, lonB float64) (meters float64) {
// Haversine ...
func Haversine(latA, lonA, latB, lonB float64) float64 {
φ1 := latA * radians
λ1 := lonA * radians
φ2 := latB * radians
λ2 := lonB * radians
Δφ := φ2 - φ1
Δλ := λ2 - λ1
a := math.Sin(Δφ/2)*math.Sin(Δφ/2) +
math.Cos(φ1)*math.Cos(φ2)*math.Sin(Δλ/2)*math.Sin(Δλ/2)
sΔφ2 := math.Sin(Δφ / 2)
sΔλ2 := math.Sin(Δλ / 2)
return sΔφ2*sΔφ2 + math.Cos(φ1)*math.Cos(φ2)*sΔλ2*sΔλ2
}
// NormalizeDistance ...
func NormalizeDistance(meters float64) float64 {
m1 := math.Mod(meters, twoPiR)
if m1 <= piR {
return m1
}
return twoPiR - m1
}
// DistanceToHaversine ...
func DistanceToHaversine(meters float64) float64 {
// convert the given distance to its haversine
sin := math.Sin(0.5 * meters / earthRadius)
return sin * sin
}
// DistanceTo return the distance in meteres between two point.
func DistanceTo(latA, lonA, latB, lonB float64) (meters float64) {
a := Haversine(latA, lonA, latB, lonB)
c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
return earthRadius * c
}
@ -60,36 +84,3 @@ func BearingTo(latA, lonA, latB, lonB float64) float64 {
return math.Mod(θ*degrees+360, 360)
}
// // SegmentIntersectsCircle ...
// func SegmentIntersectsCircle(
// startLat, startLon, endLat, endLon, centerLat, centerLon, meters float64,
// ) bool {
// // These are faster checks.
// // If they succeed there's no need do complicate things.
// if DistanceTo(startLat, startLon, centerLat, centerLon) <= meters {
// return true
// }
// if DistanceTo(endLat, endLon, centerLat, centerLon) <= meters {
// return true
// }
// // Distance between start and end
// l := DistanceTo(startLat, startLon, endLat, endLon)
// // Unit direction vector
// dLat := (endLat - startLat) / l
// dLon := (endLon - startLon) / l
// // Point of the line closest to the center
// t := dLon*(centerLon-startLon) + dLat*(centerLat-startLat)
// pLat := t*dLat + startLat
// pLon := t*dLon + startLon
// if pLon < startLon || pLon > endLon || pLat < startLat || pLat > endLat {
// // closest point is outside the segment
// return false
// }
// // Distance from the closest point to the center
// return DistanceTo(centerLat, centerLon, pLat, pLon) <= meters
// }

View File

@ -4,7 +4,19 @@
package geo
import "testing"
import (
"math"
"math/rand"
"testing"
"time"
)
func init() {
seed := time.Now().UnixNano()
//seed = 1540656736244531000
println(seed)
rand.Seed(seed)
}
func TestGeoCalc(t *testing.T) {
dist := 172853.26908429610193707048892974853515625
@ -30,3 +42,50 @@ func TestGeoCalc(t *testing.T) {
t.Fatalf("expected '%v', got '%v'", lonB, value2)
}
}
func TestHaversine(t *testing.T) {
latA := rand.Float64()*180 - 90
lonA := rand.Float64()*360 - 180
start := time.Now()
for time.Since(start) < time.Second/4 {
for i := 0; i < 1000; i++ {
latB := rand.Float64()*180 - 90
lonB := rand.Float64()*360 - 180
latC := rand.Float64()*180 - 90
lonC := rand.Float64()*360 - 180
haver1 := Haversine(latA, lonA, latB, lonB)
haver2 := Haversine(latA, lonA, latC, lonC)
meters1 := DistanceTo(latA, lonA, latB, lonB)
meters2 := DistanceTo(latA, lonA, latC, lonC)
switch {
case haver1 < haver2:
if meters1 >= meters2 {
t.Fatalf("failed")
}
case haver1 == haver2:
if meters1 != meters2 {
t.Fatalf("failed")
}
case haver1 > haver2:
if meters1 <= meters2 {
t.Fatalf("failed")
}
}
}
}
}
func TestNormalizeDistance(t *testing.T) {
start := time.Now()
for time.Since(start) < time.Second/4 {
for i := 0; i < 1000; i++ {
meters1 := rand.Float64() * 100000000
meters2 := NormalizeDistance(meters1)
dist1 := math.Floor(DistanceToHaversine(meters2) * 100000000.0)
dist2 := math.Floor(DistanceToHaversine(meters1) * 100000000.0)
if dist1 != dist2 {
t.Fatalf("expected %f, got %f", dist2, dist1)
}
}
}
}

View File

@ -31,6 +31,10 @@ func PO(x, y float64) *Point {
return NewPoint(P(x, y))
}
func MPO(points []geometry.Point) *MultiPoint {
return NewMultiPoint(points)
}
func RO(minX, minY, maxX, maxY float64) *Rect {
return NewRect(R(minX, minY, maxX, maxY))
}
@ -39,6 +43,14 @@ func LO(points []geometry.Point) *LineString {
return NewLineString(geometry.NewLine(points, nil))
}
func L(points []geometry.Point) *geometry.Line {
return geometry.NewLine(points, geometry.DefaultIndexOptions)
}
func MLO(lines []*geometry.Line) *MultiLineString {
return NewMultiLineString(lines)
}
func PPO(exterior []geometry.Point, holes [][]geometry.Point) *Polygon {
return NewPolygon(geometry.NewPoly(exterior, holes, nil))
}