WCS header support

Project.toml

@@ -10,9 +10,11 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [weakdeps]
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d"
 
 [extensions]
 FITSIOTablesExt = "Tables"
+WCSExt = "WCS"
 
 [compat]
 Aqua = "0.8"
@@ -23,6 +25,7 @@ Random = "<0.0.1, 1"
 Reexport = "0.2, 1.0"
 Tables = "1"
 Test = "<0.0.1, 1"
+WCS = "0.6.2"
 julia = "1.3"
 
 [extras]
@@ -31,6 +34,7 @@ OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d"
 
 [targets]
-test = ["Aqua", "OrderedCollections", "Random", "Tables", "Test"]
+test = ["Aqua", "OrderedCollections", "Random", "Tables", "Test", "WCS"]

docs/Project.toml

@@ -1,7 +1,9 @@
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
 FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
+WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d"
 
 [compat]
 Documenter = "1"

docs/make.jl

@@ -1,15 +1,21 @@
-using Documenter, FITSIO
+using Documenter, DocumenterInterLinks, FITSIO
+using DataFrames, WCS # Precompile package extensions
+
+links = InterLinks(
+    "WCS" => "https://juliaastro.org/WCS/stable/",
+)
 
 include("pages.jl")
 makedocs(;
-    modules = [FITSIO],
+    modules = [FITSIO, Base.get_extension(FITSIO, :WCSExt)],
     sitename = "FITSIO.jl",
     format = Documenter.HTML(
         prettyurls = get(ENV, "CI", nothing) == "true",
         canonical = "https://juliaastro.org/FITSIO/stable/",
     ),
     pages = pages,
     checkdocs = :exports,
+    plugins = [links],
 )
 
 deploydocs(;

docs/src/api.md

@@ -23,6 +23,8 @@ read_key
 write_key
 read_header
 FITSHeader
+FITSHeader(::AbstractVector{<:NamedTuple})
+FITSHeader(::WCS.WCSTransform)
 length(::FITSHeader)
 haskey(::FITSHeader, ::String)
 keys(::FITSHeader)

ext/WCSExt.jl

@@ -0,0 +1,71 @@
+module WCSExt
+
+using WCS: WCSTransform, to_header
+import FITSIO: FITSIO, FITSHeader
+
+"""
+    FITSHeader(wcs::WCS.WCSTransform)
+
+Construct a [`FITSHeader`](@ref) from a [`WCSTransform`](@extref WCS.WCSTransform) supplied by [WCS.jl](@extref).
+
+# Examples
+
+```jldoctest
+julia> using FITSIO, WCS
+
+julia> wcs = WCSTransform(2;
+           cdelt = [-0.066667, 0.066667],
+           ctype = ["RA---AIR", "DEC--AIR"],
+           crpix = [-234.75, 8.3393],
+           crval = [0., -90],
+           pv    = [(2, 1, 45.0)],
+       )
+WCSTransform(naxis=2, cdelt=[-0.066667, 0.066667], crval=[0.0, -90.0], crpix=[-234.75, 8.3393])
+
+julia> FITSHeader(wcs)
+WCSAXES = '2       '           / Number of coordinate axes
+CRPIX1  = '-234.7500'          / Pixel coordinate of reference point
+CRPIX2  = '8.3393  '           / Pixel coordinate of reference point
+CDELT1  = '-0.066667'          / [deg] Coordinate increment at reference point
+CDELT2  = '0.066667'           / [deg] Coordinate increment at reference point
+CUNIT1  = 'deg     '           / Units of coordinate increment and value
+CUNIT2  = 'deg     '           / Units of coordinate increment and value
+CTYPE1  = 'RA---AIR'           / Right ascension, Airys zenithal projection
+CTYPE2  = 'DEC--AIR'           / Declination, Airys zenithal projection
+CRVAL1  = '0.0     '           / [deg] Coordinate value at reference point
+CRVAL2  = '-90.0   '           / [deg] Coordinate value at reference point
+PV2_1   = '45.0    '           / AIR projection parameter
+LONPOLE = '180.0   '           / [deg] Native longitude of celestial pole
+LATPOLE = '-90.0   '           / [deg] Native latitude of celestial pole
+MJDREF  = '0.0     '           / [d] MJD of fiducial time
+RADESYS = 'ICRS    '           / Equatorial coordinate system
+COMMENT WCS header keyrecords produced by WCSLIB 7.7
+```
+"""
+function FITSIO.FITSHeader(wcs::WCSTransform)
+	# Split string into 80-character card images
+	card_images = Iterators.partition(to_header(wcs), 80)
+
+	# Remove any blank lines
+	is_empty = isempty ∘ strip
+	card_images = Iterators.filter(card_images) do card_image
+		!is_empty(card_image)
+	end
+
+	# Split each of those card images into their (key, value, comment) parts
+	card_image_parts = map(card_images) do card_image
+        card_image = replace(card_image, "'" => "") # Remove single quotes
+		map(strip, split(card_image, ['=', '/' ]))
+	end
+
+	# Deal with the special comment case
+    comment_card = (first ∘ pop!)(card_image_parts)
+	comment = (strip ∘ last)(split(comment_card, "COMMENT"))
+	push!(card_image_parts, ["COMMENT", "", comment])
+
+	# Store
+	k, v, c = eachcol(stack(card_image_parts; dims = 1))
+	return FITSHeader(string.(k), string.(v), string.(c))
+end
+
+end # module

src/FITSIO.jl

@@ -232,6 +232,19 @@ reading from a file. (This is similar to how an `Array` returned by
 if it was created by `read_header(::HDU)`.  You can, however, write a
 `FITSHeader` to a file using the `write(::FITS, ...)` methods that
 append a new HDU to a file.
+
+# Examples
+
+
+```jldoctest
+julia> using FITSIO
+
+julia> FITSHeader(["Key1", "Key2"], [1.0, "one"], ["Comment1", "Comment2"])
+Key1    =                  1.0 / Comment1
+Key2    = 'one     '           / Comment2
+```
+
+If [WCS.jl](@extref) is loaded, then a `FITSHeader` can also be constructed from a [`WCS.WCSTransform`](@extref).
 """
 mutable struct FITSHeader
     keys::Vector{String}

src/header.jl

@@ -6,6 +6,24 @@
 # Used here and in other files. Functions that operate on FITSFile
 # start with `fits_`.
 
+"""
+    FITSHeader(cards::AbstractVector{<:NamedTuple})
+
+Construct a [`FITSHeader`](@ref) from a vector of `NamedTuples` with the following fields: `key`, `value`, and `comment`.
+
+# Examples
+
+```jldoctest
+julia> using FITSIO
+
+julia> FITSHeader([
+           (key = "Key1", value = 1.0, comment = "Comment1"),
+           (key = "Key2", value = "one", comment = "Comment2"),
+       ])
+Key1    =                  1.0 / Comment1
+Key2    = 'one     '           / Comment2
+```
+"""
 FITSIO.FITSHeader(cards::AbstractVector{<:NamedTuple}) = FITSHeader(
     map(x -> x.key, cards),
     map(x -> x.value, cards),

test/runtests.jl

@@ -950,3 +950,5 @@ end
         end
     end
 end
+
+include("test_wcs.jl")

test/test_wcs.jl

@@ -0,0 +1,50 @@
+using WCS: WCSTransform
+
+@testset "WCS handling" begin
+    # Create sample fits data
+    img = [6 7; 8 9]
+    wcs = WCSTransform(2;
+        cdelt = [-0.066667, 0.066667],
+        ctype = ["RA---AIR", "DEC--AIR"],
+        crpix = [-234.75, 8.3393],
+        crval = [0., -90],
+        pv    = [(2, 1, 45.0)],
+    )
+    header_wcs = FITSHeader(wcs)
+
+    # Check output
+    header_default_str = """SIMPLE  =                    T / file does conform to FITS standard
+       BITPIX  =                   64 / number of bits per data pixel
+       NAXIS   =                    2 / number of data axes
+       NAXIS1  =                    2 / length of data axis 1
+       NAXIS2  =                    2 / length of data axis 2
+       EXTEND  =                    T / FITS dataset may contain extensions
+       COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronom
+       COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
+       """
+
+    header_wcs_str = """WCSAXES = '2       '           / Number of coordinate axes
+       CRPIX1  = '-234.7500'          / Pixel coordinate of reference point
+       CRPIX2  = '8.3393  '           / Pixel coordinate of reference point
+       CDELT1  = '-0.066667'          / [deg] Coordinate increment at reference point
+       CDELT2  = '0.066667'           / [deg] Coordinate increment at reference point
+       CUNIT1  = 'deg     '           / Units of coordinate increment and value
+       CUNIT2  = 'deg     '           / Units of coordinate increment and value
+       CTYPE1  = 'RA---AIR'           / Right ascension, Airys zenithal projection
+       CTYPE2  = 'DEC--AIR'           / Declination, Airys zenithal projection
+       CRVAL1  = '0.0     '           / [deg] Coordinate value at reference point
+       CRVAL2  = '-90.0   '           / [deg] Coordinate value at reference point
+       PV2_1   = '45.0    '           / AIR projection parameter
+       LONPOLE = '180.0   '           / [deg] Native longitude of celestial pole
+       LATPOLE = '-90.0   '           / [deg] Native latitude of celestial pole
+       MJDREF  = '0.0     '           / [d] MJD of fiducial time
+       RADESYS = 'ICRS    '           / Equatorial coordinate system
+       COMMENT WCS header keyrecords produced by WCSLIB 7.7"""
+
+    @test string(header_wcs) == header_wcs_str
+
+    tempnamefits() do fname
+        FITSIO.fitswrite(fname, img; header = header_wcs)
+        @test string(FITSIO.read_header(fname)) == header_default_str * header_wcs_str
+    end
+end