From 243c203592c2ca2122a5979c0593cc494e7d2a5d Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Fri, 17 Jul 2026 03:38:54 -0700 Subject: [PATCH 1/3] fix: Make same-frame isapprox periodic in longitude, apply to all coord types --- src/types.jl | 23 ++++++++++++++++------- test/runtests.jl | 7 +++++++ 2 files changed, 23 insertions(+), 7 deletions(-) diff --git a/src/types.jl b/src/types.jl index 0785432..c689868 100644 --- a/src/types.jl +++ b/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...) diff --git a/test/runtests.jl b/test/runtests.jl index 0a14535..8e8fbf5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -375,6 +375,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 From 0032dc3898b6d6e612a5bf1f2aab478625e75927 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Fri, 17 Jul 2026 04:01:01 -0700 Subject: [PATCH 2/3] feat: Add ProjectedCoords plugged in via frame_transform --- docs/src/api.md | 1 + src/SkyCoords.jl | 2 ++ src/projected.jl | 71 ++++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 52 +++++++++++++++++++++++++++++++++++ 4 files changed, 126 insertions(+) create mode 100644 src/projected.jl diff --git a/docs/src/api.md b/docs/src/api.md index 4854920..e7dddec 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -22,6 +22,7 @@ FK4NoETerms FK5Coords EclipticCoords CartesianCoords +ProjectedCoords ``` ## Conversion diff --git a/src/SkyCoords.jl b/src/SkyCoords.jl index 6196edb..dc888a8 100644 --- a/src/SkyCoords.jl +++ b/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) diff --git a/src/projected.jl b/src/projected.jl new file mode 100644 index 0000000..7b514f5 --- /dev/null +++ b/src/projected.jl @@ -0,0 +1,71 @@ +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} +end + +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...) diff --git a/test/runtests.jl b/test/runtests.jl index 8e8fbf5..2c690e8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,6 +14,7 @@ using Test import Makie import SkyCoords: lat, lon +using SkyCoords: project, origin const rng = StableRNG(2000) rad2arcsec(r) = 3600 * rad2deg(r) @@ -21,6 +22,57 @@ 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) From 7272e3b2e1d8df1a6b3054fd6a8f8b9a28a84a3d Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Sat, 18 Jul 2026 16:56:36 -0700 Subject: [PATCH 3/3] docs: Add to API docs --- docs/src/api.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/api.md b/docs/src/api.md index e7dddec..26e6e38 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -66,6 +66,7 @@ inrange separation position_angle offset +SkyCoords.project cartesian spherical ```