Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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)
@aquatiko

This comment has been minimized.

Copy link
Owner Author

@aquatiko aquatiko commented Jun 3, 2019

convert_pixel_coordinate needs some optimisation.

@aquatiko

This comment has been minimized.

Copy link
Owner Author

@aquatiko aquatiko commented Jun 3, 2019

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.

@aquatiko

This comment has been minimized.

Copy link
Owner Author

@aquatiko aquatiko commented Jun 3, 2019

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.

@aquatiko

This comment has been minimized.

Copy link
Owner Author

@aquatiko aquatiko commented Jun 9, 2019

Done!!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.