Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# aquatiko/covert_pixel_coordinate.jl

Last active Jun 9, 2019
 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 commented Jun 3, 2019 • edited

 `convert_pixel_coordinate` needs some optimisation.

### 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 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 commented Jun 9, 2019

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