// 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)) }