search

emmt / EasyFITS.jl

Using FITS files made easier for Julia
Other
8 stars 2 forks source link

Problem when writing old header with new data (structural keywords are modified) #6

Open buthanoid opened 1 year ago

buthanoid commented 1 year ago

Say I have a FITS file with some good information in its header. Say I want to create a new FITS file with this good header, but along with new data of a different shape. It breaks because the good header changes the HDU size:

julia> using EasyFITS
julia> myoldheader = FitsHeader("NAXIS" => 1, "NAXIS1" => 8, "THIS_IS" => "A Good Information")
julia> mynewdata ::Matrix{Float32} = [ 1; 2 ;; 3; 4 ]
julia> mynewfile = openfits("/tmp/toto.fits", "w!")
julia> write(mynewfile, myoldheader, mynewdata)
ERROR: number of dimensions has changed, rebuild the HDU
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] get_img_size(hdu::FitsImageHDU{Float32, 2})
   @ EasyFITS ~/Documents/EasyFITS.jl-antoine/src/images.jl:32
 [3] #write#75
   @ ~/Documents/EasyFITS.jl-antoine/src/images.jl:404 [inlined]
 [4] write
   @ ~/Documents/EasyFITS.jl-antoine/src/images.jl:398 [inlined]
 [5] write(file::FitsFile, hdr::FitsHeader, arr::Matrix{Float32})
   @ EasyFITS ~/Documents/EasyFITS.jl-antoine/src/images.jl:365

It happens from this function:

function write(file::FitsFile, hdr::Union{Header,Nothing},
               arr::AbstractArray{T,N}) where {T<:Number,N}
    write(merge!(write(file, FitsImageHDU, T, size(arr)), hdr), arr)
    return file # returns the file not the HDU
end

First an image HDU is created of size 2x2. Then the header is copied, including the NAXIS keywords which reshape the HDU. When we finally attempt to write the data, it fails because our hdu has type FitsImageHDU{Float32,2} but possess a NAXIS = 1 keyword. Let's see it step by step:

julia> using EasyFITS
julia> myoldheader = FitsHeader("NAXIS" => 1, "NAXIS1" => 8, "THIS_IS" => "A Good Information")
julia> mynewdata ::Matrix{Float32} = [ 1; 2 ;; 3; 4 ]
julia> mynewfile = openfits("/tmp/toto.fits", "w!")
julia> mynewhdu = write(mynewfile, FitsImageHDU, Float32, size(mynewdata))
FitsImageHDU{Float32,2}: num = 1, size = 2×2
julia> merge!(mynewhdu, myoldheader)
FitsImageHDU{Float32,2}: num = 1, size = 8
julia> write(mynewhdu, mynewdata)
ERROR: number of dimensions has changed, rebuild the HDU
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] get_img_size(hdu::FitsImageHDU{Float32, 2})
   @ EasyFITS ~/Documents/EasyFITS.jl-antoine/src/images.jl:32
 [3] #write#75
   @ ~/Documents/EasyFITS.jl-antoine/src/images.jl:404 [inlined]
 [4] write(hdu::FitsImageHDU{Float32, 2}, arr::Matrix{Float32})
   @ EasyFITS ~/Documents/EasyFITS.jl-antoine/src/images.jl:398

The write fails by calling this function:

julia> EasyFITS.get_img_size(mynewhdu)
ERROR: number of dimensions has changed, rebuild the HDU

This function compares CFITSIO.fits_get_img_dim with the type argument N. The first gives 1 and the second 2.\ As the number of dimensions is part of the type, and do not change when we change NAXIS, it is wrong to change NAXIS. Some silent wrong behaviour happens, see:

julia> mynewhdu.data_ndims
2
julia> mynewhdu["NAXIS"]
FitsCard: NAXIS   = 1

Whereas the function get_img_size fails, the property data_ndims still gives 2 without warning.

I propose to at least use the new BaseFITS.Parser.is_structural when merging a header in an existing HDU. And to at least print a warning when a structural keyword is written. Maybe with a keyword "nocheck" to remove the checking when performance is needed.\ For now we cannot offer a structural update by the keywords, as long as we keep the old dimension in the type of the HDU.\ The warning could direct the user to the documentation where an example of filtering the structural keywords would be found. Like this:

using EasyFITS, BaseFITS
function filter_structural_keywords(fitsheader::FitsHeader) ::FitsHeader
    return filter(!BaseFITS.Parser.is_structural, fitsheader)
end