covert_pixel_coordinate.jl

using SkyCoords: lon, lat
using FITSIO, WCS, SkyCoords, Interpolations
function wcs_to_celestial_frame(wcs::WCSTransform)
   if wcs.latpole == -1 && wcs.lonpole == -1              # or lat/lng ?
        return nothing
    end
    radesys = wcs.radesys
    
    if !isnan(wcs.equinox)
        equinox = wcs.equinox
    else
        equinox = nothing
    end
    
    xcoord = wcs.ctype[1][1:4]
    ycoord = wcs.ctype[2][1:4]
    
    if radesys == ""
        if xcoord == "RA--" && ycoord == "DEC-"
            if equinox === nothing
                radesys = "ICRS"
            elseif equinox < 1984.0
                radesys = "FK4"
            else
                radesys = "FK5"
            end  
        elseif xcoord == "GLON" && ycoord == "GLAT"
            radesys = "Gal"
        elseif xcoord == "TLON" && ycoord == "TLAT"
            radesys = "ITRS"
        end
    end
    
    return radesys
end

function convert_pixel_coordinates(input_data, wcs_in, wcs_out,shape_out)
    type_in = wcs_to_celestial_frame(wcs_in)
    type_out = wcs_to_celestial_frame(wcs_out)
    
    if type_in == type_out
        return input_data
    end
    
    img_out = fill!(zeros(shape_out), NaN)
    itp = interpolate(input_data, BSpline(Quadratic(Reflect(OnCell()))))
    
    for i in 1:shape_out[1]
        for j in 1:shape_out[2]
            pix_cord_in = [float(i); float(j)]
            world_cord_in = pix_to_world(wcs_in, pix_cord_in)

            if type_in == "ICRS"
                cord_in = ICRSCoords(deg2rad(world_cord_in[1]), deg2rad(world_cord_in[2]))
            elseif type_in == "Gal"
                cord_in = GalCoords(deg2rad(world_cord_in[1]), deg2rad(world_cord_in[2]))
            elseif type_in == "FK5"
                cord_in = FK5Coords{wcs_in.equinox}(deg2rad(world_cord_in[1]), deg2rad(world_cord_in[2]))
            else
                throw(ArgumentError("Unsupported input WCS coordinate type"))
            end

            if type_out == "ICRS"
                cord_out = convert(ICRSCoords, cord_in)
            elseif type_out == "Gal"
                cord_out = convert(GalCoords, cord_in)
            elseif type_out == "FK5"
                cord_out = convert(FK5Coords{wcs_out.equinox}, cord_in)
            else
                throw(ArgumentError("Unsupported output WCS coordinate type"))
            end

            pix_cord_out = world_to_pix(wcs_out, [rad2deg(lon(cord_out)), rad2deg(lat(cord_out))])

            if 1 <= pix_cord_out[1] <= size(input_data)[1] && 1 <= pix_cord_out[2] <= size(input_data)[2] 
                img_out[i,j] = itp(pix_cord_out[1], pix_cord_out[2])
            end
        end
    end
    
    return img_out
end
imgin = FITS("original1.fits")         # project this
imgout = FITS("original2.fits")        # into this coordinate
wcs_in = WCS.from_header(read_header(imgin[1], String))[1]
wcs_out = WCS.from_header(read_header(imgout[1], String))[1]
img = read(imgout[1])
shape_out = size(read(imgin[1]))
tff = convert_pixel_coordinates(img, wcs_in, wcs_out, shape_out)

convert_pixel_coordinate needs some optimisation.

What convert_pixel_coordinate does:
It converts all the Cartesian coordinates of input_data like (1,1) , (1,2), (2,1), (5,6) etc. from frame of wcs_in to wcs_out.
The output should be a mapping of say pixel value at (1,1) with converted pixel coordinate of (1,1) in wcs_out frame.
wcs_to_celestial_frame gives the frame of wcs given.

I'm planning to set up an interpolator using x = result[1], result[2] and y = result[3]
and then interpolate for all integer coordinates for output shape_out to form the reprojected image.

Done!!