Add round-tripping to coordinate conversion (need to optimize later)

astrofrog · Aug 21, 2014 · 55a18da · 55a18da
1 parent f9d7c62
commit 55a18da
Showing 1 changed file with 15 additions and 0 deletions.
diff --git a/reproject/interpolation/core.py b/reproject/interpolation/core.py
@@ -37,6 +37,21 @@ def get_input_pixels_celestial(wcs_in, wcs_out, shape_out):
 
     xp_in, yp_in = wcs_in.wcs_world2pix(xw_in, yw_in, 0)
 
+    # Now convert back to check that coordinates round-trip, if not then set to NaN
+
+    xw_in, yw_in = wcs_in.wcs_pix2world(xp_in, yp_in, 0)
+
+    xw_out, yw_out = convert_world_coordinates(xw_in, yw_in, wcs_in, wcs_out)
+
+    xp_out_check, yp_out_check = wcs_out.wcs_world2pix(xw_out, yw_out, 0)
+
+    print(xp_out, xp_out_check, yp_out, yp_out_check)
+
+    reset = (np.abs(xp_out - xp_out_check) > 1) | (np.abs(yp_out - yp_out_check) > 1)
+
+    xp_in[reset] = np.nan
+    yp_in[reset] = np.nan
+
     return xp_in, yp_in