diff --git a/measurement/measurement.go b/measurement/measurement.go index 16a20ee..bca7c33 100644 --- a/measurement/measurement.go +++ b/measurement/measurement.go @@ -11,7 +11,7 @@ import ( "github.com/tomchavakis/turf-go/geojson/geometry" "github.com/tomchavakis/turf-go/internal/common" "github.com/tomchavakis/turf-go/invariant" - meta "github.com/tomchavakis/turf-go/meta/coordAll" + "github.com/tomchavakis/turf-go/meta/coordAll" ) // Distance calculates the distance between two points in kilometers. This uses the Haversine formula @@ -825,20 +825,10 @@ func pointOnFeatureCollection(t interface{}, cent geometry.Point) bool { } } case geometry.Polygon: - for _, c := range gtp.Coordinates { - //todo: decide if cent is on surface - print(len(c.Coordinates)) - } - onSurface = true case geometry.MultiPolygon: - coords := gtp.Coordinates - for _, coord := range coords { - for _, pl := range coord.Coordinates { - //todo: decide if cent is on surface - print(len(pl.Coordinates)) - } + if pointOnPolygon(cent, gtp) { + onSurface = true } - onSurface = true } return onSurface || onBorderOrPoint @@ -850,3 +840,80 @@ func pointOnSegment(x float64, y float64, x1 float64, y1 float64, x2 float64, y2 var pb = math.Sqrt((x2 - x) * (x2 - x) + (y2 - y) * (y2 - y)) return ab == (ap + pb) } + +func pointOnPolygon(point geometry.Point, t interface{}) bool { + switch gpt := t.(type) { + case geometry.Polygon: + bbox, err := BBox(gpt) + if err != nil || len(bbox) == 0 || !inBBox(point, bbox) { + return false + } + return pointOnMultipolygon(point, []geometry.Polygon{gpt}) + case geometry.MultiPolygon: + bbox, err := BBox(gpt) + if err != nil || len(bbox) == 0 || !inBBox(point, bbox) { + return false + } + return pointOnMultipolygon(point, gpt.Coordinates) + } + + return false +} + +func pointOnMultipolygon(point geometry.Point, polys []geometry.Polygon) bool { + + insidePoly := false; + for i := 0; i < len(polys) && !insidePoly; i++ { + // check if it is in the outer ring first + if inRing(point, polys[i].Coordinates[0].Coordinates, false) { + inHole := false; + // check for the point in any of the holes + for j := 1; j < len(polys[i].Coordinates) && !inHole; j++ { + if inRing(point, polys[i].Coordinates[j].Coordinates, false) { + inHole = true; + } + } + if !inHole { + insidePoly = true; + } + } + } + + return insidePoly +} + +func inBBox(point geometry.Point, bbox []float64) bool { + return bbox[0] <= point.Lng && bbox[1] <= point.Lat && bbox[2] >= point.Lng && bbox[3] >= point.Lat +} + +func inRing(point geometry.Point, ring []geometry.Point, ignoreBoundary bool) bool { + isInside := false + + // remove last element if they are the same + if ring[0].Lat == ring[len(ring) - 1].Lat && ring[0].Lng == ring[len(ring) - 1].Lng { + ring = ring[:len(ring) - 1] + } + + for i, j := 0, len(ring) - 1 ; i < len(ring); i, j = i+1, i { + xi := ring[i].Lng + yi := ring[i].Lat + xj := ring[j].Lng + yj := ring[j].Lat + + onBoundary := point.Lat * (xi - xj) + yi * (xj - point.Lng) + yj * (point.Lng - xi) == 0 && + (xi - point.Lng) * (xj - point.Lng) <= 0 && + (yi - point.Lat) * (yj - point.Lat) <= 0 + + if onBoundary { + return !ignoreBoundary + } + + intersect := (yi > point.Lat) != (yj > point.Lat) && + (point.Lng < ((xj - xi) * (point.Lat - yi)) / (yj - yi) + xi) + + if intersect { + isInside = !isInside + } + } + return isInside +} diff --git a/measurement/measurement_test.go b/measurement/measurement_test.go index 6d79e15..9befd63 100644 --- a/measurement/measurement_test.go +++ b/measurement/measurement_test.go @@ -1560,3 +1560,171 @@ func TestRhumbDistance(t *testing.T) { }) } } + +func TestInBBox(t *testing.T) { + tests := map[string]struct { + point geometry.Point + bbox []float64 + want bool + }{ + "point in bbox": { + point: geometry.Point{ + Lng: 116.0, + Lat: -20.0, + }, + bbox: []float64{ + 113.0, -39.0, 154.0, -15.0, + }, + want: true, + }, + "point on bbox": { + point: geometry.Point{ + Lng: 116.0, + Lat: -20.0, + }, + bbox: []float64{ + 116.0, -20.0, 154.0, -15.0, + }, + want: true, + }, + "point not in bbox": { + point: geometry.Point{ + Lng: -20.0, + Lat: 116.0, + }, + bbox: []float64{ + 113.0, -39.0, 154.0, -15.0, + }, + want: false, + }, + } + for name, tt := range tests { + t.Run(name, func(t *testing.T) { + assert.Equal(t, tt.want, inBBox(tt.point, tt.bbox)) + }) + } +} + + +func TestInRing(t *testing.T) { + tests := map[string]struct { + point geometry.Point + ring []geometry.Point + want bool + } { + "point in ring": { + point: geometry.Point{ + Lng: 50.0, + Lat: 50.0, + }, + ring: []geometry.Point { + { + Lng: 0.0, + Lat: 0.0, + }, + { + Lng: 0.0, + Lat: 100.0, + }, + { + Lng: 100.0, + Lat: 100.0, + }, + { + Lng: 100.0, + Lat: 0.0, + }, + { + Lng: 0.0, + Lat: 0.0, + }, + }, + want: true, + }, + "point not in ring": { + point: geometry.Point{ + Lng: -50.0, + Lat: -50.0, + }, + ring: []geometry.Point { + { + Lng: 0.0, + Lat: 0.0, + }, + { + Lng: 0.0, + Lat: 100.0, + }, + { + Lng: 100.0, + Lat: 100.0, + }, + { + Lng: 100.0, + Lat: 0.0, + }, + { + Lng: 0.0, + Lat: 0.0, + }, + }, + want: false, + }, + } + + for name, tt := range tests { + t.Run(name, func(t *testing.T) { + assert.Equal(t, inRing(tt.point, tt.ring, false), tt.want) + }) + } +} + + +func TestInPolygon(t *testing.T) { + tests := map[string]struct { + point geometry.Point + polygon geometry.Polygon + want bool + } { + "point in polygon": { + point: geometry.Point{ + Lng: 2.0, + Lat: 102.0, + }, + polygon: geometry.Polygon { + Coordinates: []geometry.LineString{ + { + Coordinates: []geometry.Point{ + { + Lat: 2.0, + Lng: 102.0, + }, + { + Lat: 2.0, + Lng: 103.0, + }, + { + Lat: 3.0, + Lng: 103.0, + }, + { + Lat: 3.0, + Lng: 102.0, + }, + { + Lat: 2.0, + Lng: 102.0, + }, + }, + }, + }, + }, + want: true, + }, + } + for name, tt := range tests { + t.Run(name, func(t *testing.T) { + assert.Equal(t, pointOnPolygon(tt.point, tt.polygon), tt.want) + }) + } +}