refactor: Use FITSFiles.jl + FITSWCS.jl

Project.toml

@@ -1,20 +1,30 @@
 name = "Reproject"
 uuid = "d1dcc2e6-806e-11e9-2897-3f99785db2ae"
 authors = ["Mosè Giordano", "Rohit Kumar"]
-version = "0.3.3"
+version = "0.4.0"
 
 [workspace]
 projects = ["test", "docs"]
 
 [deps]
-FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
+FITSFiles = "358a0a88-3548-4ad6-b652-8bdbf64af8e5"
+FITSWCS = "b9a10f0d-c86c-4e01-8e2e-1f1d7c5fb3f2"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 SkyCoords = "fc659fc5-75a3-5475-a2ea-3da92c065361"
-WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d"
+
+# Temporary until releases with the required changes are registered
+# (FITSWCS is not registered yet and FITSFiles needs an unreleased version).
+# Requires Julia 1.11+ to take effect.
+[sources]
+FITSFiles = {url = "https://github.com/JuliaAstro/FITSFiles.jl", rev = "hdu-size"}
+FITSWCS = {url = "https://github.com/JuliaAstro/FITSWCS.jl"}
+# The registered SkyCoords caps its Makie extension below the pre-release
+# Makie 0.25 used by the docs; this branch widens it.
+SkyCoords = {url = "https://github.com/JuliaAstro/SkyCoords.jl", rev = "makie-v0.25"}
 
 [compat]
-FITSIO = "0.15, 0.16, 0.17"
+FITSFiles = "0.4"
+FITSWCS = "0.1"
 Interpolations = "0.13, 0.14, 0.15, 0.16"
-SkyCoords = "0.4, 1"
-WCS = "0.5, 0.6"
+SkyCoords = "0.4, 1, 2"
 julia = "1.10"

docs/Project.toml

@@ -1,13 +1,34 @@
 [deps]
-ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
+AstroAngles = "5c4adb95-c1fc-4c53-b4ea-2a94080c53d2"
+AstroImages = "fe3fc30c-9b16-11e9-1c73-17dabf39f4ad"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+ComputePipeline = "95dc2771-c249-4cd0-9c9f-1f3b4330693c"
+DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
+DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
+FITSFiles = "358a0a88-3548-4ad6-b652-8bdbf64af8e5"
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
 ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Reproject = "d1dcc2e6-806e-11e9-2897-3f99785db2ae"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+# Temporary until releases with the required changes are registered. Pkg only
+# applies the [sources] of the active project, so AstroImages' own unregistered
+# requirements from its `makie` branch (AstroAngles 0.3, the DimensionalData
+# fork, and the pre-release Makie 0.25 stack) must be mirrored here, matching
+# that branch's own [sources].
+[sources]
+AstroAngles = {url = "https://github.com/JuliaAstro/AstroAngles.jl"}
+AstroImages = {url = "https://github.com/JuliaAstro/AstroImages.jl", rev = "makie"}
+CairoMakie = {url = "https://github.com/MakieOrg/Makie.jl", rev = "ff/breaking-0.25", subdir = "CairoMakie"}
+ComputePipeline = {url = "https://github.com/MakieOrg/Makie.jl", rev = "ff/breaking-0.25", subdir = "ComputePipeline"}
+DimensionalData = {url = "https://github.com/icweaver/DimensionalData.jl", rev = "compat-makie-0.25"}
+FITSFiles = {url = "https://github.com/JuliaAstro/FITSFiles.jl", rev = "hdu-size"}
+Makie = {url = "https://github.com/MakieOrg/Makie.jl", rev = "ff/breaking-0.25", subdir = "Makie"}
 
 [compat]
+AstroImages = "0.6"
 Documenter = "1"
+DocumenterInterLinks = "1"

docs/make.jl

@@ -1,7 +1,13 @@
 using Documenter
+using DocumenterInterLinks
 using Reproject
 using Documenter.Remotes: GitHub
 
+links = InterLinks(
+    "AstroImages" => "https://juliaastro.org/AstroImages/stable/",
+    "FITSFiles" => "https://juliaastro.org/FITSFiles/stable/",
+)
+
 makedocs(;
     sitename = "Reproject.jl",
     modules = [Reproject],
@@ -17,6 +23,7 @@ makedocs(;
     repo = GitHub("JuliaAstro/Reproject.jl"),
     # Doctests are run as part of the test suite instead (test/doctest.jl)
     doctest   = false,
+    plugins   = [links],
 )
 
 deploydocs(;

docs/src/index.md

@@ -1,67 +1,33 @@
 # Reproject
 
-Implementation in [Julia](https://julialang.org/) of the
-[`reproject`](https://github.com/astropy/reproject) package by Thomas
-Robitaille, part of the Astropy project.
+Implementation in [Julia](https://julialang.org/) of the [`reproject`](https://github.com/astropy/reproject) package by Thomas Robitaille, part of the Astropy project.
 
-This package can be used to reproject astronomical images from one world coordinate to another. By reproject we mean re-gridding of images using interpolation (i.e changing the pixel resolution, orientation, coordinate system).
+This package can be used to reproject astronomical images from one world coordinate to another using the exported [`reproject`](@ref) function. By reproject, we mean re-gridding of images using interpolation (i.e., changing the pixel resolution, orientation, coordinate system).
 
-Installation
--------
-
-Reproject.jl is avilable for Julia 1.0 and later versions and can be installed with [Julia built-in package manager](https://docs.julialang.org/en/v1/stdlib/Pkg/).
+## Installation
 
 ```julia-repl
-julia> import Pkg; Pkg.add("Reproject")
+pkg> add Reproject
 ```
 
-Usage
--------
-
-After installing the package, you can start using it with
+## Usage
 
 ```julia-repl
 julia> using Reproject
 
-julia> result = reproject(input_data, output_projection)
+julia> result, footprint = reproject(input_data, output_projection)
 ```
 
-This returns a tuple of reprojected image and footprint.
-
-
-## Reprojecting Images
+To reproject astronomical images, primary requirements are image data (2D Matrix), world coordinate frame of the input image, and the required output frame in which it needs to be reprojected.
 
-To reproject astronomical images, primary requirements are image data (2D Matrix), world cordinate frame of input image and required output frame in which it needs to be reprojected.
+The image data and input frame is given together as a [`FITSFiles.HDU`](@extref), a vector of HDUs as returned by [`FITSFiles.fits`](@extref FITSFiles.fits-Tuple{IO}), or path to the FITS file in `input_data`. A keyword argument `hdu_in` can be given while using a vector of HDUs or FITS file name to select the HDU in the FITS file, either by its 1-based index or by its `EXTNAME`.
 
-The image data and input frame is given together as an [ImageHDU](http://juliaastro.org/FITSIO) or [FITS](https://github.com/JuliaAstro/FITSIO.jl) file or name of the FITS file in `input_data`. A keyword argument `hdu_in` can be given while using FITS or FITS file name to specify specific HDU in FITS file.
-
-The `output_projection` is the output world coordinate frame and needs to be a a [WCSTransform](https://github.com/JuliaAstro/WCS.jl) or an [ImageHDU](http://juliaastro.org/FITSIO) or [FITS](https://github.com/JuliaAstro/FITSIO.jl) file or name of the FITS file. A keyword argument `hdu_out` can be given while using FITS or FITS file name to specify specific HDU in FITS file. WCS information is extracted from header when ImageHDU or FITS file is given as `output_projection`.
-
-Order of Interpolation can be specified by keyword `order` (i.e., 0, 1 (default), 2). The dimensions of output image can be given by keyword `shape_out`. This can be used to change resolution.
+The `output_projection` is the output world coordinate frame. It needs to be a [`FITSWCS.WCSTransform`](https://github.com/JuliaAstro/FITSWCS.jl), a [FITSFiles.HDU](@extref), a vector of HDUs, or the path to the FITS file. A keyword argument `hdu_out` can be given while using a vector of HDUs or FITS file name to select the HDU in the FITS file, either by its 1-based index or by its `EXTNAME`. WCS information is extracted from the header when an HDU or FITS file is given as `output_projection`.
 
+The order of the interpolation used can be specified by the `order` keyword (i.e., 0, 1 (default), 2). The dimensions of the output image can be given by the `shape_out` keyword. This can be used to change resolution.
 
 ## Example
 
-```julia-repl
-julia> using Reproject, FITSIO
-
-julia> input_data = FITS("gc_msx_e.fits");
-
-julia> output_projection = FITS("gc_2mass_k.fits");
-
-julia> result, footprint = reproject(input_data, output_projection, shape_out = (1000,1000), order = 2, hdu_in = 1, hdu_out = 1);
-
-julia> close(input_data); close(output_projection)
-```
-
-Remember to `close` any `FITS` handles you open once you are done with them. An open handle keeps the file locked on Windows. Alternatively, pass the file names directly and `reproject` will open and close the files for you:
-
-```julia-repl
-julia> result, footprint = reproject("gc_msx_e.fits", "gc_2mass_k.fits", shape_out = (1000,1000), order = 2)
-```
-
-The example below runs when this documentation is built. It reprojects an MSX E-band image of the Galactic center onto the frame of a 2MASS K-band image. `gc_msx_e` and `gc_2mass_k` are the paths to the two FITS files, which come from [astropy-data](https://www.astropy.org/astropy-data/) and are fetched by the documentation build as a lazy [Pkg artifact](https://pkgdocs.julialang.org/v1/artifacts/).
-
 ```@setup gc
 using Reproject
 using Pkg.Artifacts: ensure_artifact_installed
@@ -70,39 +36,48 @@ gc_msx_e = joinpath(data, "gc_msx_e.fits")
 gc_2mass_k = joinpath(data, "gc_2mass_k.fits")
 ```
 
+Using the following data from [astropy-data](https://www.astropy.org/astropy-data/):
+
+```@example gc
+@info "File paths:" gc_msx_e gc_2mass_k
+```
+
+we reprojects an MSX E-band image of the Galactic center (`gc_msx_e`) onto the frame of a 2MASS K-band image (`gc_2mass_k`):
+
 ```@example gc
 using Reproject
 
-result, footprint = reproject(gc_msx_e, gc_2mass_k; shape_out = (1000, 1000), order = 2)
-summary(result)
+result, footprint = reproject(gc_msx_e, gc_2mass_k, shape_out = (1000,1000), order = 2)
 ```
 
-To look at the images, we apply an asinh stretch and a colormap, and flip the matrix into row-major display orientation:
+To look at the images, we use [`imview`](@extref AstroImages.imview) from [AstroImages.jl](@extref AstroImages :doc:`index`), which applies percentile limits, an asinh stretch, and a colormap for us:
 
 ```@example gc
-using FITSIO, ImageShow, ImageCore, ColorSchemes, Statistics
-
-function render(img)
-    lo, hi = quantile(filter(isfinite, vec(img)), (0.02, 0.998))
-    stretched = asinh.(10 .* clamp.((img .- lo) ./ (hi - lo), 0, 1)) ./ asinh(10)
-    return map(v -> isfinite(v) ? get(ColorSchemes.hot, v) : RGB(0, 0, 0),
-               reverse(permutedims(stretched), dims = 1))
-end
-
-input_image = FITS(gc_msx_e) do f
-    read(f[1])
-end
-render(input_image)
+using AstroImages, FITSFiles
+
+input_image = load(gc_msx_e)
 ```
 
-**Output**, reprojected onto the 2MASS frame (black corners are pixels outside the input image's footprint, where `footprint` is `false`):
+**Output**, reprojected onto the 2MASS frame and shown next to its footprint (blank corners are pixels outside the input image's footprint, where `footprint` is `false`):
 
 ```@example gc
-render(result)
+using ImageCore
+
+mosaic(imview(result), imview(footprint); nrow = 1, npad = 10)
+```
+
+For a closer look with world coordinate axes and a colorbar, we can plot the reprojected image with `implotview` and a [Makie](https://docs.makie.org/) backend. `reproject` returns a plain matrix, so we first re-attach the output frame's WCS headers with [`copyheader`](@extref AstroImages.copyheader):
+
+```@example gc
+using CairoMakie
+
+target = load(gc_2mass_k)
+fig, _ = implotview(copyheader(target, result); clims = Percent(99.6), stretch = asinhstretch, cmap = :hot)
+fig
 ```
 
 !!! note
-    See [AstroImages.jl](https://juliaastro.org/AstroImages/stable/guide/reproject/) for more on using this package with astronomical images (e.g., FITS, WCS, and plotting support).
+    See [AstroImages.jl](@extref AstroImages :doc:`guide/reproject`) for more on using this package with astronomical images (e.g., FITS, WCS, and plotting support).
 
 ## License
 

src/Reproject.jl

@@ -1,6 +1,7 @@
 module Reproject
 
-using FITSIO: FITS, ImageHDU, read_header
+using FITSFiles: FITSFiles, HDU, Image, Primary, fits
+using FITSWCS: FITSWCS, SlicedWCSTransform, WCS, WCSTransform, pixel_to_world, world_to_pixel
 using Interpolations:
     BSpline,
     Constant,
@@ -12,7 +13,6 @@ using Interpolations:
     extrapolate,
     interpolate
 using SkyCoords: SkyCoords, FK5Coords, GalCoords, ICRSCoords
-using WCS: WCS, WCSTransform, pix_to_world, world_to_pix
 
 include("parsers.jl")
 include("utils.jl")

src/core.jl

@@ -5,70 +5,41 @@ Reprojects image data to a new projection using interpolation.
 
 # Arguments
 - `input_data`: Image data which is being reprojected.
-                It can be an ImageHDU, FITS object, name of a FITS file or a tuple of image matrix and WCSTransform.
+                It can be an HDU, a vector of HDUs as returned by `FITSFiles.fits`,
+                name of a FITS file or a tuple of image matrix and WCSTransform.
 - `output_projection`: Frame in which data is reprojected.
-                       Frame can be taken from WCSTransform object, ImageHDU, FITS or name of FITS file.
-- `shape_out`: Shape of image after reprojection.
+                       Frame can be taken from WCSTransform object, HDU, vector of HDUs or name of FITS file.
+- `shape_out`: Shape of image after reprojection; defaults to the data size of the selected HDU when the output projection is given as an HDU, vector of HDUs, or name of FITS file.
 - `order`: Order of interpolation.
            0: Nearest-neighbor
            1: Linear
            2: Quadratic
-- `hdu_in`: Used to specify HDU number when giving input as FITS or name of FITS file.
-- `hdu_out:` Used to specify HDU number when giving output projection as FITS or name of FITS file.
+- `hdu_in`: Used to select the HDU when giving input as a vector of HDUs or name of FITS file, either a 1-based integer index or an `EXTNAME` string or symbol.
+- `hdu_out:` Used to select the HDU when giving output projection as a vector of HDUs or name of FITS file, either a 1-based integer index or an `EXTNAME` string or symbol.
 """
-function reproject(input_data, output_projection; shape_out = nothing, order::Int = 1, hdu_in::Int = 1, hdu_out::Int = 1)
-    if input_data isa ImageHDU || input_data isa Tuple{AbstractArray, WCSTransform}
-        array_in, wcs_out = parse_input_data(input_data)
-    else
-        array_in, wcs_out = parse_input_data(input_data, hdu_in)
-    end
-
-    if output_projection isa FITS || output_projection isa String
-        wcs_in, shape_out = parse_output_projection(output_projection, hdu_out)
-    else
-        wcs_in, shape_out = parse_output_projection(output_projection, shape_out)
-    end
+function reproject(input_data, output_projection; shape_out = nothing, order::Int = 1, hdu_in::HDUSelector = 1, hdu_out::HDUSelector = 1)
+    array_in, wcs_source = parse_input_data(input_data, hdu_in)
+    wcs_target, shape_out = parse_output_projection(output_projection, shape_out; hdu_out)
 
-    type_in = wcs_to_celestial_frame(wcs_in)
-    type_out = wcs_to_celestial_frame(wcs_out)
-
-    if type_in == type_out && shape_out === nothing
-        return array_in
-    end
+    frame_source = celestial_frame(wcs_source)
+    frame_target = celestial_frame(wcs_target)
 
     img_out = fill(NaN, shape_out)
     itp = interpolator(array_in, order)
     shape_in = size(array_in)
 
     for i in 1:shape_out[1]
         for j in 1:shape_out[2]
-            pix_coord_in = [float(i), float(j)]
-            world_coord_in = pix_to_world(wcs_in, pix_coord_in)
+            pix_coord_target = [float(i), float(j)]
+            world_coord_target = pixel_to_world(wcs_target, pix_coord_target)
 
-            if type_in == "ICRS"
-                coord_in = ICRSCoords(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2]))
-            elseif type_in == "Gal"
-                coord_in = GalCoords(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2]))
-            elseif type_in == "FK5"
-                coord_in = FK5Coords{wcs_in.equinox}(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2]))
-            else
-                throw(ArgumentError("Unsupported output WCS coordinate type"))
-            end
-
-            if type_out == "ICRS"
-                coord_out = convert(ICRSCoords, coord_in)
-            elseif type_out == "Gal"
-                coord_out = convert(GalCoords, coord_in)
-            elseif type_out == "FK5"
-                coord_out = convert(FK5Coords{wcs_out.equinox}, coord_in)
-            else
-                throw(ArgumentError("Unsupported input WCS coordinate type"))
-            end
+            coord_target = frame_target(deg2rad(world_coord_target[1]), deg2rad(world_coord_target[2]))
+            coord_source = convert(frame_source, coord_target)
 
-            pix_coord_out = world_to_pix(wcs_out, [rad2deg(SkyCoords.lon(coord_out)), rad2deg(SkyCoords.lat(coord_out))])
+            pix_coord_source = world_to_pixel(wcs_source, [rad2deg(SkyCoords.lon(coord_source)), rad2deg(SkyCoords.lat(coord_source))])
 
-            if 0.5 <= pix_coord_out[1] <= shape_in[1] + 0.5 && 0.5 <= pix_coord_out[2] <= shape_in[2] + 0.5
-                img_out[i,j] = itp(pix_coord_out[1], pix_coord_out[2])
+            if 0.5 <= pix_coord_source[1] <= shape_in[1] + 0.5 && 0.5 <= pix_coord_source[2] <= shape_in[2] + 0.5
+                img_out[i,j] = itp(pix_coord_source[1], pix_coord_source[2])
             end
         end
     end