Interpolator tweak as tweaked map_coordinates

Author: aquatiko

src/core.jl

@@ -37,7 +37,7 @@ function reproject(input_data, output_projection; shape_out = nothing, order::In
     end
 
     img_out = fill(NaN, shape_out)
-
+    array_in = pad_edges(array_in)
     itp = interpolator(array_in, order)
     shape_in = size(array_in)
 
@@ -68,8 +68,8 @@ function reproject(input_data, output_projection; shape_out = nothing, order::In
 
             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])
+            if 0.5 <= pix_coord_out[1] <= shape_in[1] - 1 && 0.5 <= pix_coord_out[2] <= shape_in[2] - 1
+                img_out[i,j] = itp(pix_coord_out[1] + 1, pix_coord_out[2] + 1)
             end
         end
     end

src/utils.jl

@@ -20,3 +20,18 @@ function wcs_to_celestial_frame(wcs::WCSTransform)
 
     return radesys
 end
+
+"""
+    pad_edges(array_in::Matrix{T}) where {T}
+
+Pads a given array and creates a border with edge elements.
+"""
+function pad_edges(array_in::Matrix{T}) where {T}
+    image = zeros(size(array_in)[1] + 2, size(array_in)[2] + 2)
+    image[2:end-1,2:end-1] = array_in
+    image[2:end-1,1] = array_in[:,1]
+    image[2:end-1,end] = array_in[:,end]
+    image[1,:] = image[2,:]
+    image[end,:] = image[end-1,:]
+    return image
+end

test/core.jl

@@ -25,25 +25,12 @@ end
     hdu1 = astropy.io.fits.open(joinpath("data", "gc_2mass_k.fits"))[1]
     hdu2 = astropy.io.fits.open(joinpath("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
+    @test isapprox(reproject(imgin, imgout, order = 0)[1]', rp.reproject_interp(hdu2, hdu1.header, order = 0)[1], nans = true, atol = 1e-4)
+    @test isapprox(reproject(imgout, imgin, order = 0)[1]', rp.reproject_interp(hdu1, hdu2.header, order = 0)[1], nans = true, atol = 1e-3)
+    @test isapprox(reproject(imgin, imgout, order = 1)[1]', rp.reproject_interp(hdu2, hdu1.header, order = 1)[1], nans = true, atol = 1e-4)
+    @test isapprox(reproject(imgin, imgout, order = 2)[1]', rp.reproject_interp(hdu2, hdu1.header, order = 2)[1], nans = true, atol = 1e-4)
+    @test isapprox(reproject(imgin[1], imgout[1], shape_out = (1000,1000))[1]',
+            rp.reproject_interp(hdu2, astropy.wcs.WCS(hdu1.header), shape_out = (1000,1000))[1], nans = true, atol = 1e-4)
 
     wcs = WCSTransform(2; ctype = ["RA---AIR", "DEC--AIR"], radesys = "UNK")
     @test_throws ArgumentError reproject(imgin, wcs, shape_out = (100,100))