diff --git a/Project.toml b/Project.toml index 8a7facb..dc93c3e 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"} +# 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" diff --git a/docs/Project.toml b/docs/Project.toml index e61ca3a..5e35f60 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 = "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" 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..6ec3ba5 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 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 diff --git a/src/Reproject.jl b/src/Reproject.jl index 7a31db4..5274ff7 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, 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") diff --git a/src/core.jl b/src/core.jl index 5c9de50..b7cce82 100644 --- a/src/core.jl +++ b/src/core.jl @@ -5,36 +5,24 @@ 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) @@ -42,33 +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 = 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 diff --git a/src/parsers.jl b/src/parsers.jl index 3f1c65a..ff3ac6d 100644 --- a/src/parsers.jl +++ b/src/parsers.jl @@ -1,85 +1,68 @@ """ - parse_input_data(input_data::ImageHDU) - parse_input_data(input_data::Tuple{AbstractArray, WCSTransform}) - parse_input_data(input_data::String, hdu_in) - parse_input_data(input_data::FITS, hdu_in) + HDUSelector -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. -- `hdu_in`: used to set HDU to use when more than one HDU is present. +An HDU within a FITS file, selected either by its 1-based position (an integer) or by its `EXTNAME` (a string or symbol). """ -function parse_input_data(input_data::ImageHDU) - return read(input_data), WCS.from_header(read_header(input_data, String))[1] -end +const HDUSelector = Union{Integer, AbstractString, Symbol} -function parse_input_data(input_data::Tuple{AbstractArray, WCSTransform}) - return input_data[1], input_data[2] -end +""" + select_hdu(x, hdu::HDUSelector) -function parse_input_data(input_data::String, hdu_in) - return FITS(input_data) do f - parse_input_data(f, hdu_in) - end -end +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] -function parse_input_data(input_data::FITS, hdu_in) - return parse_input_data(input_data[hdu_in]) -end +""" + parse_input_data(input_data, hdu_in::HDUSelector = 1) +Parse input data and return an Array and WCS object. -# TODO: extend support for passing FITSHeader when FITSHeader to WCSTransform support is possible. +# 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::ImageHDU; shape_out) - parse_output_projection(output_projection::String, hdu_number) - parse_output_projection(output_projection::FITS, hdu_number) + 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 ImageHDU or WCSTransform. -- `shape_out`: shape of the output image. -- `hdu_number`: specifies HDU number when file name 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::ImageHDU, shape_out) - wcs_out = WCS.from_header(read_header(output_projection, String))[1] +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) - return FITS(output_projection) do f - parse_output_projection(f, 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 -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 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 - - return wcs_out, shape_out end diff --git a/src/utils.jl b/src/utils.jl index a7cd8c4..146ea0e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,22 +1,32 @@ +# TODO: Consider upstreaming this as an extension to FITSWCS.jl """ - wcs_to_celestial_frame(wcs::WCSTransform) + celestial_frame(wcs::WCSTransform) + celestial_frame(wcs::SlicedWCSTransform) -Returns the reference frame of a WCSTransform. -The reference frame supported in Julia are FK5, ICRS and Galactic. -""" -function wcs_to_celestial_frame(wcs::WCSTransform) - radesys = wcs.radesys +The `SkyCoords` coordinate type corresponding to the celestial reference frame of a WCS transform: `ICRSCoords`, `FK5Coords{equinox}`, or `GalCoords`. - xcoord = wcs.ctype[1][1:4] - ycoord = wcs.ctype[2][1:4] +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. - if radesys == "" - if xcoord == "GLON" && ycoord == "GLAT" - radesys = "Gal" - elseif xcoord == "TLON" && ycoord == "TLAT" - radesys = "ITRS" - end +Throws an `ArgumentError` for frames without a `SkyCoords` counterpart (e.g., FK4 or ITRS). +""" +function celestial_frame(wcs::WCSTransform) + radesys = wcs.radesys + 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 + frame = isempty(radesys) ? join(wcs.ctype, ", ") : radesys + throw(ArgumentError("WCS celestial frame ($frame) is not supported by SkyCoords")) +end - return radesys +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/Project.toml b/test/Project.toml index 221b475..b308009 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,12 +1,19 @@ [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" +SkyCoords = "fc659fc5-75a3-5475-a2ea-3da92c065361" 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"} diff --git a/test/core.jl b/test/core.jl index 72e1b8d..61bc9c5 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,24 @@ @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") + # 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)) - 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..9c0d726 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,58 +42,68 @@ 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, [Card("EXTNAME", "SCI"); 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 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; hdu_out = "SCI") + @test result[1] isa WCSTransform + @test result[2] isa Tuple + end 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) + 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 = 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..cab1fc4 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,31 +1,44 @@ -using Reproject: wcs_to_celestial_frame -@testset "wcs to celestial frame" begin - wcs1 = WCSTransform(2; - ctype = ["RA---AIR", "DEC--AIR"], - ) - wcs2 = WCSTransform(2; - ctype = ["RA---AIR", "DEC--AIR"], - equinox = 1888.67 - ) - wcs3 = WCSTransform(2; - ctype = ["RA---AIR", "DEC--AIR"], - equinox = 2000 - ) - wcs4 = WCSTransform(2; - ctype = ["GLON--", "GLAT--"], - ) - wcs5 = WCSTransform(2; - ctype = ["TLON", "TLAT"], - ) - wcs6 = WCSTransform(2; - ctype = ["RA---AIR", "DEC--AIR"], - radesys = "UNK" - ) +using Reproject: celestial_frame +using SkyCoords: ICRSCoords, GalCoords, FK5Coords - @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" +@testset "celestial frame from WCS" begin + wcs1 = WCS(2; + ctype = ["RA---AIR", "DEC--AIR"], + ) + wcs2 = WCS(2; + ctype = ["RA---AIR", "DEC--AIR"], + equinox = 1888.67 + ) + wcs3 = WCS(2; + ctype = ["RA---AIR", "DEC--AIR"], + equinox = 2000 + ) + wcs4 = WCS(2; + ctype = ["GLON--", "GLAT--"], + ) + wcs5 = WCS(2; + ctype = ["TLON", "TLAT"], + ) + wcs6 = WCS(2; + ctype = ["RA---AIR", "DEC--AIR"], + radesys = "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 + +@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