diff --git a/internal/geo/geojson.go b/internal/geo/geojson.go new file mode 100644 index 0000000..3baea7c --- /dev/null +++ b/internal/geo/geojson.go @@ -0,0 +1,107 @@ +package geo + +import ( + "encoding/json" + "fmt" +) + +type geometry struct { + Type string `json:"type"` + Coordinates json.RawMessage `json:"coordinates"` +} + +// ContainsPoint reports whether a GeoJSON Polygon or MultiPolygon contains p. +// GeoJSON coordinate order is [longitude, latitude]. +func ContainsPoint(raw []byte, p Point) (bool, error) { + if len(raw) == 0 { + return false, fmt.Errorf("geojson geometry is empty") + } + + var geom geometry + if err := json.Unmarshal(raw, &geom); err != nil { + return false, fmt.Errorf("decode geojson geometry: %w", err) + } + + switch geom.Type { + case "Polygon": + polygon, err := decodePolygon(geom.Coordinates) + if err != nil { + return false, fmt.Errorf("decode polygon: %w", err) + } + return polygonContainsPoint(polygon, p), nil + case "MultiPolygon": + multiPolygon, err := decodeMultiPolygon(geom.Coordinates) + if err != nil { + return false, fmt.Errorf("decode multipolygon: %w", err) + } + for _, polygon := range multiPolygon { + if polygonContainsPoint(polygon, p) { + return true, nil + } + } + return false, nil + case "": + return false, fmt.Errorf("geojson geometry type is required") + default: + return false, fmt.Errorf("unsupported geojson geometry type %q", geom.Type) + } +} + +func decodePolygon(raw json.RawMessage) (Polygon, error) { + var coords [][][]float64 + if err := json.Unmarshal(raw, &coords); err != nil { + return nil, err + } + return polygonFromCoordinates(coords) +} + +func decodeMultiPolygon(raw json.RawMessage) ([]Polygon, error) { + var coords [][][][]float64 + if err := json.Unmarshal(raw, &coords); err != nil { + return nil, err + } + if len(coords) == 0 { + return nil, fmt.Errorf("multipolygon has no polygons") + } + + out := make([]Polygon, 0, len(coords)) + for i, polygonCoords := range coords { + polygon, err := polygonFromCoordinates(polygonCoords) + if err != nil { + return nil, fmt.Errorf("polygons[%d]: %w", i, err) + } + out = append(out, polygon) + } + return out, nil +} + +func polygonFromCoordinates(coords [][][]float64) (Polygon, error) { + if len(coords) == 0 { + return nil, fmt.Errorf("polygon has no rings") + } + + polygon := make(Polygon, 0, len(coords)) + for i, ringCoords := range coords { + ring, err := ringFromCoordinates(ringCoords) + if err != nil { + return nil, fmt.Errorf("rings[%d]: %w", i, err) + } + polygon = append(polygon, ring) + } + return polygon, nil +} + +func ringFromCoordinates(coords [][]float64) (Ring, error) { + if len(coords) == 0 { + return nil, fmt.Errorf("ring has no points") + } + + ring := make(Ring, 0, len(coords)) + for i, pair := range coords { + if len(pair) < 2 { + return nil, fmt.Errorf("points[%d] has %d values, need longitude and latitude", i, len(pair)) + } + ring = append(ring, Point{Longitude: pair[0], Latitude: pair[1]}) + } + return ring, nil +} diff --git a/internal/geo/point.go b/internal/geo/point.go new file mode 100644 index 0000000..6a740c4 --- /dev/null +++ b/internal/geo/point.go @@ -0,0 +1,105 @@ +package geo + +import "math" + +const epsilon = 1e-9 + +// Point is a geographic coordinate in decimal degrees. +type Point struct { + Longitude float64 + Latitude float64 +} + +// Ring is one GeoJSON linear ring. +type Ring []Point + +// Polygon is a GeoJSON polygon. The first ring is the exterior ring; subsequent +// rings are holes. +type Polygon []Ring + +func polygonContainsPoint(polygon Polygon, p Point) bool { + if len(polygon) == 0 { + return false + } + if pointOnRing(polygon[0], p) { + return true + } + if !ringContainsPoint(polygon[0], p) { + return false + } + for _, hole := range polygon[1:] { + if pointOnRing(hole, p) { + return true + } + if ringContainsPoint(hole, p) { + return false + } + } + return true +} + +func ringContainsPoint(ring Ring, p Point) bool { + inside := false + n := len(ring) + if n == 0 { + return false + } + + for i, j := 0, n-1; i < n; j, i = i, i+1 { + a := ring[j] + b := ring[i] + if pointOnSegment(p, a, b) { + return true + } + + intersects := (a.Latitude > p.Latitude) != (b.Latitude > p.Latitude) + if intersects { + x := (b.Longitude-a.Longitude)*(p.Latitude-a.Latitude)/(b.Latitude-a.Latitude) + a.Longitude + if almostEqual(x, p.Longitude) { + return true + } + if x > p.Longitude { + inside = !inside + } + } + } + return inside +} + +func pointOnRing(ring Ring, p Point) bool { + n := len(ring) + if n == 0 { + return false + } + for i, j := 0, n-1; i < n; j, i = i, i+1 { + if pointOnSegment(p, ring[j], ring[i]) { + return true + } + } + return false +} + +func pointOnSegment(p, a, b Point) bool { + cross := (p.Latitude-a.Latitude)*(b.Longitude-a.Longitude) - (p.Longitude-a.Longitude)*(b.Latitude-a.Latitude) + if math.Abs(cross) > epsilon { + return false + } + + minLon, maxLon := minMax(a.Longitude, b.Longitude) + minLat, maxLat := minMax(a.Latitude, b.Latitude) + return p.Longitude >= minLon-epsilon && + p.Longitude <= maxLon+epsilon && + p.Latitude >= minLat-epsilon && + p.Latitude <= maxLat+epsilon +} + +func minMax(a, b float64) (float64, float64) { + if a < b { + return a, b + } + return b, a +} + +func almostEqual(a, b float64) bool { + return math.Abs(a-b) <= epsilon +} diff --git a/internal/geo/point_test.go b/internal/geo/point_test.go new file mode 100644 index 0000000..62d8f32 --- /dev/null +++ b/internal/geo/point_test.go @@ -0,0 +1,181 @@ +package geo + +import ( + "encoding/json" + "strings" + "testing" +) + +const squarePolygon = `{ + "type": "Polygon", + "coordinates": [[ + [-91.0, 38.0], + [-90.0, 38.0], + [-90.0, 39.0], + [-91.0, 39.0], + [-91.0, 38.0] + ]] +}` + +func TestContainsPointInsideSimplePolygon(t *testing.T) { + got, err := ContainsPoint([]byte(squarePolygon), Point{Longitude: -90.5, Latitude: 38.5}) + if err != nil { + t.Fatalf("ContainsPoint() error = %v", err) + } + if !got { + t.Fatalf("ContainsPoint() = false, want true") + } +} + +func TestContainsPointAcceptsRawMessage(t *testing.T) { + got, err := ContainsPoint(json.RawMessage(squarePolygon), Point{Longitude: -90.5, Latitude: 38.5}) + if err != nil { + t.Fatalf("ContainsPoint() error = %v", err) + } + if !got { + t.Fatalf("ContainsPoint() = false, want true") + } +} + +func TestContainsPointOutsideSimplePolygon(t *testing.T) { + got, err := ContainsPoint([]byte(squarePolygon), Point{Longitude: -89.5, Latitude: 38.5}) + if err != nil { + t.Fatalf("ContainsPoint() error = %v", err) + } + if got { + t.Fatalf("ContainsPoint() = true, want false") + } +} + +func TestContainsPointOnBoundary(t *testing.T) { + got, err := ContainsPoint([]byte(squarePolygon), Point{Longitude: -91.0, Latitude: 38.5}) + if err != nil { + t.Fatalf("ContainsPoint() error = %v", err) + } + if !got { + t.Fatalf("ContainsPoint() = false, want true") + } +} + +func TestContainsPointInHoleReturnsFalse(t *testing.T) { + const polygonWithHole = `{ + "type": "Polygon", + "coordinates": [ + [[0,0],[10,0],[10,10],[0,10],[0,0]], + [[4,4],[6,4],[6,6],[4,6],[4,4]] + ] +}` + + got, err := ContainsPoint([]byte(polygonWithHole), Point{Longitude: 5, Latitude: 5}) + if err != nil { + t.Fatalf("ContainsPoint() error = %v", err) + } + if got { + t.Fatalf("ContainsPoint() = true, want false") + } +} + +func TestContainsPointOnHoleBoundaryReturnsTrue(t *testing.T) { + const polygonWithHole = `{ + "type": "Polygon", + "coordinates": [ + [[0,0],[10,0],[10,10],[0,10],[0,0]], + [[4,4],[6,4],[6,6],[4,6],[4,4]] + ] +}` + + got, err := ContainsPoint([]byte(polygonWithHole), Point{Longitude: 4, Latitude: 5}) + if err != nil { + t.Fatalf("ContainsPoint() error = %v", err) + } + if !got { + t.Fatalf("ContainsPoint() = false, want true") + } +} + +func TestContainsPointInsideOneMultiPolygonMember(t *testing.T) { + const multiPolygon = `{ + "type": "MultiPolygon", + "coordinates": [ + [[[0,0],[1,0],[1,1],[0,1],[0,0]]], + [[[10,10],[12,10],[12,12],[10,12],[10,10]]] + ] +}` + + got, err := ContainsPoint([]byte(multiPolygon), Point{Longitude: 11, Latitude: 11}) + if err != nil { + t.Fatalf("ContainsPoint() error = %v", err) + } + if !got { + t.Fatalf("ContainsPoint() = false, want true") + } +} + +func TestContainsPointUsesLongitudeLatitudeOrder(t *testing.T) { + const narrowPolygon = `{ + "type": "Polygon", + "coordinates": [[ + [-91.0, 38.0], + [-90.0, 38.0], + [-90.0, 39.0], + [-91.0, 39.0], + [-91.0, 38.0] + ]] +}` + + got, err := ContainsPoint([]byte(narrowPolygon), Point{Longitude: -90.5, Latitude: 38.5}) + if err != nil { + t.Fatalf("ContainsPoint() error = %v", err) + } + if !got { + t.Fatalf("ContainsPoint() = false, want true") + } + + got, err = ContainsPoint([]byte(narrowPolygon), Point{Longitude: 38.5, Latitude: -90.5}) + if err != nil { + t.Fatalf("ContainsPoint() reversed error = %v", err) + } + if got { + t.Fatalf("ContainsPoint() with reversed coordinate values = true, want false") + } +} + +func TestContainsPointUnsupportedGeometryError(t *testing.T) { + _, err := ContainsPoint([]byte(`{"type":"Point","coordinates":[-90,38]}`), Point{Longitude: -90, Latitude: 38}) + if err == nil { + t.Fatalf("ContainsPoint() error = nil, want error") + } + if !strings.Contains(err.Error(), `unsupported geojson geometry type "Point"`) { + t.Fatalf("ContainsPoint() error = %q", err) + } +} + +func TestContainsPointInvalidJSONError(t *testing.T) { + _, err := ContainsPoint([]byte(`{"type":"Polygon"`), Point{Longitude: -90, Latitude: 38}) + if err == nil { + t.Fatalf("ContainsPoint() error = nil, want error") + } + if !strings.Contains(err.Error(), "decode geojson geometry") { + t.Fatalf("ContainsPoint() error = %q", err) + } +} + +func TestContainsPointMalformedCoordinatesError(t *testing.T) { + _, err := ContainsPoint([]byte(`{"type":"Polygon","coordinates":[[[1]]]}`), Point{Longitude: 1, Latitude: 1}) + if err == nil { + t.Fatalf("ContainsPoint() error = nil, want error") + } + if !strings.Contains(err.Error(), "need longitude and latitude") { + t.Fatalf("ContainsPoint() error = %q", err) + } +} + +func TestContainsPointEmptyRingError(t *testing.T) { + _, err := ContainsPoint([]byte(`{"type":"Polygon","coordinates":[[]]}`), Point{Longitude: 1, Latitude: 1}) + if err == nil { + t.Fatalf("ContainsPoint() error = nil, want error") + } + if !strings.Contains(err.Error(), "ring has no points") { + t.Fatalf("ContainsPoint() error = %q", err) + } +}