fix: Follow the FITS standard when RADESYS is absent, and support pv in the keyword constructor

Project.toml

@@ -16,7 +16,7 @@ FITSWCSFITSFilesExt = "FITSFiles"
 FITSWCSFITSIOExt = "FITSIO"
 
 [compat]
-FITSFiles = "0.2, 0.3"
+FITSFiles = "0.2, 0.3, 0.4"
 FITSIO = "0.17"
 StaticArrays = "1"
 julia = "1.10"

src/parsing.jl

@@ -1055,10 +1055,20 @@ function WCS(header::AbstractDict; fobj = nothing, alt::Char = ' ', minerr::Real
     end
 
     # ── Reference-frame metadata (pass-through, not used internally) ──────────
-    radesys = String(get(header, "RADESYS$(alt_str)",
-                     get(header, "RADESYS", "ICRS")))
     equinox = Float64(get(header, "EQUINOX$(alt_str)",
                       get(header, "EQUINOX", NaN)))
+    radesys_card = get(header, "RADESYS$(alt_str)", get(header, "RADESYS", nothing))
+    if radesys_card !== nothing
+        radesys = String(radesys_card)
+    elseif lon_axis > 0 && first(parse_ctype(ctype[lon_axis])) in ("RA", "ELON", "HLON")
+        # FITS WCS Paper II, Section 3.1: for equatorial and ecliptic systems an
+        # absent RADESYS defaults from EQUINOX — FK4 before 1984.0, FK5 from
+        # 1984.0 on, and ICRS when EQUINOX is also absent.
+        radesys = isnan(equinox) ? "ICRS" : (equinox < 1984.0 ? "FK4" : "FK5")
+    else
+        # RADESYS does not apply to non-equatorial/ecliptic coordinate systems.
+        radesys = ""
+    end
     wcsname = String(get(header, "WCSNAME$(alt_str)",
                      get(header, "WCSNAME", "")))
 
@@ -1087,7 +1097,9 @@ end
 Construct a transform from WCS.jl-style keyword vectors and matrices.
 
 Supported keywords are `crpix`, `crval`, `cdelt`, `ctype`, `cunit`, `pc`, `cd`,
-`crota`, `lonpole`, `latpole`, `radesys`, `equinox`, `wcsname`.
+`crota`, `pv`, `lonpole`, `latpole`, `radesys`, `equinox`, `wcsname`.
+`pv` takes a vector of `(axis, parameter, value)` tuples that map to the
+corresponding `PVi_m` keywords.
 More specialized wcslib fields remain unsupported in this pure-Julia constructor.
 """
 function WCS(naxis::Integer; preserve_units::Bool = false, kwds...)
@@ -1124,6 +1136,10 @@ function _constructor_keyword_to_header!(header::Dict{String,Any},
         _put_matrix_values!(header, "PC", value, naxis)
     elseif key === :cd
         _put_matrix_values!(header, "CD", value, naxis)
+    elseif key === :pv
+        for (i, m, v) in value
+            header["PV$(Int(i))_$(Int(m))"] = Float64(v)
+        end
     elseif key === :lonpole
         header["LONPOLE"] = value
     elseif key === :latpole

test/runtests.jl

@@ -3143,6 +3143,63 @@ end
             @test world2 ≈ world  atol=1e-10
         end
     end
+@testset "RADESYS defaulting" begin
+    # FITS WCS Paper II, Section 3.1: for equatorial and ecliptic systems an
+    # absent RADESYS defaults from EQUINOX — FK4 before 1984.0, FK5 from
+    # 1984.0 on, and ICRS when EQUINOX is also absent.
+    @test WCS(2; ctype = ["RA---AIR", "DEC--AIR"]).radesys == "ICRS"
+    @test WCS(2; ctype = ["RA---AIR", "DEC--AIR"], equinox = 1888.67).radesys == "FK4"
+    @test WCS(2; ctype = ["RA---AIR", "DEC--AIR"], equinox = 2000).radesys == "FK5"
+    @test WCS(2; ctype = ["ELON-AIT", "ELAT-AIT"], equinox = 2000).radesys == "FK5"
+
+    # RADESYS does not apply to non-equatorial/ecliptic coordinate systems.
+    @test WCS(2; ctype = ["GLON-CAR", "GLAT-CAR"]).radesys == ""
+    @test WCS(2; ctype = ["X", "Y"]).radesys == ""
+
+    # An explicit RADESYS always wins.
+    @test WCS(2; ctype = ["RA---AIR", "DEC--AIR"], radesys = "UNK").radesys == "UNK"
+    @test WCS(2; ctype = ["RA---AIR", "DEC--AIR"], radesys = "ICRS", equinox = 1950).radesys == "ICRS"
+
+    # The header-dict parsing path applies the same defaulting.
+    hdr = Dict{String,Any}(
+        "NAXIS" => 2,
+        "CTYPE1" => "RA---TAN", "CTYPE2" => "DEC--TAN",
+        "CRPIX1" => 1.0, "CRPIX2" => 1.0,
+        "CRVAL1" => 0.0, "CRVAL2" => 0.0,
+        "CDELT1" => 1.0, "CDELT2" => 1.0,
+        "EQUINOX" => 1950.0,
+    )
+    @test WCS(hdr).radesys == "FK4"
+    hdr["RADESYS"] = "FK5"
+    @test WCS(hdr).radesys == "FK5"
+end
+
+@testset "pv keyword constructor" begin
+    wcs = WCS(2;
+        cdelt = [-0.066667, 0.066667],
+        ctype = ["RA---AIR", "DEC--AIR"],
+        crpix = [-234.75, 8.3393],
+        crval = [0.0, -90.0],
+        pv = [(2, 1, 45.0)],
+    )
+    @test wcs isa WCSTransform
+    @test wcs.projection isa FITSWCS.AIR
+    @test wcs.projection.theta_b ≈ 45.0
+
+    # Equivalent to spelling out the corresponding PVi_m header keyword.
+    hdr = Dict{String,Any}(
+        "NAXIS" => 2,
+        "CTYPE1" => "RA---AIR", "CTYPE2" => "DEC--AIR",
+        "CRPIX1" => -234.75, "CRPIX2" => 8.3393,
+        "CRVAL1" => 0.0, "CRVAL2" => -90.0,
+        "CDELT1" => -0.066667, "CDELT2" => 0.066667,
+        "PV2_1" => 45.0,
+    )
+    fromheader = WCS(hdr)
+    p = [10.0, 20.0]
+    @test pixel_to_world(wcs, p) ≈ pixel_to_world(fromheader, p)
+end
+
 end  # @testset "FITSWCS"
 
 # Regression tests against wcslib reference values (Milestone 5).