example: ProjectedCoords #43 with new proposed coords system

docs/src/api.md

@@ -22,6 +22,7 @@ FK4NoETerms
 FK5Coords
 EclipticCoords
 CartesianCoords
+ProjectedCoords
 ```
 
 ## Conversion
@@ -65,6 +66,7 @@ inrange
 separation
 position_angle
 offset
+SkyCoords.project
 cartesian
 spherical
 ```

src/SkyCoords.jl

@@ -14,6 +14,7 @@ export AbstractSkyCoords,
     FK5Coords,
     EclipticCoords,
     CartesianCoords,
+    ProjectedCoords,
     separation,
     position_angle,
     offset,
@@ -22,6 +23,7 @@ export AbstractSkyCoords,
 
 include("types.jl")
 include("cartesian.jl")
+include("projected.jl")
 
 # -----------------------------------------------------------------------------
 # Helper functions: Create rotation matrix about a given axis (x, y, z)

src/projected.jl

@@ -0,0 +1,79 @@
+abstract type AbstractProjectedCoords <: AbstractSkyCoords end
+
+"""
+    ProjectedCoords{TC <: AbstractSkyCoords, T <: Real} <: AbstractSkyCoords
+
+Represent coordinate `c` as a flat-sky offset from `origin` using a
+small-field-of-view approximation.
+
+`offset[1]` is the longitude difference scaled by `cos(lat(origin))`, and
+`offset[2]` is the latitude difference, both in radians.
+
+A projected coordinate is a _representation_ around a point, not a coordinate frame of its own.
+Its frame is the origin's. Conversions out of a `ProjectedCoords` (spherical or [`CartesianCoords`](@ref) targets alike)
+work like any other coordinate. Conversions into one are not possible via `convert`, since they require an origin value.
+Use [`project`](@ref) instead.
+"""
+struct ProjectedCoords{TC <: AbstractSkyCoords, T <: Real} <: AbstractProjectedCoords
+    origin::TC
+    offset::SVector{2, T}
+
+    # Typed like the spherical types' inner constructors so the two-argument
+    # form stays disjoint from extension constructors taking other argument
+    # pairs
+    ProjectedCoords{TC, T}(origin::AbstractSkyCoords, offset::AbstractVector{<:Real}) where {TC <: AbstractSkyCoords, T <: Real} =
+        new(origin, offset)
+end
+ProjectedCoords(origin::TC, offset::SVector{2, T}) where {TC <: AbstractSkyCoords, T <: Real} =
+    ProjectedCoords{TC, T}(origin, offset)
+
+origin(c::ProjectedCoords) = c.origin
+lon(c::AbstractProjectedCoords) = lon(origin(c)) + c.offset[1] / cos(lat(origin(c)))
+lat(c::AbstractProjectedCoords) = lat(origin(c)) + c.offset[2]
+
+"""
+    project(origin::AbstractSkyCoords, c::AbstractSkyCoords) -> ProjectedCoords
+
+Represent `c` as a [`ProjectedCoords`](@ref) offset from `origin`.
+
+`c` is first converted to `origin`'s frame, so the offset is expressed in that frame.
+"""
+function project(origin::AbstractSkyCoords, c::AbstractSkyCoords)
+    cc = convert(typeof(origin), c)
+    Δlon = rem2pi(lon(cc) - lon(origin), RoundNearest)
+    offset = SVector(Δlon * cos(lat(origin)), lat(cc) - lat(origin))
+    return ProjectedCoords(origin, offset)
+end
+
+# A projected coordinate's frame is its origin's, so the frame-change
+# primitive delegates to it; every `convert` out of a ProjectedCoords then
+# works generically (`coords2cart` sees the projected lon/lat above).
+frame_transform(::Type{T}, ::Type{<:ProjectedCoords{TC}}, v) where {T <: AbstractSkyCoords, TC <: AbstractSkyCoords} =
+    frame_transform(T, TC, v)
+# Disambiguate against FK4Coords's own generic-source frame_transform method
+frame_transform(::Type{<:FK4Coords{e}}, ::Type{<:ProjectedCoords{TC}}, v) where {e, TC <: AbstractSkyCoords} =
+    frame_transform(FK4Coords{e}, TC, v)
+
+# A ProjectedCoords target can never be satisfied by `convert`. The origin
+# is a value, not a type parameter. Every conversion pathway checks its
+# target frame through `_checkframe`, so one override turns them all into a
+# clear error (instances already matching the target still convert by identity without reaching this).
+_checkframe(::Type{TC}) where {TC <: AbstractProjectedCoords} = throw(ArgumentError(
+    "Cannot `convert` into a projected coordinate type; construct one with " *
+    "`project(origin, c)` instead.",
+))
+
+# The Cartesian representation is tagged with the origin's frame.
+# `constructorof(typeof(c))` (the generic tag choice) must keep returning
+# `ProjectedCoords` for `setproperties` reconstruction, but the frame this vector lives in is the origin's.
+cartesian(c::AbstractProjectedCoords) = CartesianCoords{constructorof(typeof(origin(c)))}(coords2cart(c))
+
+# Equality cannot fall back to the generic frame-tag + lon/lat comparison in
+# types.jl: the bare `ProjectedCoords` tag erases the origin's frame, so raw
+# lon/lat numbers from different origin frames would compare equal. Compare
+# the origin (frame included) and the offset instead, with a matching `hash`.
+Base.:(==)(a::AbstractProjectedCoords, b::AbstractProjectedCoords) =
+    origin(a) == origin(b) && a.offset == b.offset
+Base.hash(c::AbstractProjectedCoords, h::UInt) = hash(c.offset, hash(origin(c), h))
+Base.isapprox(a::ProjectedCoords, b::ProjectedCoords; kwargs...) =
+    isapprox(origin(a), origin(b); kwargs...) && isapprox(a.offset, b.offset; kwargs...)

src/types.jl

@@ -188,10 +188,19 @@ end
 # `Dict`/`Set` as required by the `hash` contract.
 Base.:(==)(a::AbstractSkyCoords, b::AbstractSkyCoords) = constructorof(typeof(a)) == constructorof(typeof(b)) && lonlat(a) == lonlat(b)
 Base.hash(c::AbstractSkyCoords, h::UInt) = hash(lonlat(c), hash(constructorof(typeof(c)), h))
-Base.isapprox(a::ICRSCoords, b::ICRSCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
-Base.isapprox(a::GalCoords, b::GalCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
-Base.isapprox(a::SuperGalCoords, b::SuperGalCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
-Base.isapprox(a::FK4Coords{e}, b::FK4Coords{e}; kwargs...) where {e} = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
-Base.isapprox(a::FK4NoETerms{e}, b::FK4NoETerms{e}; kwargs...) where {e} = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
-Base.isapprox(a::FK5Coords{e}, b::FK5Coords{e}; kwargs...) where {e} = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
-Base.isapprox(a::EclipticCoords{e}, b::EclipticCoords{e}; kwargs...) where {e} = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
+
+# Same-frame approximate equality compares (lon, lat), with the longitude
+# difference taken in (-π, π] so that nearly-equal points on either side of
+# the lon = 0 wrap still compare ≈, e.g., `ICRSCoords(eps(), 1) ≈ ICRSCoords(-eps(), 1)`.
+_isapprox_lonlat(a, b; kwargs...) = isapprox(
+    SVector(lon(a), lat(a)),
+    SVector(lon(a) + rem2pi(lon(b) - lon(a), RoundNearest), lat(b));
+    kwargs...
+)
+Base.isapprox(a::ICRSCoords, b::ICRSCoords; kwargs...) = _isapprox_lonlat(a, b; kwargs...)
+Base.isapprox(a::GalCoords, b::GalCoords; kwargs...) = _isapprox_lonlat(a, b; kwargs...)
+Base.isapprox(a::SuperGalCoords, b::SuperGalCoords; kwargs...) = _isapprox_lonlat(a, b; kwargs...)
+Base.isapprox(a::FK4Coords{e}, b::FK4Coords{e}; kwargs...) where {e} = _isapprox_lonlat(a, b; kwargs...)
+Base.isapprox(a::FK4NoETerms{e}, b::FK4NoETerms{e}; kwargs...) where {e} = _isapprox_lonlat(a, b; kwargs...)
+Base.isapprox(a::FK5Coords{e}, b::FK5Coords{e}; kwargs...) where {e} = _isapprox_lonlat(a, b; kwargs...)
+Base.isapprox(a::EclipticCoords{e}, b::EclipticCoords{e}; kwargs...) where {e} = _isapprox_lonlat(a, b; kwargs...)

test/runtests.jl

@@ -14,13 +14,65 @@ using Test
 import Makie
 
 import SkyCoords: lat, lon
+using SkyCoords: project, origin
 
 const rng = StableRNG(2000)
 rad2arcsec(r) = 3600 * rad2deg(r)
 
 # tests against astropy.coordinates
 include("astropy.jl")
 
+@testset "projected coords" begin
+    c0 = ICRSCoords(0.1, -0.2)
+    c1 = ICRSCoords(0.1 + 1.0e-5, -0.2 + 3.0e-5)
+    cp = project(c0, c1)::ProjectedCoords
+    @test origin(cp) == c0
+    @test cp.offset[1] ≈ 0.98 * 1.0e-5  rtol = 1.0e-4
+    @test cp.offset[2] ≈ 3.0e-5
+    @test convert(ICRSCoords, cp) ≈ c1
+    @test convert(GalCoords, cp) ≈ convert(GalCoords, c1)
+    @test cp == cp
+    @test cp ≈ cp
+
+    # Equality compares origin (frame included) and offset. The same raw
+    # lon/lat numbers around an origin in a different frame are a different point on the sky
+    @test cp == ProjectedCoords(c0, cp.offset)
+    @test hash(cp) == hash(ProjectedCoords(c0, cp.offset))
+    @test cp != ProjectedCoords(GalCoords(0.1, -0.2), cp.offset)
+    @test cp != ProjectedCoords(c0, zero(cp.offset))
+    @test cp != c1
+
+    # Same-type separation works through the projected lon/lat directly
+    @test separation(cp, project(c0, c0)) ≈ separation(c1, c0)
+
+    # Projected coords are first-class in the conversion lattice.
+    # Cartesian targets work through `frame_transform` and tag with the origin's frame.
+    cart = @inferred cartesian(cp)
+    @test cart isa CartesianCoords{ICRSCoords, Float64}
+    @test cart ≈ cartesian(convert(ICRSCoords, cp))
+    @test convert(CartesianCoords, cp) == cart
+    @test convert(CartesianCoords{GalCoords}, cp) ≈ cartesian(convert(GalCoords, c1))
+    @test convert(CartesianCoords{GalCoords, Float32}, cp) isa CartesianCoords{GalCoords, Float32}
+    @test spherical(cart) isa ICRSCoords
+
+    # A non-rotational origin frame works too (E-terms via frame_transform)
+    f0 = FK4Coords{1950}(0.3, 0.1)
+    f1 = offset(f0, 1.0e-4, 0.5)
+    fp = project(f0, f1)
+    @test convert(FK4Coords{1950}, fp) ≈ f1
+    @test convert(ICRSCoords, fp) ≈ convert(ICRSCoords, f1)
+    @test convert(CartesianCoords{ICRSCoords}, fp) ≈ cartesian(convert(ICRSCoords, f1))
+
+    # Nested projections chain through their origins
+    @test convert(ICRSCoords, ProjectedCoords(cp, zero(cp.offset))) ≈ c1
+
+    # Converting into a projected type needs an origin value, which a type
+    # cannot carry. Instances already matching the target stay identity.
+    @test_throws ArgumentError convert(typeof(cp), c1)
+    @test_throws ArgumentError convert(typeof(cp), cartesian(cp))
+    @test convert(typeof(cp), cp) === cp
+end
+
 # Test separation between coordinates and conversion with mixed floating types.
 @testset "Separation" begin
     c1 = ICRSCoords(ℯ, pi / 2)
@@ -408,6 +460,13 @@ end
         @test c1 ≈ c2  rtol = 1.0e-3
         @test c1 ≈ c4  rtol = 1.0e-3
 
+         # longitude comparison is periodic: points on either side of the
+         # lon = 0 wrap are still ≈
+         @test T(eps(), 1) ≈ T(0, 1)
+         @test T(eps(), 1) ≈ T(-eps(), 1)
+         @test !(T(π, 1) ≈ T(-π + 1.0e-3, 1))
+         @test T(π, 1) ≈ T(-π + 1.0e-3, 1)  atol = 1.0e-2
+
         # `==` implies equal hashes, so value-equal coordinates of different
         # element types collapse in a Set; c2 and c4 stay distinct because
         # 2.001 rounds to different values in Float32 and Float64