From 06baabd5ac59bb6974246dcacfa03d31d0bb9796 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 23 Jul 2026 08:38:51 -0700 Subject: [PATCH 1/7] refactor: FITSFiles.jl + FITSWCS.jl up --- Project.toml | 20 ++++++--- docs/Project.toml | 27 +++++++++++-- docs/make.jl | 7 ++++ docs/src/index.md | 101 +++++++++++++++++----------------------------- src/Reproject.jl | 4 +- src/core.jl | 58 +++++++++++++------------- src/parsers.jl | 46 +++++++++------------ test/Project.toml | 10 ++++- test/core.jl | 34 +++++++--------- test/parsers.jl | 83 +++++++++++++++++-------------------- test/runtests.jl | 2 +- test/utils.jl | 12 +++--- 12 files changed, 200 insertions(+), 204 deletions(-) diff --git a/Project.toml b/Project.toml index 8a7facb..89126f7 100644 --- a/Project.toml +++ b/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", rev = "radesys-defaults"} +# 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" julia = "1.10" diff --git a/docs/Project.toml b/docs/Project.toml index e61ca3a..7bfb22d 100644 --- a/docs/Project.toml +++ b/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 = "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" diff --git a/docs/make.jl b/docs/make.jl index 3988e5d..4c2c940 100644 --- a/docs/make.jl +++ b/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(; diff --git a/docs/src/index.md b/docs/src/index.md index 8d8504f..d4b7fcd 100644 --- a/docs/src/index.md +++ b/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 specify the specific HDU in the FITS file. -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 specify the specific HDU in the FITS file. 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 diff --git a/src/Reproject.jl b/src/Reproject.jl index 2072c8b..1f707a1 100644 --- a/src/Reproject.jl +++ b/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, WCS, WCSTransform, pixel_to_world, world_to_pixel using Interpolations: BSpline, Constant, @@ -10,7 +11,6 @@ using Interpolations: Quadratic, interpolate using SkyCoords: SkyCoords, FK5Coords, GalCoords, ICRSCoords -using WCS: WCS, WCSTransform, pix_to_world, world_to_pix include("parsers.jl") include("utils.jl") diff --git a/src/core.jl b/src/core.jl index 5552bc5..785bab8 100644 --- a/src/core.jl +++ b/src/core.jl @@ -5,25 +5,26 @@ 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. + Frame can be taken from WCSTransform object, HDU, vector of HDUs or name of FITS file. - `shape_out`: Shape of image after reprojection. - `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 specify HDU number when giving input as a vector of HDUs or name of FITS file. +- `hdu_out:` Used to specify HDU number when giving output projection as a vector of HDUs or name of FITS file. """ 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} + if input_data isa HDU || 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 + if output_projection isa AbstractVector{<:HDU} || 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) @@ -36,6 +37,26 @@ function reproject(input_data, output_projection; shape_out = nothing, order::In return array_in end + frame_in = if type_in == "ICRS" + ICRSCoords + elseif type_in == "Gal" + GalCoords + elseif type_in == "FK5" + FK5Coords{wcs_in.equinox} + else + throw(ArgumentError("Unsupported output WCS coordinate type")) + end + + frame_out = if type_out == "ICRS" + ICRSCoords + elseif type_out == "Gal" + GalCoords + elseif type_out == "FK5" + FK5Coords{wcs_out.equinox} + else + throw(ArgumentError("Unsupported input WCS coordinate type")) + end + img_out = fill(NaN, shape_out) array_in = pad_edges(array_in) itp = interpolator(array_in, order) @@ -44,29 +65,12 @@ function reproject(input_data, output_projection; shape_out = nothing, order::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) - - 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 + world_coord_in = pixel_to_world(wcs_in, pix_coord_in) - 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_in = frame_in(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2])) + coord_out = convert(frame_out, coord_in) - pix_coord_out = world_to_pix(wcs_out, [rad2deg(SkyCoords.lon(coord_out)), rad2deg(SkyCoords.lat(coord_out))]) + pix_coord_out = world_to_pixel(wcs_out, [rad2deg(SkyCoords.lon(coord_out)), rad2deg(SkyCoords.lat(coord_out))]) if 0.5 <= pix_coord_out[1] <= shape_in[1] - 1.5 && 0.5 <= pix_coord_out[2] <= shape_in[2] - 1.5 img_out[i,j] = itp(pix_coord_out[1] + 1, pix_coord_out[2] + 1) diff --git a/src/parsers.jl b/src/parsers.jl index 3f1c65a..7eff9b0 100644 --- a/src/parsers.jl +++ b/src/parsers.jl @@ -1,18 +1,18 @@ """ - parse_input_data(input_data::ImageHDU) + parse_input_data(input_data::HDU) parse_input_data(input_data::Tuple{AbstractArray, WCSTransform}) parse_input_data(input_data::String, hdu_in) - parse_input_data(input_data::FITS, hdu_in) + parse_input_data(input_data::AbstractVector{<:HDU}, hdu_in) Parse input data and returns an Array and WCS object. # Arguments - `input_data`: image to reproject which can be name of a FITS file, - an ImageHDU or a FITS file. + an HDU or a vector of HDUs as returned by `FITSFiles.fits`. - `hdu_in`: used to set HDU to use when more than one HDU is present. """ -function parse_input_data(input_data::ImageHDU) - return read(input_data), WCS.from_header(read_header(input_data, String))[1] +function parse_input_data(input_data::HDU) + return input_data.data, WCS(input_data) end function parse_input_data(input_data::Tuple{AbstractArray, WCSTransform}) @@ -20,30 +20,25 @@ function parse_input_data(input_data::Tuple{AbstractArray, WCSTransform}) end function parse_input_data(input_data::String, hdu_in) - return FITS(input_data) do f - parse_input_data(f, hdu_in) - end + return parse_input_data(fits(input_data), hdu_in) end -function parse_input_data(input_data::FITS, hdu_in) +function parse_input_data(input_data::AbstractVector{<:HDU}, hdu_in) return parse_input_data(input_data[hdu_in]) end -# TODO: extend support for passing FITSHeader when FITSHeader to WCSTransform support is possible. - - """ parse_output_projection(output_projection::WCSTransform, shape_out) - parse_output_projection(output_projection::ImageHDU; shape_out) + parse_output_projection(output_projection::HDU, shape_out) parse_output_projection(output_projection::String, hdu_number) - parse_output_projection(output_projection::FITS, hdu_number) + parse_output_projection(output_projection::AbstractVector{<:HDU}, hdu_number) Parse output projection and returns a WCS object and shape of output. # Arguments - `output_projection`: WCS information about the image to be reprojected which can be - name of a FITS file, an ImageHDU or WCSTransform. + name of a FITS file, an HDU, a vector of HDUs or a WCSTransform. - `shape_out`: shape of the output image. - `hdu_number`: specifies HDU number when file name is given as input. """ @@ -55,8 +50,8 @@ function parse_output_projection(output_projection::WCSTransform, shape_out) return output_projection, shape_out end -function parse_output_projection(output_projection::ImageHDU, shape_out) - wcs_out = WCS.from_header(read_header(output_projection, String))[1] +function parse_output_projection(output_projection::HDU, shape_out) + wcs_out = WCS(output_projection) if shape_out === nothing shape_out = size(output_projection) end @@ -67,19 +62,14 @@ function parse_output_projection(output_projection::ImageHDU, shape_out) end function parse_output_projection(output_projection::String, hdu_number) - return FITS(output_projection) do f - parse_output_projection(f, hdu_number) - end + return parse_output_projection(fits(output_projection), hdu_number) end -function parse_output_projection(output_projection::FITS, hdu_number) - wcs_out = WCS.from_header(read_header(output_projection[hdu_number], String))[1] - - if output_projection[hdu_number] isa ImageHDU - shape_out = size(output_projection[hdu_number]) - else - throw(ArgumentError("Given FITS file doesn't have ImageHDU")) +function parse_output_projection(output_projection::AbstractVector{<:HDU}, hdu_number) + hdu = output_projection[hdu_number] + if !(hdu isa HDU{Primary} || hdu isa HDU{Image}) + throw(ArgumentError("Given FITS file doesn't have an image HDU")) end - return wcs_out, shape_out + return WCS(hdu), size(hdu) end diff --git a/test/Project.toml b/test/Project.toml index 221b475..bcbd178 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,12 +1,18 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb" +FITSFiles = "358a0a88-3548-4ad6-b652-8bdbf64af8e5" +FITSWCS = "b9a10f0d-c86c-4e01-8e2e-1f1d7c5fb3f2" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Reproject = "d1dcc2e6-806e-11e9-2897-3f99785db2ae" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d" [compat] PythonCall = "0.8, 0.9" + +# Temporary until releases with the required changes are registered; keep in +# sync with the root Project.toml. +[sources] +FITSFiles = {url = "https://github.com/JuliaAstro/FITSFiles.jl", rev = "hdu-size"} +FITSWCS = {url = "https://github.com/JuliaAstro/FITSWCS.jl", rev = "radesys-defaults"} diff --git a/test/core.jl b/test/core.jl index 72e1b8d..1180428 100644 --- a/test/core.jl +++ b/test/core.jl @@ -3,8 +3,8 @@ gc_msx_e = joinpath(data, "gc_msx_e.fits") gc_2mass_k = joinpath(data, "gc_2mass_k.fits") - imgin = FITS(gc_msx_e) # project this - imgout = FITS(gc_2mass_k) # into this coordinate + imgin = fits(gc_msx_e) # project this + imgout = fits(gc_2mass_k) # into this coordinate hdu1 = astropy.io.fits.open(gc_2mass_k)[0] hdu2 = astropy.io.fits.open(gc_msx_e)[0] @@ -31,26 +31,20 @@ @test isapprox(reproject(imgin[1], imgout[1], shape_out = (1000,1000))[1]', pyconvert(Matrix, rp.reproject_interp(hdu2, astropy.wcs.WCS(hdu1.header), shape_out = (1000,1000))[0]), nans = true, rtol = 1e-7) - wcs = WCSTransform(2; ctype = ["RA---AIR", "DEC--AIR"], radesys = "UNK") + wcs = WCS(2; ctype = ["RA---AIR", "DEC--AIR"], radesys = "UNK") @test_throws ArgumentError reproject(imgin, wcs, shape_out = (100,100)) - fname = tempname() * ".fits" - f = FITS(fname, "w") - inhdr = FITSHeader(["CTYPE1", "CTYPE2", "RADESYS", "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT", - "HISTORY"], - ["RA---TAN", "DEC--TAN", "UNK", 1.0, 1, true, "string value", nothing, nothing], - ["", - "", - "", - "floating point keyword", - "", - "boolean keyword", - "string value", - "this is a comment", - "this is a history"]) + inhdr = [Card("CTYPE1", "RA---TAN"), + Card("CTYPE2", "DEC--TAN"), + Card("RADESYS", "UNK"), + Card("FLTKEY", 1.0, "floating point keyword"), + Card("INTKEY", 1), + Card("BOOLKEY", true, "boolean keyword"), + Card("STRKEY", "string value", "string value"), + Card("COMMENT", "this is a comment"), + Card("HISTORY", "this is a history")] indata = reshape(Float32[1:100;], 5, 20) - write(f, indata; header=inhdr) - @test_throws ArgumentError reproject(f[1], imgin, shape_out = (100,100)) - close(f) + hdu = HDU(indata, inhdr) + @test_throws ArgumentError reproject(hdu, imgin, shape_out = (100,100)) end diff --git a/test/parsers.jl b/test/parsers.jl index 3cbb74b..62b6629 100644 --- a/test/parsers.jl +++ b/test/parsers.jl @@ -1,44 +1,40 @@ using Reproject: parse_input_data, parse_output_projection @testset "input parser" begin fname = tempname() * ".fits" - f = FITS(fname, "w") - inhdr = FITSHeader(["FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT", - "HISTORY"], - [1.0, 1, true, "string value", nothing, nothing], - ["floating point keyword", - "", - "boolean keyword", - "string value", - "this is a comment", - "this is a history"]) + inhdr = [Card("FLTKEY", 1.0, "floating point keyword"), + Card("INTKEY", 1), + Card("BOOLKEY", true, "boolean keyword"), + Card("STRKEY", "string value", "string value"), + Card("COMMENT", "this is a comment"), + Card("HISTORY", "this is a history")] indata = reshape(Float32[1:100;], 5, 20) - write(f, indata; header=inhdr) + hdu = HDU(indata, inhdr) + write(fname, [hdu]) - @testset "ImageHDU type" begin - result = parse_input_data(f[1]) + @testset "HDU type" begin + result = parse_input_data(hdu) @test result[1] isa Array @test result[2] isa WCSTransform end @testset "data matrix and WCSTransform tuple" begin - 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)]) + wcs = WCS(2; + cdelt = [-0.066667, 0.066667], + ctype = ["RA---AIR", "DEC--AIR"], + crpix = [-234.75, 8.3393], + crval = [0., -90], + pv = [(2, 1, 45.0)]) result = parse_input_data((indata, wcs)) @test result[1] isa Array @test result[2] isa WCSTransform end @testset "Single HDU FITS file" begin - result = parse_input_data(f, 1) + result = parse_input_data(fits(fname), 1) @test result[1] isa Array @test result[2] isa WCSTransform end - close(f) @testset "String filename input" begin result = parse_input_data(fname, 1) @@ -46,16 +42,13 @@ using Reproject: parse_input_data, parse_output_projection @test result[2] isa WCSTransform end - f = FITS(fname, "w") - write(f, indata; header=inhdr) - write(f, indata; header=inhdr) + write(fname, [hdu, HDU(Image, indata, inhdr)]) @testset "Multiple HDU FITS file" begin - result = parse_input_data(f, 2) + result = parse_input_data(fits(fname), 2) @test result[1] isa Array @test result[2] isa WCSTransform - close(f) result = parse_input_data(fname, 1) @test result[1] isa Array @test result[2] isa WCSTransform @@ -64,27 +57,23 @@ end @testset "output parser" begin fname = tempname() * ".fits" - f = FITS(fname, "w") - inhdr = FITSHeader(["FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT", - "HISTORY"], - [1.0, 1, true, "string value", nothing, nothing], - ["floating point keyword", - "", - "boolean keyword", - "string value", - "this is a comment", - "this is a history"]) + inhdr = [Card("FLTKEY", 1.0, "floating point keyword"), + Card("INTKEY", 1), + Card("BOOLKEY", true, "boolean keyword"), + Card("STRKEY", "string value", "string value"), + Card("COMMENT", "this is a comment"), + Card("HISTORY", "this is a history")] indata = reshape(Float32[1:100;], 5, 20) - write(f, indata; header=inhdr) + hdu = HDU(indata, inhdr) + write(fname, [hdu]) - @testset "ImageHDU type" begin - result = parse_output_projection(f[1], (12,12)) + @testset "HDU type" begin + result = parse_output_projection(hdu, (12,12)) @test result[1] isa WCSTransform @test result[2] isa Tuple - @test_throws DomainError parse_output_projection(f[1], ()) + @test_throws DomainError parse_output_projection(hdu, ()) end - close(f) @testset "String filename" begin result = parse_output_projection(fname, 1) @@ -92,12 +81,12 @@ end @test result[2] isa Tuple end - 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)]) + wcs = WCS(2; + cdelt = [-0.066667, 0.066667], + ctype = ["RA---AIR", "DEC--AIR"], + crpix = [-234.75, 8.3393], + crval = [0., -90], + pv = [(2, 1, 45.0)]) @testset "WCSTransform input" begin result = parse_output_projection(wcs, (12,12)) diff --git a/test/runtests.jl b/test/runtests.jl index b98024a..0b485cb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,7 @@ using Reproject using Test using PythonCall -using FITSIO, WCS +using FITSFiles, FITSWCS using Pkg.Artifacts: ensure_artifact_installed rp = pyimport("reproject") diff --git a/test/utils.jl b/test/utils.jl index 7ca58f1..0961a23 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,23 +1,23 @@ using Reproject: wcs_to_celestial_frame @testset "wcs to celestial frame" begin - wcs1 = WCSTransform(2; + wcs1 = WCS(2; ctype = ["RA---AIR", "DEC--AIR"], ) - wcs2 = WCSTransform(2; + wcs2 = WCS(2; ctype = ["RA---AIR", "DEC--AIR"], equinox = 1888.67 ) - wcs3 = WCSTransform(2; + wcs3 = WCS(2; ctype = ["RA---AIR", "DEC--AIR"], equinox = 2000 ) - wcs4 = WCSTransform(2; + wcs4 = WCS(2; ctype = ["GLON--", "GLAT--"], ) - wcs5 = WCSTransform(2; + wcs5 = WCS(2; ctype = ["TLON", "TLAT"], ) - wcs6 = WCSTransform(2; + wcs6 = WCS(2; ctype = ["RA---AIR", "DEC--AIR"], radesys = "UNK" ) From 49851c8c459163fd33d2f2f772420519d01e0c83 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 23 Jul 2026 13:39:02 -0700 Subject: [PATCH 2/7] refactor: Use padding built into Interpolations.jl instead of rolling our own --- src/Reproject.jl | 2 ++ src/core.jl | 17 ++++++++++------- src/utils.jl | 15 --------------- 3 files changed, 12 insertions(+), 22 deletions(-) diff --git a/src/Reproject.jl b/src/Reproject.jl index 1f707a1..44beaf9 100644 --- a/src/Reproject.jl +++ b/src/Reproject.jl @@ -5,10 +5,12 @@ using FITSWCS: FITSWCS, WCS, WCSTransform, pixel_to_world, world_to_pixel using Interpolations: BSpline, Constant, + Flat, InPlace, Linear, OnCell, Quadratic, + extrapolate, interpolate using SkyCoords: SkyCoords, FK5Coords, GalCoords, ICRSCoords diff --git a/src/core.jl b/src/core.jl index 785bab8..d859a43 100644 --- a/src/core.jl +++ b/src/core.jl @@ -58,7 +58,6 @@ function reproject(input_data, output_projection; shape_out = nothing, order::In end img_out = fill(NaN, shape_out) - array_in = pad_edges(array_in) itp = interpolator(array_in, order) shape_in = size(array_in) @@ -72,8 +71,8 @@ function reproject(input_data, output_projection; shape_out = nothing, order::In pix_coord_out = world_to_pixel(wcs_out, [rad2deg(SkyCoords.lon(coord_out)), rad2deg(SkyCoords.lat(coord_out))]) - if 0.5 <= pix_coord_out[1] <= shape_in[1] - 1.5 && 0.5 <= pix_coord_out[2] <= shape_in[2] - 1.5 - img_out[i,j] = itp(pix_coord_out[1] + 1, pix_coord_out[2] + 1) + 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]) end end end @@ -83,12 +82,16 @@ end """ - interpolator(array_in, order::Int) + interpolator(array_in, order::Int; padding = Flat()) Returns an interpolator with the given array and order of interpolation. + +`padding` is the Interpolations.jl boundary condition used for evaluations that +fall in the outer half-pixel border of the array, replacing the edge padding +that was previously applied to a copy of the input. """ -function interpolator(array_in::AbstractArray, order::Int) - if order == 0 +function interpolator(array_in::AbstractArray, order::Int; padding = Flat()) + if order == 0 itp = interpolate(array_in, BSpline(Constant())) elseif order == 1 itp = interpolate(array_in, BSpline(Linear())) @@ -96,5 +99,5 @@ function interpolator(array_in::AbstractArray, order::Int) itp = interpolate(array_in, BSpline(Quadratic(InPlace(OnCell())))) end - return itp + return extrapolate(itp, padding) end diff --git a/src/utils.jl b/src/utils.jl index 43a585d..a7cd8c4 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -20,18 +20,3 @@ function wcs_to_celestial_frame(wcs::WCSTransform) return radesys end - -""" - pad_edges(array_in::Matrix{T}) where {T} - -Pads a given array and creates a border with edge elements. -""" -function pad_edges(array_in::Matrix{T}) where {T} - image = Matrix{T}(undef, size(array_in)[1] + 2, size(array_in)[2] + 2) - image[2:end-1,2:end-1] = array_in - image[2:end-1,1] = array_in[:,1] - image[2:end-1,end] = array_in[:,end] - image[1,:] = image[2,:] - image[end,:] = image[end-1,:] - return image -end From cfc9fcd7e97561a9cba3cfda4b70f39e79669b58 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 23 Jul 2026 16:39:32 -0700 Subject: [PATCH 3/7] feat: More ways to select HDU --- docs/src/index.md | 4 ++-- src/core.jl | 6 +++--- src/parsers.jl | 28 +++++++++++++++++----------- test/parsers.jl | 16 +++++++++++++++- 4 files changed, 37 insertions(+), 17 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index d4b7fcd..6ec3ba5 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -20,9 +20,9 @@ julia> result, footprint = reproject(input_data, output_projection) 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. -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 specify the specific HDU in the FITS file. +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 `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 specify the specific HDU in the FITS file. WCS information is extracted from the header when an HDU or FITS file is given as `output_projection`. +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. diff --git a/src/core.jl b/src/core.jl index d859a43..5751ba9 100644 --- a/src/core.jl +++ b/src/core.jl @@ -14,10 +14,10 @@ Reprojects image data to a new projection using interpolation. 0: Nearest-neighbor 1: Linear 2: Quadratic -- `hdu_in`: Used to specify HDU number when giving input as a vector of HDUs or name of FITS file. -- `hdu_out:` Used to specify HDU number when giving output projection as a vector of HDUs 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) +function reproject(input_data, output_projection; shape_out = nothing, order::Int = 1, hdu_in::HDUSelector = 1, hdu_out::HDUSelector = 1) if input_data isa HDU || input_data isa Tuple{AbstractArray, WCSTransform} array_in, wcs_out = parse_input_data(input_data) else diff --git a/src/parsers.jl b/src/parsers.jl index 7eff9b0..8956a4d 100644 --- a/src/parsers.jl +++ b/src/parsers.jl @@ -1,8 +1,15 @@ +""" + HDUSelector + +An HDU within a FITS file, selected either by its 1-based position (an integer) or by its `EXTNAME` (a string or symbol). +""" +const HDUSelector = Union{Integer, AbstractString, Symbol} + """ parse_input_data(input_data::HDU) parse_input_data(input_data::Tuple{AbstractArray, WCSTransform}) - parse_input_data(input_data::String, hdu_in) - parse_input_data(input_data::AbstractVector{<:HDU}, hdu_in) + parse_input_data(input_data::String, hdu_in::HDUSelector) + parse_input_data(input_data::AbstractVector{<:HDU}, hdu_in::HDUSelector) Parse input data and returns an Array and WCS object. @@ -19,11 +26,11 @@ function parse_input_data(input_data::Tuple{AbstractArray, WCSTransform}) return input_data[1], input_data[2] end -function parse_input_data(input_data::String, hdu_in) +function parse_input_data(input_data::String, hdu_in::HDUSelector) return parse_input_data(fits(input_data), hdu_in) end -function parse_input_data(input_data::AbstractVector{<:HDU}, hdu_in) +function parse_input_data(input_data::AbstractVector{<:HDU}, hdu_in::HDUSelector) return parse_input_data(input_data[hdu_in]) end @@ -31,16 +38,15 @@ end """ parse_output_projection(output_projection::WCSTransform, shape_out) parse_output_projection(output_projection::HDU, shape_out) - parse_output_projection(output_projection::String, hdu_number) - parse_output_projection(output_projection::AbstractVector{<:HDU}, hdu_number) + parse_output_projection(output_projection::String, hdu_number::HDUSelector) + parse_output_projection(output_projection::AbstractVector{<:HDU}, hdu_number::HDUSelector) Parse output projection and returns a WCS object and shape of output. # Arguments -- `output_projection`: WCS information about the image to be reprojected which can be - name of a FITS file, an HDU, a vector of HDUs or a WCSTransform. +- `output_projection`: WCS information about the image to be reprojected which can be name of a FITS file, an HDU, a vector of HDUs or a WCSTransform. - `shape_out`: shape of the output image. -- `hdu_number`: specifies HDU number when file name is given as input. +- `hdu_number`: specifies the HDU when a file name or vector of HDUs is given as input. """ function parse_output_projection(output_projection::WCSTransform, shape_out) if length(shape_out) == 0 @@ -61,11 +67,11 @@ function parse_output_projection(output_projection::HDU, shape_out) return wcs_out, shape_out end -function parse_output_projection(output_projection::String, hdu_number) +function parse_output_projection(output_projection::String, hdu_number::HDUSelector) return parse_output_projection(fits(output_projection), hdu_number) end -function parse_output_projection(output_projection::AbstractVector{<:HDU}, hdu_number) +function parse_output_projection(output_projection::AbstractVector{<:HDU}, hdu_number::HDUSelector) hdu = output_projection[hdu_number] if !(hdu isa HDU{Primary} || hdu isa HDU{Image}) throw(ArgumentError("Given FITS file doesn't have an image HDU")) diff --git a/test/parsers.jl b/test/parsers.jl index 62b6629..e72b5ab 100644 --- a/test/parsers.jl +++ b/test/parsers.jl @@ -42,7 +42,7 @@ using Reproject: parse_input_data, parse_output_projection @test result[2] isa WCSTransform end - write(fname, [hdu, HDU(Image, indata, inhdr)]) + write(fname, [hdu, HDU(Image, indata, [Card("EXTNAME", "SCI"); inhdr])]) @testset "Multiple HDU FITS file" begin result = parse_input_data(fits(fname), 2) @@ -53,6 +53,20 @@ using Reproject: parse_input_data, parse_output_projection @test result[1] isa Array @test result[2] isa WCSTransform end + + @testset "HDU selection by EXTNAME" begin + result = parse_input_data(fname, "SCI") + @test result[1] isa Array + @test result[2] isa WCSTransform + + result = parse_input_data(fits(fname), :SCI) + @test result[1] isa Array + @test result[2] isa WCSTransform + + result = parse_output_projection(fname, "SCI") + @test result[1] isa WCSTransform + @test result[2] isa Tuple + end end @testset "output parser" begin From e9c43e9488cb422ed0ea0c5d92fb90ba2a8fd78c Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 23 Jul 2026 17:24:30 -0700 Subject: [PATCH 4/7] refactor: Simplify parsing --- src/core.jl | 57 +++++++------------------------ src/parsers.jl | 85 ++++++++++++++++++++--------------------------- src/utils.jl | 30 ++++++++--------- test/Project.toml | 1 + test/core.jl | 4 +++ test/parsers.jl | 7 ++-- test/utils.jl | 50 +++++++++++++++------------- 7 files changed, 99 insertions(+), 135 deletions(-) diff --git a/src/core.jl b/src/core.jl index 5751ba9..b7cce82 100644 --- a/src/core.jl +++ b/src/core.jl @@ -9,7 +9,7 @@ Reprojects image data to a new projection using interpolation. 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, HDU, vector of HDUs or name of FITS file. -- `shape_out`: Shape of image after reprojection. +- `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 @@ -18,44 +18,11 @@ Reprojects image data to a new projection using interpolation. - `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::HDUSelector = 1, hdu_out::HDUSelector = 1) - if input_data isa HDU || 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 + array_in, wcs_source = parse_input_data(input_data, hdu_in) + wcs_target, shape_out = parse_output_projection(output_projection, shape_out; hdu_out) - if output_projection isa AbstractVector{<:HDU} || 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 - - 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_in = if type_in == "ICRS" - ICRSCoords - elseif type_in == "Gal" - GalCoords - elseif type_in == "FK5" - FK5Coords{wcs_in.equinox} - else - throw(ArgumentError("Unsupported output WCS coordinate type")) - end - - frame_out = if type_out == "ICRS" - ICRSCoords - elseif type_out == "Gal" - GalCoords - elseif type_out == "FK5" - FK5Coords{wcs_out.equinox} - else - throw(ArgumentError("Unsupported input WCS coordinate type")) - end + frame_source = celestial_frame(wcs_source) + frame_target = celestial_frame(wcs_target) img_out = fill(NaN, shape_out) itp = interpolator(array_in, order) @@ -63,16 +30,16 @@ function reproject(input_data, output_projection; shape_out = nothing, order::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 = pixel_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) - coord_in = frame_in(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2])) - coord_out = convert(frame_out, coord_in) + 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_pixel(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 diff --git a/src/parsers.jl b/src/parsers.jl index 8956a4d..ff3ac6d 100644 --- a/src/parsers.jl +++ b/src/parsers.jl @@ -6,76 +6,63 @@ An HDU within a FITS file, selected either by its 1-based position (an integer) const HDUSelector = Union{Integer, AbstractString, Symbol} """ - parse_input_data(input_data::HDU) - parse_input_data(input_data::Tuple{AbstractArray, WCSTransform}) - parse_input_data(input_data::String, hdu_in::HDUSelector) - parse_input_data(input_data::AbstractVector{<:HDU}, hdu_in::HDUSelector) + select_hdu(x, hdu::HDUSelector) -Parse input data and returns an Array and WCS object. +Select an HDU from a FITS file name or a vector of HDUs, deferring to `FITSFiles.fits` and its indexing. Anything else (an HDU, a WCSTransform, a data/transform tuple) already refers to a single item and passes through unchanged. +""" +select_hdu(x, ::HDUSelector) = x +select_hdu(path::AbstractString, hdu::HDUSelector) = fits(path)[hdu] +select_hdu(hdus::AbstractVector{<:HDU}, hdu::HDUSelector) = hdus[hdu] -# Arguments -- `input_data`: image to reproject which can be name of a FITS file, - an HDU or a vector of HDUs as returned by `FITSFiles.fits`. -- `hdu_in`: used to set HDU to use when more than one HDU is present. """ -function parse_input_data(input_data::HDU) - return input_data.data, WCS(input_data) -end + parse_input_data(input_data, hdu_in::HDUSelector = 1) -function parse_input_data(input_data::Tuple{AbstractArray, WCSTransform}) - return input_data[1], input_data[2] -end +Parse input data and return an Array and WCS object. -function parse_input_data(input_data::String, hdu_in::HDUSelector) - return parse_input_data(fits(input_data), hdu_in) -end - -function parse_input_data(input_data::AbstractVector{<:HDU}, hdu_in::HDUSelector) - return parse_input_data(input_data[hdu_in]) +# Arguments +- `input_data`: image to reproject which can be the name of a FITS file, an HDU, a vector of HDUs as returned by `FITSFiles.fits`, or a tuple of an image matrix and a WCSTransform. +- `hdu_in`: selects the HDU when `input_data` is a file name or a vector of HDUs. +""" +function parse_input_data(input_data, hdu_in::HDUSelector = 1) + return data_and_wcs(select_hdu(input_data, hdu_in)) end +data_and_wcs(hdu::HDU) = hdu.data, WCS(hdu) +data_and_wcs(input_data::Tuple{AbstractArray, WCSTransform}) = input_data """ - parse_output_projection(output_projection::WCSTransform, shape_out) - parse_output_projection(output_projection::HDU, shape_out) - parse_output_projection(output_projection::String, hdu_number::HDUSelector) - parse_output_projection(output_projection::AbstractVector{<:HDU}, hdu_number::HDUSelector) + parse_output_projection(output_projection, shape_out = nothing; hdu_out::HDUSelector = 1) -Parse output projection and returns a WCS object and shape of output. +Parse the output projection and return a WCS object and the shape of the output. # Arguments -- `output_projection`: WCS information about the image to be reprojected which can be name of a FITS file, an HDU, a vector of HDUs or a WCSTransform. -- `shape_out`: shape of the output image. -- `hdu_number`: specifies the HDU when a file name or vector of HDUs is given as input. +- `output_projection`: WCS information about the image to be reprojected which can be the name of a FITS file, an HDU, a vector of HDUs, or a WCSTransform. +- `shape_out`: shape of the output image. Defaults to the data size of the selected HDU when one is given. +- `hdu_out`: selects the HDU when `output_projection` is a file name or a vector of HDUs. """ -function parse_output_projection(output_projection::WCSTransform, shape_out) - if length(shape_out) == 0 - throw(DomainError(shape_out, "The shape of the output image should not be an empty tuple")) - end - - return output_projection, shape_out +function parse_output_projection(output_projection, shape_out = nothing; hdu_out::HDUSelector = 1) + return wcs_and_shape(select_hdu(output_projection, hdu_out), shape_out) end -function parse_output_projection(output_projection::HDU, shape_out) - wcs_out = WCS(output_projection) +function wcs_and_shape(wcs::WCSTransform, shape_out) if shape_out === nothing - shape_out = size(output_projection) + throw(ArgumentError("`shape_out` must be given when the output projection is a WCSTransform")) end - if length(shape_out) == 0 - throw(DomainError(shape_out, "The shape of the output image should not be an empty tuple")) - end - return wcs_out, shape_out + validate_shape(shape_out) + return wcs, shape_out end -function parse_output_projection(output_projection::String, hdu_number::HDUSelector) - return parse_output_projection(fits(output_projection), hdu_number) -end - -function parse_output_projection(output_projection::AbstractVector{<:HDU}, hdu_number::HDUSelector) - hdu = output_projection[hdu_number] +function wcs_and_shape(hdu::HDU, shape_out) if !(hdu isa HDU{Primary} || hdu isa HDU{Image}) throw(ArgumentError("Given FITS file doesn't have an image HDU")) end + shape = something(shape_out, size(hdu)) + validate_shape(shape) + return WCS(hdu), shape +end - return WCS(hdu), size(hdu) +function validate_shape(shape_out) + if length(shape_out) == 0 + throw(DomainError(shape_out, "The shape of the output image should not be an empty tuple.")) + end end diff --git a/src/utils.jl b/src/utils.jl index a7cd8c4..8220057 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,22 +1,20 @@ +# TODO: Consider upstreaming this as an extension to FITSWCS.jl """ - wcs_to_celestial_frame(wcs::WCSTransform) + celestial_frame(wcs::WCSTransform) -Returns the reference frame of a WCSTransform. -The reference frame supported in Julia are FK5, ICRS and Galactic. +The `SkyCoords` coordinate type corresponding to the celestial reference frame of a WCS transform: `ICRSCoords`, `FK5Coords{equinox}`, or `GalCoords`. + +Throws an `ArgumentError` for frames without a `SkyCoords` counterpart (e.g., FK4 or ITRS). """ -function wcs_to_celestial_frame(wcs::WCSTransform) +function celestial_frame(wcs::WCSTransform) radesys = wcs.radesys - - xcoord = wcs.ctype[1][1:4] - ycoord = wcs.ctype[2][1:4] - - if radesys == "" - if xcoord == "GLON" && ycoord == "GLAT" - radesys = "Gal" - elseif xcoord == "TLON" && ycoord == "TLAT" - radesys = "ITRS" - end + if radesys == "ICRS" + return ICRSCoords + elseif radesys == "FK5" + return FK5Coords{wcs.equinox} + elseif radesys == "" && startswith(wcs.ctype[1], "GLON") && startswith(wcs.ctype[2], "GLAT") + return GalCoords end - - return radesys + frame = isempty(radesys) ? join(wcs.ctype, ", ") : radesys + throw(ArgumentError("WCS celestial frame ($frame) is not supported by SkyCoords")) end diff --git a/test/Project.toml b/test/Project.toml index bcbd178..e35e1a3 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -6,6 +6,7 @@ FITSWCS = "b9a10f0d-c86c-4e01-8e2e-1f1d7c5fb3f2" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Reproject = "d1dcc2e6-806e-11e9-2897-3f99785db2ae" +SkyCoords = "fc659fc5-75a3-5475-a2ea-3da92c065361" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] diff --git a/test/core.jl b/test/core.jl index 1180428..61bc9c5 100644 --- a/test/core.jl +++ b/test/core.jl @@ -31,6 +31,10 @@ @test isapprox(reproject(imgin[1], imgout[1], shape_out = (1000,1000))[1]', pyconvert(Matrix, rp.reproject_interp(hdu2, astropy.wcs.WCS(hdu1.header), shape_out = (1000,1000))[0]), nans = true, rtol = 1e-7) + # shape_out used to be silently ignored when the output projection was + # given as a file name or vector of HDUs + @test size(reproject(gc_msx_e, gc_2mass_k, shape_out = (100,100), order = 0)[1]) == (100, 100) + wcs = WCS(2; ctype = ["RA---AIR", "DEC--AIR"], radesys = "UNK") @test_throws ArgumentError reproject(imgin, wcs, shape_out = (100,100)) diff --git a/test/parsers.jl b/test/parsers.jl index e72b5ab..9c0d726 100644 --- a/test/parsers.jl +++ b/test/parsers.jl @@ -63,7 +63,7 @@ using Reproject: parse_input_data, parse_output_projection @test result[1] isa Array @test result[2] isa WCSTransform - result = parse_output_projection(fname, "SCI") + result = parse_output_projection(fname; hdu_out = "SCI") @test result[1] isa WCSTransform @test result[2] isa Tuple end @@ -90,9 +90,12 @@ end end @testset "String filename" begin - result = parse_output_projection(fname, 1) + result = parse_output_projection(fname) @test result[1] isa WCSTransform @test result[2] isa Tuple + + result = parse_output_projection(fname, (12, 12)) + @test result[2] == (12, 12) end wcs = WCS(2; diff --git a/test/utils.jl b/test/utils.jl index 0961a23..7ae04ef 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,31 +1,35 @@ -using Reproject: wcs_to_celestial_frame -@testset "wcs to celestial frame" begin +using Reproject: celestial_frame +using SkyCoords: ICRSCoords, GalCoords, FK5Coords + +@testset "celestial frame from WCS" begin wcs1 = WCS(2; - ctype = ["RA---AIR", "DEC--AIR"], - ) + ctype = ["RA---AIR", "DEC--AIR"], + ) wcs2 = WCS(2; - ctype = ["RA---AIR", "DEC--AIR"], - equinox = 1888.67 - ) + ctype = ["RA---AIR", "DEC--AIR"], + equinox = 1888.67 + ) wcs3 = WCS(2; - ctype = ["RA---AIR", "DEC--AIR"], - equinox = 2000 - ) + ctype = ["RA---AIR", "DEC--AIR"], + equinox = 2000 + ) wcs4 = WCS(2; - ctype = ["GLON--", "GLAT--"], - ) + ctype = ["GLON--", "GLAT--"], + ) wcs5 = WCS(2; - ctype = ["TLON", "TLAT"], - ) + ctype = ["TLON", "TLAT"], + ) wcs6 = WCS(2; - ctype = ["RA---AIR", "DEC--AIR"], - radesys = "UNK" - ) + ctype = ["RA---AIR", "DEC--AIR"], + radesys = "UNK" + ) - @test wcs_to_celestial_frame(wcs1) == "ICRS" - @test wcs_to_celestial_frame(wcs2) == "FK4" - @test wcs_to_celestial_frame(wcs3) == "FK5" - @test wcs_to_celestial_frame(wcs4) == "Gal" - @test wcs_to_celestial_frame(wcs5) == "ITRS" - @test wcs_to_celestial_frame(wcs6) == "UNK" + @test celestial_frame(wcs1) == ICRSCoords + # FK4 (pre-1984 equinox) has no SkyCoords counterpart + @test_throws ArgumentError celestial_frame(wcs2) + @test celestial_frame(wcs3) == FK5Coords{2000.0} + @test celestial_frame(wcs4) == GalCoords + # ITRS (TLON/TLAT) has no SkyCoords counterpart + @test_throws ArgumentError celestial_frame(wcs5) + @test_throws ArgumentError celestial_frame(wcs6) end From 39f7eb16ee4825f36bcdc3e2845747f9284bc35b Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Tue, 13 Oct 2026 01:17:39 -0700 Subject: [PATCH 5/7] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 89126f7..a6e8720 100644 --- a/Project.toml +++ b/Project.toml @@ -26,5 +26,5 @@ SkyCoords = {url = "https://github.com/JuliaAstro/SkyCoords.jl", rev = "makie-v0 FITSFiles = "0.4" FITSWCS = "0.1" Interpolations = "0.13, 0.14, 0.15, 0.16" -SkyCoords = "0.4, 1" +SkyCoords = "0.4, 1, 2" julia = "1.10" From cd9ef60d29c0ac1426f2bcc87deedbb6d60ef6aa Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Tue, 13 Oct 2026 02:47:35 -0700 Subject: [PATCH 6/7] feat: Add `SlicedWCSTransform` support --- Project.toml | 2 +- src/Reproject.jl | 2 +- src/utils.jl | 12 ++++++++++++ test/utils.jl | 9 +++++++++ 4 files changed, 23 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index a6e8720..dc93c3e 100644 --- a/Project.toml +++ b/Project.toml @@ -17,7 +17,7 @@ SkyCoords = "fc659fc5-75a3-5475-a2ea-3da92c065361" # 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", rev = "radesys-defaults"} +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"} diff --git a/src/Reproject.jl b/src/Reproject.jl index 44beaf9..5274ff7 100644 --- a/src/Reproject.jl +++ b/src/Reproject.jl @@ -1,7 +1,7 @@ module Reproject using FITSFiles: FITSFiles, HDU, Image, Primary, fits -using FITSWCS: FITSWCS, WCS, WCSTransform, pixel_to_world, world_to_pixel +using FITSWCS: FITSWCS, SlicedWCSTransform, WCS, WCSTransform, pixel_to_world, world_to_pixel using Interpolations: BSpline, Constant, diff --git a/src/utils.jl b/src/utils.jl index 8220057..146ea0e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,9 +1,12 @@ # TODO: Consider upstreaming this as an extension to FITSWCS.jl """ celestial_frame(wcs::WCSTransform) + celestial_frame(wcs::SlicedWCSTransform) The `SkyCoords` coordinate type corresponding to the celestial reference frame of a WCS transform: `ICRSCoords`, `FK5Coords{equinox}`, or `GalCoords`. +A `SlicedWCSTransform` (e.g., of a cutout, or of a single channel of a data cube) takes the frame of its parent, and must keep exactly the two celestial axes of its parent. + Throws an `ArgumentError` for frames without a `SkyCoords` counterpart (e.g., FK4 or ITRS). """ function celestial_frame(wcs::WCSTransform) @@ -18,3 +21,12 @@ function celestial_frame(wcs::WCSTransform) frame = isempty(radesys) ? join(wcs.ctype, ", ") : radesys throw(ArgumentError("WCS celestial frame ($frame) is not supported by SkyCoords")) end + +function celestial_frame(wcs::SlicedWCSTransform) + p = wcs.parent + if Set(wcs.world_keep) != Set((p.lon_axis, p.lat_axis)) + kept = join(p.ctype[wcs.world_keep], ", ") + throw(ArgumentError("Sliced WCS must keep exactly the two celestial axes of its parent, but keeps: $kept")) + end + return celestial_frame(p) +end diff --git a/test/utils.jl b/test/utils.jl index 7ae04ef..cab1fc4 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -33,3 +33,12 @@ using SkyCoords: ICRSCoords, GalCoords, FK5Coords @test_throws ArgumentError celestial_frame(wcs5) @test_throws ArgumentError celestial_frame(wcs6) end + +@testset "celestial frame from sliced WCS" begin + cube = WCS(3; ctype = ["GLON-TAN", "GLAT-TAN", "VRAD"]) + + # Cutout of a single channel: keeps both celestial axes + @test celestial_frame(slice_wcs(cube, 10:20, 5:15, 3)) == GalCoords + # Position-velocity slice: keeps longitude and velocity + @test_throws ArgumentError celestial_frame(slice_wcs(cube, :, 7, :)) +end From 276625f0b2cc0c5539be97aba4b16483722b1350 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Tue, 13 Oct 2026 03:39:13 -0700 Subject: [PATCH 7/7] typo --- docs/Project.toml | 2 +- test/Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 7bfb22d..5e35f60 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -24,7 +24,7 @@ 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 = "makie-0.25"} +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"} diff --git a/test/Project.toml b/test/Project.toml index e35e1a3..b308009 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -16,4 +16,4 @@ PythonCall = "0.8, 0.9" # sync with the root Project.toml. [sources] FITSFiles = {url = "https://github.com/JuliaAstro/FITSFiles.jl", rev = "hdu-size"} -FITSWCS = {url = "https://github.com/JuliaAstro/FITSWCS.jl", rev = "radesys-defaults"} +FITSWCS = {url = "https://github.com/JuliaAstro/FITSWCS.jl"}