diff --git a/types/geo/doc.go b/types/geo/doc.go new file mode 100644 index 000000000..749c63080 --- /dev/null +++ b/types/geo/doc.go @@ -0,0 +1,6 @@ +// Copyright (c) Tailscale Inc & AUTHORS +// SPDX-License-Identifier: BSD-3-Clause + +// Package geo provides functionality to represent and process geographical +// locations on a spherical Earth. +package geo diff --git a/types/geo/point.go b/types/geo/point.go new file mode 100644 index 000000000..d7160ac59 --- /dev/null +++ b/types/geo/point.go @@ -0,0 +1,279 @@ +// Copyright (c) Tailscale Inc & AUTHORS +// SPDX-License-Identifier: BSD-3-Clause + +package geo + +import ( + "encoding/binary" + "errors" + "fmt" + "math" + "strconv" +) + +// ErrBadPoint indicates that the point is malformed. +var ErrBadPoint = errors.New("not a valid point") + +// Point represents a pair of latitude and longitude coordinates. +type Point struct { + lat Degrees + // lng180 is the longitude offset by +180° so the zero value is invalid + // and +0+0/ is Point{lat: +0.0, lng180: +180.0}. + lng180 Degrees +} + +// MakePoint returns a Point representing a given latitude and longitude on +// a WGS 84 ellipsoid. The Coordinate Reference System is EPSG:4326. +// Latitude is wrapped to [-90°, +90°] and longitude to (-180°, +180°]. +func MakePoint(latitude, longitude Degrees) Point { + lat, lng := float64(latitude), float64(longitude) + + switch { + case math.IsNaN(lat) || math.IsInf(lat, 0): + // don’t wrap + case lat < -90 || lat > 90: + // Latitude wraps by flipping the longitude + lat = math.Mod(lat, 360.0) + switch { + case lat == 0.0: + lat = 0.0 // -0.0 == 0.0, but -0° is not valid + case lat < -270.0: + lat = +360.0 + lat + case lat < -90.0: + lat = -180.0 - lat + lng += 180.0 + case lat > +270.0: + lat = -360.0 + lat + case lat > +90.0: + lat = +180.0 - lat + lng += 180.0 + } + } + + switch { + case lat == -90.0 || lat == +90.0: + // By convention, the north and south poles have longitude 0°. + lng = 0 + case math.IsNaN(lng) || math.IsInf(lng, 0): + // don’t wrap + case lng <= -180.0 || lng > 180.0: + // Longitude wraps around normally + lng = math.Mod(lng, 360.0) + switch { + case lng == 0.0: + lng = 0.0 // -0.0 == 0.0, but -0° is not valid + case lng <= -180.0: + lng = +360.0 + lng + case lng > +180.0: + lng = -360.0 + lng + } + } + + return Point{ + lat: Degrees(lat), + lng180: Degrees(lng + 180.0), + } +} + +// Valid reports if p is a valid point. +func (p Point) Valid() bool { + return !p.IsZero() +} + +// LatLng reports the latitude and longitude. +func (p Point) LatLng() (lat, lng Degrees, err error) { + if p.IsZero() { + return 0 * Degree, 0 * Degree, ErrBadPoint + } + return p.lat, p.lng180 - 180.0*Degree, nil +} + +// LatLng reports the latitude and longitude in float64. If err is nil, then lat +// and lng will never both be 0.0 to disambiguate between an empty struct and +// Null Island (0° 0°). +func (p Point) LatLngFloat64() (lat, lng float64, err error) { + dlat, dlng, err := p.LatLng() + if err != nil { + return 0.0, 0.0, err + } + if dlat == 0.0 && dlng == 0.0 { + // dlng must survive conversion to float32. + dlng = math.SmallestNonzeroFloat32 + } + return float64(dlat), float64(dlng), err +} + +// SphericalAngleTo returns the angular distance from p to q, calculated on a +// spherical Earth. +func (p Point) SphericalAngleTo(q Point) (Radians, error) { + pLat, pLng, pErr := p.LatLng() + qLat, qLng, qErr := q.LatLng() + switch { + case pErr != nil && qErr != nil: + return 0.0, fmt.Errorf("spherical distance from %v to %v: %w", p, q, errors.Join(pErr, qErr)) + case pErr != nil: + return 0.0, fmt.Errorf("spherical distance from %v: %w", p, pErr) + case qErr != nil: + return 0.0, fmt.Errorf("spherical distance to %v: %w", q, qErr) + } + // The spherical law of cosines is accurate enough for close points when + // using float64. + // + // The haversine formula is an alternative, but it is poorly behaved + // when points are on opposite sides of the sphere. + rLat, rLng := float64(pLat.Radians()), float64(pLng.Radians()) + sLat, sLng := float64(qLat.Radians()), float64(qLng.Radians()) + cosA := math.Sin(rLat)*math.Sin(sLat) + + math.Cos(rLat)*math.Cos(sLat)*math.Cos(rLng-sLng) + return Radians(math.Acos(cosA)), nil +} + +// DistanceTo reports the great-circle distance between p and q, in meters. +func (p Point) DistanceTo(q Point) (Distance, error) { + r, err := p.SphericalAngleTo(q) + if err != nil { + return 0, err + } + return DistanceOnEarth(r.Turns()), nil +} + +// String returns a space-separated pair of latitude and longitude, in decimal +// degrees. Positive latitudes are in the northern hemisphere, and positive +// longitudes are east of the prime meridian. If p was not initialized, this +// will return "nowhere". +func (p Point) String() string { + lat, lng, err := p.LatLng() + if err != nil { + if err == ErrBadPoint { + return "nowhere" + } + panic(err) + } + + return lat.String() + " " + lng.String() +} + +// AppendBinary implements [encoding.BinaryAppender]. The output consists of two +// float32s in big-endian byte order: latitude and longitude offset by 180°. +// If p is not a valid, the output will be an 8-byte zero value. +func (p Point) AppendBinary(b []byte) ([]byte, error) { + end := binary.BigEndian + b = end.AppendUint32(b, math.Float32bits(float32(p.lat))) + b = end.AppendUint32(b, math.Float32bits(float32(p.lng180))) + return b, nil +} + +// MarshalBinary implements [encoding.BinaryMarshaller]. The output matches that +// of calling [Point.AppendBinary]. +func (p Point) MarshalBinary() ([]byte, error) { + var b [8]byte + return p.AppendBinary(b[:0]) +} + +// UnmarshalBinary implements [encoding.BinaryUnmarshaler]. It expects input +// that was formatted by [Point.AppendBinary]: in big-endian byte order, a +// float32 representing latitude followed by a float32 representing longitude +// offset by 180°. If latitude and longitude fall outside valid ranges, then +// an error is returned. +func (p *Point) UnmarshalBinary(data []byte) error { + if len(data) < 8 { // Two uint32s are 8 bytes long + return fmt.Errorf("%w: not enough data: %q", ErrBadPoint, data) + } + + end := binary.BigEndian + lat := Degrees(math.Float32frombits(end.Uint32(data[0:]))) + if lat < -90*Degree || lat > 90*Degree { + return fmt.Errorf("%w: latitude outside [-90°, +90°]: %s", ErrBadPoint, lat) + } + lng180 := Degrees(math.Float32frombits(end.Uint32(data[4:]))) + if lng180 != 0 && (lng180 < 0*Degree || lng180 > 360*Degree) { + // lng180 == 0 is OK: the zero value represents invalid points. + lng := lng180 - 180*Degree + return fmt.Errorf("%w: longitude outside (-180°, +180°]: %s", ErrBadPoint, lng) + } + + p.lat = lat + p.lng180 = lng180 + return nil +} + +// AppendText implements [encoding.TextAppender]. The output is a point +// formatted as OGC Well-Known Text, as "POINT (longitude latitude)" where +// longitude and latitude are in decimal degrees. If p is not valid, the output +// will be "POINT EMPTY". +func (p Point) AppendText(b []byte) ([]byte, error) { + if p.IsZero() { + b = append(b, []byte("POINT EMPTY")...) + return b, nil + } + + lat, lng, err := p.LatLng() + if err != nil { + return b, err + } + + b = append(b, []byte("POINT (")...) + b = strconv.AppendFloat(b, float64(lng), 'f', -1, 64) + b = append(b, ' ') + b = strconv.AppendFloat(b, float64(lat), 'f', -1, 64) + b = append(b, ')') + return b, nil +} + +// MarshalText implements [encoding.TextMarshaller]. The output matches that +// of calling [Point.AppendText]. +func (p Point) MarshalText() ([]byte, error) { + var b [8]byte + return p.AppendText(b[:0]) +} + +// MarshalUint64 produces the same output as MashalBinary, encoded in a uint64. +func (p Point) MarshalUint64() (uint64, error) { + b, err := p.MarshalBinary() + return binary.NativeEndian.Uint64(b), err +} + +// UnmarshalUint64 expects input formatted by MarshalUint64. +func (p *Point) UnmarshalUint64(v uint64) error { + b := binary.NativeEndian.AppendUint64(nil, v) + return p.UnmarshalBinary(b) +} + +// IsZero reports if p is the zero value. +func (p Point) IsZero() bool { + return p == Point{} +} + +// EqualApprox reports if p and q are approximately equal: that is the absolute +// difference of both latitude and longitude are less than tol. If tol is +// negative, then tol defaults to a reasonably small number (10⁻⁵). If tol is +// zero, then p and q must be exactly equal. +func (p Point) EqualApprox(q Point, tol float64) bool { + if tol == 0 { + return p == q + } + + if p.IsZero() && q.IsZero() { + return true + } else if p.IsZero() || q.IsZero() { + return false + } + + plat, plng, err := p.LatLng() + if err != nil { + panic(err) + } + qlat, qlng, err := q.LatLng() + if err != nil { + panic(err) + } + + if tol < 0 { + tol = 1e-5 + } + + dlat := float64(plat) - float64(qlat) + dlng := float64(plng) - float64(qlng) + return ((dlat < 0 && -dlat < tol) || (dlat >= 0 && dlat < tol)) && + ((dlng < 0 && -dlng < tol) || (dlng >= 0 && dlng < tol)) +} diff --git a/types/geo/point_test.go b/types/geo/point_test.go new file mode 100644 index 000000000..308c1a183 --- /dev/null +++ b/types/geo/point_test.go @@ -0,0 +1,541 @@ +// Copyright (c) Tailscale Inc & AUTHORS +// SPDX-License-Identifier: BSD-3-Clause + +package geo_test + +import ( + "fmt" + "math" + "testing" + "testing/quick" + + "tailscale.com/types/geo" +) + +func TestPointZero(t *testing.T) { + var zero geo.Point + + if got := zero.IsZero(); !got { + t.Errorf("IsZero() got %t", got) + } + + if got := zero.Valid(); got { + t.Errorf("Valid() got %t", got) + } + + wantErr := geo.ErrBadPoint.Error() + if _, _, err := zero.LatLng(); err.Error() != wantErr { + t.Errorf("LatLng() err %q, want %q", err, wantErr) + } + + wantStr := "nowhere" + if got := zero.String(); got != wantStr { + t.Errorf("String() got %q, want %q", got, wantStr) + } + + wantB := []byte{0, 0, 0, 0, 0, 0, 0, 0} + if b, err := zero.MarshalBinary(); err != nil { + t.Errorf("MarshalBinary() err %q, want nil", err) + } else if string(b) != string(wantB) { + t.Errorf("MarshalBinary got %q, want %q", b, wantB) + } + + wantI := uint64(0x00000000) + if i, err := zero.MarshalUint64(); err != nil { + t.Errorf("MarshalUint64() err %q, want nil", err) + } else if i != wantI { + t.Errorf("MarshalUint64 got %v, want %v", i, wantI) + } +} + +func TestPoint(t *testing.T) { + for _, tt := range []struct { + name string + lat geo.Degrees + lng geo.Degrees + wantLat geo.Degrees + wantLng geo.Degrees + wantString string + wantText string + }{ + { + name: "null-island", + lat: +0.0, + lng: +0.0, + wantLat: +0.0, + wantLng: +0.0, + wantString: "+0° +0°", + wantText: "POINT (0 0)", + }, + { + name: "north-pole", + lat: +90.0, + lng: +0.0, + wantLat: +90.0, + wantLng: +0.0, + wantString: "+90° +0°", + wantText: "POINT (0 90)", + }, + { + name: "south-pole", + lat: -90.0, + lng: +0.0, + wantLat: -90.0, + wantLng: +0.0, + wantString: "-90° +0°", + wantText: "POINT (0 -90)", + }, + { + name: "north-pole-weird-longitude", + lat: +90.0, + lng: +1.0, + wantLat: +90.0, + wantLng: +0.0, + wantString: "+90° +0°", + wantText: "POINT (0 90)", + }, + { + name: "south-pole-weird-longitude", + lat: -90.0, + lng: +1.0, + wantLat: -90.0, + wantLng: +0.0, + wantString: "-90° +0°", + wantText: "POINT (0 -90)", + }, + { + name: "almost-north", + lat: +89.0, + lng: +0.0, + wantLat: +89.0, + wantLng: +0.0, + wantString: "+89° +0°", + wantText: "POINT (0 89)", + }, + { + name: "past-north", + lat: +91.0, + lng: +0.0, + wantLat: +89.0, + wantLng: +180.0, + wantString: "+89° +180°", + wantText: "POINT (180 89)", + }, + { + name: "almost-south", + lat: -89.0, + lng: +0.0, + wantLat: -89.0, + wantLng: +0.0, + wantString: "-89° +0°", + wantText: "POINT (0 -89)", + }, + { + name: "past-south", + lat: -91.0, + lng: +0.0, + wantLat: -89.0, + wantLng: +180.0, + wantString: "-89° +180°", + wantText: "POINT (180 -89)", + }, + { + name: "antimeridian-north", + lat: +180.0, + lng: +0.0, + wantLat: +0.0, + wantLng: +180.0, + wantString: "+0° +180°", + wantText: "POINT (180 0)", + }, + { + name: "antimeridian-south", + lat: -180.0, + lng: +0.0, + wantLat: +0.0, + wantLng: +180.0, + wantString: "+0° +180°", + wantText: "POINT (180 0)", + }, + { + name: "almost-antimeridian-north", + lat: +179.0, + lng: +0.0, + wantLat: +1.0, + wantLng: +180.0, + wantString: "+1° +180°", + wantText: "POINT (180 1)", + }, + { + name: "past-antimeridian-north", + lat: +181.0, + lng: +0.0, + wantLat: -1.0, + wantLng: +180.0, + wantString: "-1° +180°", + wantText: "POINT (180 -1)", + }, + { + name: "almost-antimeridian-south", + lat: -179.0, + lng: +0.0, + wantLat: -1.0, + wantLng: +180.0, + wantString: "-1° +180°", + wantText: "POINT (180 -1)", + }, + { + name: "past-antimeridian-south", + lat: -181.0, + lng: +0.0, + wantLat: +1.0, + wantLng: +180.0, + wantString: "+1° +180°", + wantText: "POINT (180 1)", + }, + { + name: "circumnavigate-north", + lat: +360.0, + lng: +1.0, + wantLat: +0.0, + wantLng: +1.0, + wantString: "+0° +1°", + wantText: "POINT (1 0)", + }, + { + name: "circumnavigate-south", + lat: -360.0, + lng: +1.0, + wantLat: +0.0, + wantLng: +1.0, + wantString: "+0° +1°", + wantText: "POINT (1 0)", + }, + { + name: "almost-circumnavigate-north", + lat: +359.0, + lng: +1.0, + wantLat: -1.0, + wantLng: +1.0, + wantString: "-1° +1°", + wantText: "POINT (1 -1)", + }, + { + name: "past-circumnavigate-north", + lat: +361.0, + lng: +1.0, + wantLat: +1.0, + wantLng: +1.0, + wantString: "+1° +1°", + wantText: "POINT (1 1)", + }, + { + name: "almost-circumnavigate-south", + lat: -359.0, + lng: +1.0, + wantLat: +1.0, + wantLng: +1.0, + wantString: "+1° +1°", + wantText: "POINT (1 1)", + }, + { + name: "past-circumnavigate-south", + lat: -361.0, + lng: +1.0, + wantLat: -1.0, + wantLng: +1.0, + wantString: "-1° +1°", + wantText: "POINT (1 -1)", + }, + { + name: "antimeridian-east", + lat: +0.0, + lng: +180.0, + wantLat: +0.0, + wantLng: +180.0, + wantString: "+0° +180°", + wantText: "POINT (180 0)", + }, + { + name: "antimeridian-west", + lat: +0.0, + lng: -180.0, + wantLat: +0.0, + wantLng: +180.0, + wantString: "+0° +180°", + wantText: "POINT (180 0)", + }, + { + name: "almost-antimeridian-east", + lat: +0.0, + lng: +179.0, + wantLat: +0.0, + wantLng: +179.0, + wantString: "+0° +179°", + wantText: "POINT (179 0)", + }, + { + name: "past-antimeridian-east", + lat: +0.0, + lng: +181.0, + wantLat: +0.0, + wantLng: -179.0, + wantString: "+0° -179°", + wantText: "POINT (-179 0)", + }, + { + name: "almost-antimeridian-west", + lat: +0.0, + lng: -179.0, + wantLat: +0.0, + wantLng: -179.0, + wantString: "+0° -179°", + wantText: "POINT (-179 0)", + }, + { + name: "past-antimeridian-west", + lat: +0.0, + lng: -181.0, + wantLat: +0.0, + wantLng: +179.0, + wantString: "+0° +179°", + wantText: "POINT (179 0)", + }, + { + name: "montreal", + lat: +45.508888, + lng: -73.561668, + wantLat: +45.508888, + wantLng: -73.561668, + wantString: "+45.508888° -73.561668°", + wantText: "POINT (-73.561668 45.508888)", + }, + { + name: "canada", + lat: 57.550480044655636, + lng: -98.41680517868062, + wantLat: 57.550480044655636, + wantLng: -98.41680517868062, + wantString: "+57.550480044655636° -98.41680517868062°", + wantText: "POINT (-98.41680517868062 57.550480044655636)", + }, + } { + t.Run(tt.name, func(t *testing.T) { + p := geo.MakePoint(tt.lat, tt.lng) + + lat, lng, err := p.LatLng() + if !approx(lat, tt.wantLat) { + t.Errorf("MakePoint: lat %v, want %v", lat, tt.wantLat) + } + if !approx(lng, tt.wantLng) { + t.Errorf("MakePoint: lng %v, want %v", lng, tt.wantLng) + } + if err != nil { + t.Fatalf("LatLng: err %q, expected nil", err) + } + + if got := p.String(); got != tt.wantString { + t.Errorf("String: got %q, wantString %q", got, tt.wantString) + } + + txt, err := p.MarshalText() + if err != nil { + t.Errorf("Text: err %q, expected nil", err) + } else if string(txt) != tt.wantText { + t.Errorf("Text: got %q, wantText %q", txt, tt.wantText) + } + + b, err := p.MarshalBinary() + if err != nil { + t.Fatalf("MarshalBinary: err %q, expected nil", err) + } + + var q geo.Point + if err := q.UnmarshalBinary(b); err != nil { + t.Fatalf("UnmarshalBinary: err %q, expected nil", err) + } + if !q.EqualApprox(p, -1) { + t.Errorf("UnmarshalBinary: roundtrip failed: %#v != %#v", q, p) + } + + i, err := p.MarshalUint64() + if err != nil { + t.Fatalf("MarshalUint64: err %q, expected nil", err) + } + + var r geo.Point + if err := r.UnmarshalUint64(i); err != nil { + t.Fatalf("UnmarshalUint64: err %r, expected nil", err) + } + if !q.EqualApprox(r, -1) { + t.Errorf("UnmarshalUint64: roundtrip failed: %#v != %#v", r, p) + } + }) + } +} + +func TestPointMarshalBinary(t *testing.T) { + roundtrip := func(p geo.Point) error { + b, err := p.MarshalBinary() + if err != nil { + return fmt.Errorf("marshal: %v", err) + } + var q geo.Point + if err := q.UnmarshalBinary(b); err != nil { + return fmt.Errorf("unmarshal: %v", err) + } + if q != p { + return fmt.Errorf("%#v != %#v", q, p) + } + return nil + } + + t.Run("nowhere", func(t *testing.T) { + var nowhere geo.Point + if err := roundtrip(nowhere); err != nil { + t.Errorf("roundtrip: %v", err) + } + }) + + t.Run("quick-check", func(t *testing.T) { + f := func(lat geo.Degrees, lng geo.Degrees) (ok bool) { + pt := geo.MakePoint(lat, lng) + if err := roundtrip(pt); err != nil { + t.Errorf("roundtrip: %v", err) + } + return !t.Failed() + } + if err := quick.Check(f, nil); err != nil { + t.Error(err) + } + }) +} + +func TestPointMarshalUint64(t *testing.T) { + t.Skip("skip") + roundtrip := func(p geo.Point) error { + i, err := p.MarshalUint64() + if err != nil { + return fmt.Errorf("marshal: %v", err) + } + var q geo.Point + if err := q.UnmarshalUint64(i); err != nil { + return fmt.Errorf("unmarshal: %v", err) + } + if q != p { + return fmt.Errorf("%#v != %#v", q, p) + } + return nil + } + + t.Run("nowhere", func(t *testing.T) { + var nowhere geo.Point + if err := roundtrip(nowhere); err != nil { + t.Errorf("roundtrip: %v", err) + } + }) + + t.Run("quick-check", func(t *testing.T) { + f := func(lat geo.Degrees, lng geo.Degrees) (ok bool) { + if err := roundtrip(geo.MakePoint(lat, lng)); err != nil { + t.Errorf("roundtrip: %v", err) + } + return !t.Failed() + } + if err := quick.Check(f, nil); err != nil { + t.Error(err) + } + }) +} + +func TestPointSphericalAngleTo(t *testing.T) { + const earthRadius = 6371.000 // volumetric mean radius (km) + const kmToRad = 1 / earthRadius + for _, tt := range []struct { + name string + x geo.Point + y geo.Point + want geo.Radians + wantErr string + }{ + { + name: "same-point-null-island", + x: geo.MakePoint(0, 0), + y: geo.MakePoint(0, 0), + want: 0.0 * geo.Radian, + }, + { + name: "same-point-north-pole", + x: geo.MakePoint(+90, 0), + y: geo.MakePoint(+90, +90), + want: 0.0 * geo.Radian, + }, + { + name: "same-point-south-pole", + x: geo.MakePoint(-90, 0), + y: geo.MakePoint(-90, -90), + want: 0.0 * geo.Radian, + }, + { + name: "north-pole-to-south-pole", + x: geo.MakePoint(+90, 0), + y: geo.MakePoint(-90, -90), + want: math.Pi * geo.Radian, + }, + { + name: "toronto-to-montreal", + x: geo.MakePoint(+43.6532, -79.3832), + y: geo.MakePoint(+45.5019, -73.5674), + want: 504.26 * kmToRad * geo.Radian, + }, + { + name: "sydney-to-san-francisco", + x: geo.MakePoint(-33.8727, +151.2057), + y: geo.MakePoint(+37.7749, -122.4194), + want: 11948.18 * kmToRad * geo.Radian, + }, + { + name: "new-york-to-paris", + x: geo.MakePoint(+40.7128, -74.0060), + y: geo.MakePoint(+48.8575, +2.3514), + want: 5837.15 * kmToRad * geo.Radian, + }, + { + name: "seattle-to-tokyo", + x: geo.MakePoint(+47.6061, -122.3328), + y: geo.MakePoint(+35.6764, +139.6500), + want: 7700.00 * kmToRad * geo.Radian, + }, + } { + t.Run(tt.name, func(t *testing.T) { + got, err := tt.x.SphericalAngleTo(tt.y) + if tt.wantErr == "" && err != nil { + t.Fatalf("err %q, expected nil", err) + } + if tt.wantErr != "" && (err == nil || err.Error() != tt.wantErr) { + t.Fatalf("err %q, expected %q", err, tt.wantErr) + } + if tt.wantErr != "" { + return + } + + if !approx(got, tt.want) { + t.Errorf("x to y: got %v, want %v", got, tt.want) + } + + // Distance should be commutative + got, err = tt.y.SphericalAngleTo(tt.x) + if err != nil { + t.Fatalf("err %q, expected nil", err) + } + if !approx(got, tt.want) { + t.Errorf("y to x: got %v, want %v", got, tt.want) + } + t.Logf("x to y: %v km", got/kmToRad) + }) + } +} + +func approx[T ~float64](x, y T) bool { + return math.Abs(float64(x)-float64(y)) <= 1e-5 +} diff --git a/types/geo/quantize.go b/types/geo/quantize.go new file mode 100644 index 000000000..18ec11f9f --- /dev/null +++ b/types/geo/quantize.go @@ -0,0 +1,106 @@ +// Copyright (c) Tailscale Inc & AUTHORS +// SPDX-License-Identifier: BSD-3-Clause + +package geo + +import ( + "math" + "sync" +) + +// MinSeparation is the minimum separation between two points after quantizing. +// [Point.Quantize] guarantees that two points will either be snapped to exactly +// the same point, which conflates multiple positions together, or that the two +// points will be far enough apart that successfully performing most reverse +// lookups would be highly improbable. +const MinSeparation = 50_000 * Meter + +// Latitude +var ( + // numSepsEquatorToPole is the number of separations between a point on + // the equator to a point on a pole, that satisfies [minPointSep]. In + // other words, the number of separations between 0° and +90° degrees + // latitude. + numSepsEquatorToPole = int(math.Floor(float64( + earthPolarCircumference / MinSeparation / 4))) + + // latSep is the number of degrees between two adjacent latitudinal + // points. In other words, the next point going straight north of + // 0° would be latSep°. + latSep = Degrees(90.0 / float64(numSepsEquatorToPole)) +) + +// snapToLat returns the number of the nearest latitudinal separation to +// lat. A positive result is north of the equator, a negative result is south, +// and zero is the equator itself. For example, a result of -1 would mean a +// point that is [latSep] south of the equator. +func snapToLat(lat Degrees) int { + return int(math.Round(float64(lat / latSep))) +} + +// lngSep is a lookup table for the number of degrees between two adjacent +// longitudinal separations. where the index corresponds to the absolute value +// of the latitude separation. The first value corresponds to the equator and +// the last value corresponds to the separation before the pole. There is no +// value for the pole itself, because longitude has no meaning there. +// +// [lngSep] is calculated on init, which is so quick and will be used so often +// that the startup cost is negligible. +var lngSep = sync.OnceValue(func() []Degrees { + lut := make([]Degrees, numSepsEquatorToPole) + + // i ranges from the equator to a pole + for i := range len(lut) { + // lat ranges from [0°, 90°], because the southern hemisphere is + // a reflection of the northern one. + lat := Degrees(i) * latSep + ratio := math.Cos(float64(lat.Radians())) + circ := Distance(ratio) * earthEquatorialCircumference + num := int(math.Floor(float64(circ / MinSeparation))) + // We define lut[0] as 0°, lut[len(lut)] to be the north pole, + // which means -lut[len(lut)] is the south pole. + lut[i] = Degrees(360.0 / float64(num)) + } + return lut +}) + +// snapToLatLng returns the number of the nearest latitudinal separation to lat, +// and the nearest longitudinal separation to lng. +func snapToLatLng(lat, lng Degrees) (Degrees, Degrees) { + latN := snapToLat(lat) + + // absolute index into lngSep + n := latN + if n < 0 { + n = -latN + } + + lngSep := lngSep() + if n < len(lngSep) { + sep := lngSep[n] + lngN := int(math.Round(float64(lng / sep))) + return Degrees(latN) * latSep, Degrees(lngN) * sep + } + if latN < 0 { // south pole + return -90 * Degree, 0 * Degree + } else { // north pole + return +90 * Degree, 0 * Degree + } +} + +// Quantize returns a new [Point] after throwing away enough location data in p +// so that it would be difficult to distinguish a node among all the other nodes +// in its general vicinity. One caveat is that if there’s only one point in an +// obscure location, someone could triangulate the node using additional data. +// +// This method is stable: given the same p, it will always return the same +// result. It is equivalent to snapping to points on Earth that are at least +// [MinSeparation] apart. +func (p Point) Quantize() Point { + if p.IsZero() { + return p + } + + lat, lng := snapToLatLng(p.lat, p.lng180-180) + return MakePoint(lat, lng) +} diff --git a/types/geo/quantize_test.go b/types/geo/quantize_test.go new file mode 100644 index 000000000..3c707e303 --- /dev/null +++ b/types/geo/quantize_test.go @@ -0,0 +1,130 @@ +// Copyright (c) Tailscale Inc & AUTHORS +// SPDX-License-Identifier: BSD-3-Clause + +package geo_test + +import ( + "testing" + "testing/quick" + + "tailscale.com/types/geo" +) + +func TestPointAnonymize(t *testing.T) { + t.Run("nowhere", func(t *testing.T) { + var zero geo.Point + p := zero.Quantize() + want := zero.Valid() + if got := p.Valid(); got != want { + t.Fatalf("zero.Valid %t, want %t", got, want) + } + }) + + t.Run("separation", func(t *testing.T) { + // Walk from the south pole to the north pole and check that each + // point on the latitude is approximately MinSeparation apart. + const southPole = -90 * geo.Degree + const northPole = 90 * geo.Degree + const dateLine = 180 * geo.Degree + + llat := southPole + for lat := llat; lat <= northPole; lat += 0x1p-4 { + last := geo.MakePoint(llat, 0) + cur := geo.MakePoint(lat, 0) + anon := cur.Quantize() + switch l, g, err := anon.LatLng(); { + case err != nil: + t.Fatal(err) + case lat == southPole: + // initialize llng, to the first snapped longitude + llat = l + goto Lng + case g != 0: + t.Fatalf("%v is west or east of %v", anon, last) + case l < llat: + t.Fatalf("%v is south of %v", anon, last) + case l == llat: + continue + case l > llat: + switch dist, err := last.DistanceTo(anon); { + case err != nil: + t.Fatal(err) + case dist == 0.0: + continue + case dist < geo.MinSeparation: + t.Logf("lat=%v last=%v cur=%v anon=%v", lat, last, cur, anon) + t.Fatalf("%v is too close to %v", anon, last) + default: + llat = l + } + } + + Lng: + llng := dateLine + for lng := llng; lng <= dateLine && lng >= -dateLine; lng -= 0x1p-3 { + last := geo.MakePoint(llat, llng) + cur := geo.MakePoint(lat, lng) + anon := cur.Quantize() + switch l, g, err := anon.LatLng(); { + case err != nil: + t.Fatal(err) + case lng == dateLine: + // initialize llng, to the first snapped longitude + llng = g + continue + case l != llat: + t.Fatalf("%v is north or south of %v", anon, last) + case g != llng: + const tolerance = geo.MinSeparation * 0x1p-9 + switch dist, err := last.DistanceTo(anon); { + case err != nil: + t.Fatal(err) + case dist < tolerance: + continue + case dist < (geo.MinSeparation - tolerance): + t.Logf("lat=%v lng=%v last=%v cur=%v anon=%v", lat, lng, last, cur, anon) + t.Fatalf("%v is too close to %v: %v", anon, last, dist) + default: + llng = g + } + + } + } + } + if llat == southPole { + t.Fatal("llat never incremented") + } + }) + + t.Run("quick-check", func(t *testing.T) { + f := func(lat, lng geo.Degrees) bool { + p := geo.MakePoint(lat, lng) + q := p.Quantize() + t.Logf("quantize %v = %v", p, q) + + lat, lng, err := q.LatLng() + if err != nil { + t.Errorf("err %v, want nil", err) + return !t.Failed() + } + + if lat < -90*geo.Degree || lat > 90*geo.Degree { + t.Errorf("lat outside [-90°, +90°]: %v", lat) + } + if lng < -180*geo.Degree || lng > 180*geo.Degree { + t.Errorf("lng outside [-180°, +180°], %v", lng) + } + + if dist, err := p.DistanceTo(q); err != nil { + t.Error(err) + } else if dist > (geo.MinSeparation * 2) { + t.Errorf("moved too far: %v", dist) + } + + return !t.Failed() + } + if err := quick.Check(f, nil); err != nil { + t.Fatal(err) + } + }) +} diff --git a/types/geo/units.go b/types/geo/units.go new file mode 100644 index 000000000..76a4c02f7 --- /dev/null +++ b/types/geo/units.go @@ -0,0 +1,191 @@ +// Copyright (c) Tailscale Inc & AUTHORS +// SPDX-License-Identifier: BSD-3-Clause + +package geo + +import ( + "math" + "strconv" + "strings" + "unicode" +) + +const ( + Degree Degrees = 1 + Radian Radians = 1 + Turn Turns = 1 + Meter Distance = 1 +) + +// Degrees represents a latitude or longitude, in decimal degrees. +type Degrees float64 + +// ParseDegrees parses s as decimal degrees. +func ParseDegrees(s string) (Degrees, error) { + s = strings.TrimSuffix(s, "°") + f, err := strconv.ParseFloat(s, 64) + return Degrees(f), err +} + +// MustParseDegrees parses s as decimal degrees, but panics on error. +func MustParseDegrees(s string) Degrees { + d, err := ParseDegrees(s) + if err != nil { + panic(err) + } + return d +} + +// String implements the [Stringer] interface. The output is formatted in +// decimal degrees, prefixed by either the appropriate + or - sign, and suffixed +// by a ° degree symbol. +func (d Degrees) String() string { + b, _ := d.AppendText(nil) + b = append(b, []byte("°")...) + return string(b) +} + +// AppendText implements [encoding.TextAppender]. The output is formatted in +// decimal degrees, prefixed by either the appropriate + or - sign. +func (d Degrees) AppendText(b []byte) ([]byte, error) { + b = d.AppendZeroPaddedText(b, 0) + return b, nil +} + +// AppendZeroPaddedText appends d formatted as decimal degrees to b. The number of +// integer digits will be zero-padded to nint. +func (d Degrees) AppendZeroPaddedText(b []byte, nint int) []byte { + n := float64(d) + + if math.IsInf(n, 0) || math.IsNaN(n) { + return strconv.AppendFloat(b, n, 'f', -1, 64) + } + + sign := byte('+') + if math.Signbit(n) { + sign = '-' + n = -n + } + b = append(b, sign) + + pad := nint - 1 + for nn := n / 10; nn >= 1 && pad > 0; nn /= 10 { + pad-- + } + for range pad { + b = append(b, '0') + } + return strconv.AppendFloat(b, n, 'f', -1, 64) +} + +// Radians converts d into radians. +func (d Degrees) Radians() Radians { + return Radians(d * math.Pi / 180.0) +} + +// Turns converts d into a number of turns. +func (d Degrees) Turns() Turns { + return Turns(d / 360.0) +} + +// Radians represents a latitude or longitude, in radians. +type Radians float64 + +// ParseRadians parses s as radians. +func ParseRadians(s string) (Radians, error) { + s = strings.TrimSuffix(s, "rad") + s = strings.TrimRightFunc(s, unicode.IsSpace) + f, err := strconv.ParseFloat(s, 64) + return Radians(f), err +} + +// MustParseRadians parses s as radians, but panics on error. +func MustParseRadians(s string) Radians { + r, err := ParseRadians(s) + if err != nil { + panic(err) + } + return r +} + +// String implements the [Stringer] interface. +func (r Radians) String() string { + return strconv.FormatFloat(float64(r), 'f', -1, 64) + " rad" +} + +// Degrees converts r into decimal degrees. +func (r Radians) Degrees() Degrees { + return Degrees(r * 180.0 / math.Pi) +} + +// Turns converts r into a number of turns. +func (r Radians) Turns() Turns { + return Turns(r / 2 / math.Pi) +} + +// Turns represents a number of complete revolutions around a sphere. +type Turns float64 + +// String implements the [Stringer] interface. +func (o Turns) String() string { + return strconv.FormatFloat(float64(o), 'f', -1, 64) +} + +// Degrees converts t into decimal degrees. +func (o Turns) Degrees() Degrees { + return Degrees(o * 360.0) +} + +// Radians converts t into radians. +func (o Turns) Radians() Radians { + return Radians(o * 2 * math.Pi) +} + +// Distance represents a great-circle distance in meters. +type Distance float64 + +// ParseDistance parses s as distance in meters. +func ParseDistance(s string) (Distance, error) { + s = strings.TrimSuffix(s, "m") + s = strings.TrimRightFunc(s, unicode.IsSpace) + f, err := strconv.ParseFloat(s, 64) + return Distance(f), err +} + +// MustParseDistance parses s as distance in meters, but panics on error. +func MustParseDistance(s string) Distance { + d, err := ParseDistance(s) + if err != nil { + panic(err) + } + return d +} + +// String implements the [Stringer] interface. +func (d Distance) String() string { + return strconv.FormatFloat(float64(d), 'f', -1, 64) + "m" +} + +// DistanceOnEarth converts t turns into the great-circle distance, in meters. +func DistanceOnEarth(t Turns) Distance { + return Distance(t) * EarthMeanCircumference +} + +// Earth Fact Sheet +// https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html +const ( + // EarthMeanRadius is the volumetric mean radius of the Earth. + EarthMeanRadius = 6_371_000 * Meter + // EarthMeanCircumference is the volumetric mean circumference of the Earth. + EarthMeanCircumference = 2 * math.Pi * EarthMeanRadius + + // earthEquatorialRadius is the equatorial radius of the Earth. + earthEquatorialRadius = 6_378_137 * Meter + // earthEquatorialCircumference is the equatorial circumference of the Earth. + earthEquatorialCircumference = 2 * math.Pi * earthEquatorialRadius + + // earthPolarRadius is the polar radius of the Earth. + earthPolarRadius = 6_356_752 * Meter + // earthPolarCircumference is the polar circumference of the Earth. + earthPolarCircumference = 2 * math.Pi * earthPolarRadius +) diff --git a/types/geo/units_test.go b/types/geo/units_test.go new file mode 100644 index 000000000..b6f724ce0 --- /dev/null +++ b/types/geo/units_test.go @@ -0,0 +1,395 @@ +// Copyright (c) Tailscale Inc & AUTHORS +// SPDX-License-Identifier: BSD-3-Clause + +package geo_test + +import ( + "math" + "strings" + "testing" + + "tailscale.com/types/geo" +) + +func TestDegrees(t *testing.T) { + for _, tt := range []struct { + name string + degs geo.Degrees + wantStr string + wantText string + wantPad string + wantRads geo.Radians + wantTurns geo.Turns + }{ + { + name: "zero", + degs: 0.0 * geo.Degree, + wantStr: "+0°", + wantText: "+0", + wantPad: "+000", + wantRads: 0.0 * geo.Radian, + wantTurns: 0 * geo.Turn, + }, + { + name: "quarter-turn", + degs: 90.0 * geo.Degree, + wantStr: "+90°", + wantText: "+90", + wantPad: "+090", + wantRads: 0.5 * math.Pi * geo.Radian, + wantTurns: 0.25 * geo.Turn, + }, + { + name: "half-turn", + degs: 180.0 * geo.Degree, + wantStr: "+180°", + wantText: "+180", + wantPad: "+180", + wantRads: 1.0 * math.Pi * geo.Radian, + wantTurns: 0.5 * geo.Turn, + }, + { + name: "full-turn", + degs: 360.0 * geo.Degree, + wantStr: "+360°", + wantText: "+360", + wantPad: "+360", + wantRads: 2.0 * math.Pi * geo.Radian, + wantTurns: 1.0 * geo.Turn, + }, + { + name: "negative-zero", + degs: geo.MustParseDegrees("-0.0"), + wantStr: "-0°", + wantText: "-0", + wantPad: "-000", + wantRads: 0 * geo.Radian * -1, + wantTurns: 0 * geo.Turn * -1, + }, + { + name: "small-degree", + degs: -1.2003 * geo.Degree, + wantStr: "-1.2003°", + wantText: "-1.2003", + wantPad: "-001.2003", + wantRads: -0.020949187011687936 * geo.Radian, + wantTurns: -0.0033341666666666663 * geo.Turn, + }, + } { + t.Run(tt.name, func(t *testing.T) { + if got := tt.degs.String(); got != tt.wantStr { + t.Errorf("String got %q, want %q", got, tt.wantStr) + } + + d, err := geo.ParseDegrees(tt.wantStr) + if err != nil { + t.Fatalf("ParseDegrees err %q, want nil", err.Error()) + } + if d != tt.degs { + t.Errorf("ParseDegrees got %q, want %q", d, tt.degs) + } + + b, err := tt.degs.AppendText(nil) + if err != nil { + t.Fatalf("AppendText err %q, want nil", err.Error()) + } + if string(b) != tt.wantText { + t.Errorf("AppendText got %q, want %q", b, tt.wantText) + } + + b = tt.degs.AppendZeroPaddedText(nil, 3) + if string(b) != tt.wantPad { + t.Errorf("AppendZeroPaddedText got %q, want %q", b, tt.wantPad) + } + + r := tt.degs.Radians() + if r != tt.wantRads { + t.Errorf("Radian got %v, want %v", r, tt.wantRads) + } + if d := r.Degrees(); d != tt.degs { // Roundtrip + t.Errorf("Degrees got %v, want %v", d, tt.degs) + } + + o := tt.degs.Turns() + if o != tt.wantTurns { + t.Errorf("Turns got %v, want %v", o, tt.wantTurns) + } + }) + } +} + +func TestRadians(t *testing.T) { + for _, tt := range []struct { + name string + rads geo.Radians + wantStr string + wantText string + wantDegs geo.Degrees + wantTurns geo.Turns + }{ + { + name: "zero", + rads: 0.0 * geo.Radian, + wantStr: "0 rad", + wantDegs: 0.0 * geo.Degree, + wantTurns: 0 * geo.Turn, + }, + { + name: "quarter-turn", + rads: 0.5 * math.Pi * geo.Radian, + wantStr: "1.5707963267948966 rad", + wantDegs: 90.0 * geo.Degree, + wantTurns: 0.25 * geo.Turn, + }, + { + name: "half-turn", + rads: 1.0 * math.Pi * geo.Radian, + wantStr: "3.141592653589793 rad", + wantDegs: 180.0 * geo.Degree, + wantTurns: 0.5 * geo.Turn, + }, + { + name: "full-turn", + rads: 2.0 * math.Pi * geo.Radian, + wantStr: "6.283185307179586 rad", + wantDegs: 360.0 * geo.Degree, + wantTurns: 1.0 * geo.Turn, + }, + { + name: "negative-zero", + rads: geo.MustParseRadians("-0"), + wantStr: "-0 rad", + wantDegs: 0 * geo.Degree * -1, + wantTurns: 0 * geo.Turn * -1, + }, + } { + t.Run(tt.name, func(t *testing.T) { + if got := tt.rads.String(); got != tt.wantStr { + t.Errorf("String got %q, want %q", got, tt.wantStr) + } + + r, err := geo.ParseRadians(tt.wantStr) + if err != nil { + t.Fatalf("ParseDegrees err %q, want nil", err.Error()) + } + if r != tt.rads { + t.Errorf("ParseDegrees got %q, want %q", r, tt.rads) + } + + d := tt.rads.Degrees() + if d != tt.wantDegs { + t.Errorf("Degrees got %v, want %v", d, tt.wantDegs) + } + if r := d.Radians(); r != tt.rads { // Roundtrip + t.Errorf("Radians got %v, want %v", r, tt.rads) + } + + o := tt.rads.Turns() + if o != tt.wantTurns { + t.Errorf("Turns got %v, want %v", o, tt.wantTurns) + } + }) + } +} + +func TestTurns(t *testing.T) { + for _, tt := range []struct { + name string + turns geo.Turns + wantStr string + wantText string + wantDegs geo.Degrees + wantRads geo.Radians + }{ + { + name: "zero", + turns: 0.0, + wantStr: "0", + wantDegs: 0.0 * geo.Degree, + wantRads: 0 * geo.Radian, + }, + { + name: "quarter-turn", + turns: 0.25, + wantStr: "0.25", + wantDegs: 90.0 * geo.Degree, + wantRads: 0.5 * math.Pi * geo.Radian, + }, + { + name: "half-turn", + turns: 0.5, + wantStr: "0.5", + wantDegs: 180.0 * geo.Degree, + wantRads: 1.0 * math.Pi * geo.Radian, + }, + { + name: "full-turn", + turns: 1.0, + wantStr: "1", + wantDegs: 360.0 * geo.Degree, + wantRads: 2.0 * math.Pi * geo.Radian, + }, + { + name: "negative-zero", + turns: geo.Turns(math.Copysign(0, -1)), + wantStr: "-0", + wantDegs: 0 * geo.Degree * -1, + wantRads: 0 * geo.Radian * -1, + }, + } { + t.Run(tt.name, func(t *testing.T) { + if got := tt.turns.String(); got != tt.wantStr { + t.Errorf("String got %q, want %q", got, tt.wantStr) + } + + d := tt.turns.Degrees() + if d != tt.wantDegs { + t.Errorf("Degrees got %v, want %v", d, tt.wantDegs) + } + if o := d.Turns(); o != tt.turns { // Roundtrip + t.Errorf("Turns got %v, want %v", o, tt.turns) + } + + r := tt.turns.Radians() + if r != tt.wantRads { + t.Errorf("Turns got %v, want %v", r, tt.wantRads) + } + }) + } +} + +func TestDistance(t *testing.T) { + for _, tt := range []struct { + name string + dist geo.Distance + wantStr string + }{ + { + name: "zero", + dist: 0.0 * geo.Meter, + wantStr: "0m", + }, + { + name: "random", + dist: 4 * geo.Meter, + wantStr: "4m", + }, + { + name: "light-second", + dist: 299_792_458 * geo.Meter, + wantStr: "299792458m", + }, + { + name: "planck-length", + dist: 1.61625518e-35 * geo.Meter, + wantStr: "0.0000000000000000000000000000000000161625518m", + }, + { + name: "negative-zero", + dist: geo.Distance(math.Copysign(0, -1)), + wantStr: "-0m", + }, + } { + t.Run(tt.name, func(t *testing.T) { + if got := tt.dist.String(); got != tt.wantStr { + t.Errorf("String got %q, want %q", got, tt.wantStr) + } + + r, err := geo.ParseDistance(tt.wantStr) + if err != nil { + t.Fatalf("ParseDegrees err %q, want nil", err.Error()) + } + if r != tt.dist { + t.Errorf("ParseDegrees got %q, want %q", r, tt.dist) + } + }) + } +} + +func TestDistanceOnEarth(t *testing.T) { + for _, tt := range []struct { + name string + here geo.Point + there geo.Point + want geo.Distance + wantErr string + }{ + { + name: "no-points", + here: geo.Point{}, + there: geo.Point{}, + wantErr: "not a valid point", + }, + { + name: "not-here", + here: geo.Point{}, + there: geo.MakePoint(0, 0), + wantErr: "not a valid point", + }, + { + name: "not-there", + here: geo.MakePoint(0, 0), + there: geo.Point{}, + wantErr: "not a valid point", + }, + { + name: "null-island", + here: geo.MakePoint(0, 0), + there: geo.MakePoint(0, 0), + want: 0 * geo.Meter, + }, + { + name: "equator-to-south-pole", + here: geo.MakePoint(0, 0), + there: geo.MakePoint(-90, 0), + want: geo.EarthMeanCircumference / 4, + }, + { + name: "north-pole-to-south-pole", + here: geo.MakePoint(+90, 0), + there: geo.MakePoint(-90, 0), + want: geo.EarthMeanCircumference / 2, + }, + { + name: "meridian-to-antimeridian", + here: geo.MakePoint(0, 0), + there: geo.MakePoint(0, -180), + want: geo.EarthMeanCircumference / 2, + }, + { + name: "positive-to-negative-antimeridian", + here: geo.MakePoint(0, 180), + there: geo.MakePoint(0, -180), + want: 0 * geo.Meter, + }, + { + name: "toronto-to-montreal", + here: geo.MakePoint(+43.70011, -79.41630), + there: geo.MakePoint(+45.50884, -73.58781), + want: 503_200 * geo.Meter, + }, + { + name: "montreal-to-san-francisco", + here: geo.MakePoint(+45.50884, -73.58781), + there: geo.MakePoint(+37.77493, -122.41942), + want: 4_082_600 * geo.Meter, + }, + } { + t.Run(tt.name, func(t *testing.T) { + got, err := tt.here.DistanceTo(tt.there) + if tt.wantErr == "" && err != nil { + t.Fatalf("err %q, want nil", err) + } + if tt.wantErr != "" && !strings.Contains(err.Error(), tt.wantErr) { + t.Fatalf("err %q, want %q", err, tt.wantErr) + } + + approx := func(x, y geo.Distance) bool { + return math.Abs(float64(x)-float64(y)) <= 10 + } + if !approx(got, tt.want) { + t.Fatalf("got %v, want %v", got, tt.want) + } + }) + } +}