Reproject-core (#3)

* reproject-core function added

* docs and logic fixes

* CI error fixes

* incresing coverage

Author: aquatiko

Project.toml

@@ -5,6 +5,8 @@ version = "0.1.0"
 
 [deps]
 FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+SkyCoords = "fc659fc5-75a3-5475-a2ea-3da92c065361"
 WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d"
 
 [compat]

src/Reproject.jl

@@ -1,8 +1,13 @@
 module Reproject
 
-using FITSIO, WCS
+using FITSIO, WCS, Interpolations, SkyCoords
+using SkyCoords: lat,lon
 
 include("parsers.jl")
 include("utils.jl")
+include("core.jl")
+
+export
+    reproject
 
 end # module

src/core.jl

@@ -0,0 +1,96 @@
+"""
+    reproject(input_data, output_projection; shape_out = nothing, order = 1, hdu_in = 1, hdu_out = 1)
+
+Reprojects image data to a new projection using interpolation.
+
+# Arguments
+- `input_data`: Image data which is being reprojected.
+                It can be an ImageHDU, FITS object or name of a FITS file.
+- `output_projection`: Frame in which data is reprojected.
+                       Frame can be taken from WCSTransform object, ImageHDU, FITS or name of FITS file.
+- `shape_out`: Shape of image after reprojection.
+- `order`: Order of interpolation.
+           0: Nearest-neighbor
+           1: Linear
+           2: Quadratic
+- `hdu_in`: Used to specify HDU number when giving input as FITS or name of FITS file.
+- `hud_out:` Used to specify HDU number when giving output projection as FITS or name of FITS file.
+"""
+function reproject(input_data, output_projection; shape_out = nothing, order::Int = 1, hdu_in::Int = 1, hdu_out::Int = 1)
+    if input_data isa ImageHDU
+        array_in, wcs_out = parse_input_data(input_data)
+    else
+        array_in, wcs_out = parse_input_data(input_data, hdu_in)
+    end
+    
+    if output_projection isa FITS || output_projection isa String
+        wcs_in, shape_out = parse_output_projection(output_projection, hdu_out)
+    else
+        wcs_in, shape_out = parse_output_projection(output_projection, shape_out)
+    end
+    
+    type_in = wcs_to_celestial_frame(wcs_in)
+    type_out = wcs_to_celestial_frame(wcs_out)
+    
+    if type_in == type_out && shape_out === nothing
+        return array_in
+    end
+    
+    img_out = fill(NaN, shape_out)
+    
+    itp = interpolator(array_in, order)
+    shape_in = size(array_in)
+    
+    for i in 1:shape_out[1]
+        for j in 1:shape_out[2]
+            pix_coord_in = [float(i), float(j)]
+            world_coord_in = pix_to_world(wcs_in, pix_coord_in)
+
+            if type_in == "ICRS"
+                coord_in = ICRSCoords(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2]))
+            elseif type_in == "Gal"
+                coord_in = GalCoords(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2]))
+            elseif type_in == "FK5"
+                coord_in = FK5Coords{wcs_in.equinox}(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2]))
+            else
+                throw(ArgumentError("Unsupported output WCS coordinate type"))
+            end
+
+            if type_out == "ICRS"
+                coord_out = convert(ICRSCoords, coord_in)
+            elseif type_out == "Gal"
+                coord_out = convert(GalCoords, coord_in)
+            elseif type_out == "FK5"
+                coord_out = convert(FK5Coords{wcs_out.equinox}, coord_in)
+            else
+                throw(ArgumentError("Unsupported input WCS coordinate type"))
+            end
+
+            pix_coord_out = world_to_pix(wcs_out, [rad2deg(lon(coord_out)), rad2deg(lat(coord_out))])
+
+            if 1 <= pix_coord_out[1] <= shape_in[1] && 1 <= pix_coord_out[2] <= shape_in[2] 
+                img_out[i,j] = itp(pix_coord_out[1], pix_coord_out[2])
+            end
+        end
+    end
+    
+    return img_out, (!isnan).(img_out)
+end
+
+
+"""
+    interpolator(array_in, order::Int)
+
+Returns an interpolator with the given array and order of interpolation.
+"""
+function interpolator(array_in::AbstractArray, order::Int)
+   if order == 0
+        itp = interpolate(array_in, BSpline(Constant()))
+    elseif order == 1
+        itp = interpolate(array_in, BSpline(Linear()))
+    else
+        itp = interpolate(array_in, BSpline(Quadratic(InPlace(OnCell()))))
+    end
+    
+    return itp
+end

test/core.jl

@@ -0,0 +1,59 @@
+@testset "reproject-core" begin
+    if !isfile("data/gc_2mass_k.fits")
+        mkpath("data")
+        download("https://astropy.stsci.edu/data/galactic_center/gc_2mass_k.fits", "data/gc_2mass_k.fits")
+    end
+    if !isfile("data/gc_msx_e.fits")
+        mkpath("data")
+        download("https://astropy.stsci.edu/data/galactic_center/gc_msx_e.fits", "data/gc_msx_e.fits")
+    end
+    
+    imgin = FITS("data/gc_msx_e.fits")         # project this
+    imgout = FITS("data/gc_2mass_k.fits")        # into this coordinate
+    
+    hdu1 = astropy.io.fits.open("data/gc_2mass_k.fits")[1]
+    hdu2 = astropy.io.fits.open("data/gc_msx_e.fits")[1]
+    
+    function check_diff(arr1, arr2)
+        diff = 1e-3
+        count = 0
+        for i in 1 : size(arr1)[1] - 1
+            for j in 1 : size(arr2)[2] - 1
+                if abs(arr1[i,j] - arr2[j,i]) > diff
+                    count+=1
+                end
+            end
+        end
+        return  iszero(count)
+    end
+    
+    @test check_diff(reproject(imgin, imgout, order = 0)[1], rp.reproject_interp(hdu2, hdu1.header, order = 0)[1]) == true
+    @test check_diff(reproject(imgout, imgin, order = 0)[1], rp.reproject_interp(hdu1, hdu2.header, order = 0)[1]) == true
+    @test check_diff(reproject(imgin, imgout, order = 1)[1], rp.reproject_interp(hdu2, hdu1.header, order = 1)[1]) == true
+    @test check_diff(reproject(imgin, imgout, order = 2)[1], rp.reproject_interp(hdu2, hdu1.header, order = 2)[1]) == true
+    @test check_diff(reproject(imgin[1], imgout[1], shape_out = (1000,1000))[1],
+            rp.reproject_interp(hdu2, astropy.wcs.WCS(hdu1.header), shape_out = (1000,1000))[1]) == true
+    
+    wcs = WCSTransform(2; ctype = ["RA---AIR", "DEC--AIR"], radesys = "UNK")
+    @test_throws ArgumentError reproject(imgin, wcs, shape_out = (100,100))
+
+    fname = tempname() * ".fits"
+    f = FITS(fname, "w")
+    inhdr = FITSHeader(["CTYPE1", "CTYPE2", "RADESYS", "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT",
+                        "HISTORY"],
+                       ["RA---TAN", "DEC--TAN", "UNK", 1.0, 1, true, "string value", nothing, nothing],
+                       ["",
+                        "",
+                        "",
+                        "floating point keyword",
+                        "",
+                        "boolean keyword",
+                        "string value",
+                        "this is a comment",
+                        "this is a history"])
+
+    indata = reshape(Float32[1:100;], 5, 20)
+    write(f, indata; header=inhdr)
+    @test_throws ArgumentError reproject(f[1], imgin, shape_out = (100,100))
+    close(f)
+end

test/runtests.jl

@@ -7,8 +7,10 @@ ENV["PYTHON"]=""
 Conda.add_channel("astropy")
 Conda.add("reproject")
 rp = pyimport("reproject")
+astropy = pyimport("astropy")
 
 @testset "Reproject.jl" begin
     include("parsers.jl")
     include("utils.jl")
+    include("core.jl")
 end