JuliaAstro / FITSIO.jl

Flexible Image Transport System (FITS) file support for Julia
http://juliaastro.org/FITSIO.jl/
MIT License
55 stars 29 forks source link

Modifying an image in a FITS file #116

Closed evgenyneu closed 5 years ago

evgenyneu commented 5 years ago

Is there a way to modify image data in a .fits file without touching anything else?

I have a .fits file that contains image data as well as several extensions (tables). I want to replace the image data while keeping everything else (headers/extensions etc.) unchanged. Is there a way to do it?

Right now, I'm creating a new file and copying all the extensions to it:

hdul_data = FITS("input.fits")

# Load the image data
image_data = read(hdul_data[1])

# Modify data
image_data[10, 10] = 42

# Save the data array into the output fits file
data_out = FITS("output.fits", "w");

for hdu in hdul_data
  if typeof(hdu)==FITSIO.ImageHDU
    # Write the image to output file
    write(data_out, image_data; header=header_data)
  elseif typeof(hdu)==FITSIO.TableHDU
    # Copy the table into the output file
    table_data = Dict{String,Array}()

    for colname in FITSIO.colnames(hdu)
      col_data = read(hdu, colname)
      table_data[colname] = col_data
    end

    header = read_header(hdu)
    ext_name = header["EXTNAME"]
    write(data_out, table_data; hdutype=FITSIO.TableHDU, header=header, name=ext_name)
  end
end

This is not ideal, because I'm not sure if the new file will be exactly the same as the old one (except for the changed image). I would prefer to just modify the image data. These are example of how this could be done in Python and Fortran.

Python example:

from astropy.io import fits

hdul_data = fits.open("data.fits")

# Load the image data
image_data = hdul_data["PRIMARY"].data

# Change image data
image_data[5,5]=123

# Save the data array into the output fits file
hdul_data.writeto("output.fits")

Fortran example (cfitsio):

! Get an unused Logical Unit Number to use to open the FITS file.
call ftgiou(unit_data,status)

! Open the FITS file, with read-only access.
readwrite=0
call ftopen(unit_data,filename_data,readwrite,blocksize,status)

! Load the image data
group=1
nullval=0
call ftg2di(unit_data, group, nullval, max_x, naxes_data(1), naxes_data(2), image_data, anynull, status)

! Modify the image
image_data(10,10) = 42

! Create the output fits file
! -------------

! Get an unused Logical Unit Number to use to create the FITS file.
call ftgiou(unit_output, status)

! Create the output file
blocksize=1
call ftinit(unit_output, filename_output, blocksize, status)

! Copy the entire data file to the output file
call ftcpfl(unit_data, unit_output, 1, 1, 1, status)

! Move back to the primary array
call ftmahd(unit_output, 1, hdutype, status)

! Save the data array into the output fits file
group=1
nullval=0
call ftp2di(unit_output, group, max_x, naxes_data(1), naxes_data(2), image_data, status)
kbarbary commented 5 years ago

This should be possible by opening the file in read-write mode, and then writing the image data to the appropriate extension. However, it is not currently possible. We need to create a new method write(::ImageHDU, ::Array) to write to an existing Image HDU instead of creating a new one.

evgenyneu commented 5 years ago

@kbarbary thank you, that's good to know. Can I attempt to submit a pull request with this addition?

kbarbary commented 5 years ago

That would be great! If you need a starting place, look at the write(::FITS, ::Array, ...) method in image.jl. the implementation should be similar, just without the creation of the new extension.

giordano commented 5 years ago

Fixed by #117