Merge 2e2a5fa into ed4c2d2

.travis.yml

@@ -101,11 +101,11 @@ matrix:
         # time.
 
         - os: linux
-          env: PYTHON_VERSION=3.5 NUMPY_VERSION=1.13 ASTROPY_VERSION=3.0
+          env: PYTHON_VERSION=3.5 NUMPY_VERSION=1.16 ASTROPY_VERSION=3.0
         - os: linux
-          env: PYTHON_VERSION=3.6 NUMPY_VERSION=1.16
+          env: PYTHON_VERSION=3.6 NUMPY_VERSION=1.17
         - os: linux
-          env: NUMPY_VERSION=1.17
+          env: NUMPY_VERSION=1.18
 
         # Try numpy pre-release
         - os: linux

CHANGES.rst

@@ -1,6 +1,7 @@
 0.6 (unreleased)
 ----------------
 - Update wcs fitting to match Astropy (and use Astropy when available) [#29]
+- Limit the number of pixels used for WCS fitting to 100 [#30]
 
 
 0.5 (2020-01-13)

astrocut/cube_cut.py

@@ -2,30 +2,29 @@
 
 """This module implements the cutout functionality."""
 
+import os
+import warnings
+
+from time import time
+from itertools import product
+
 import numpy as np
 import astropy.units as u
 
 from astropy.io import fits
 from astropy.coordinates import SkyCoord
 from astropy import wcs
 
-from itertools import product
+from . import __version__ 
+from .exceptions import InputWarning, TypeWarning, InvalidQueryError
 
 # Note: Use the astropy function if available
 import astropy
-if astropy.utils.minversion(astropy,"4.0.1"):
+if astropy.utils.minversion(astropy, "4.0.2"):
     from astropy.wcs.utils import fit_wcs_from_points
 else:
     from .utils.wcs_fitting import fit_wcs_from_points
 
-from time import time
-
-import os
-import warnings
-
-from . import __version__
-from .exceptions import InputWarning, TypeWarning, InvalidQueryError
-
 
 class CutoutFactory():
     """
@@ -272,17 +271,49 @@ def _fit_cutout_wcs(self, cutout_wcs, cutout_shape):
         """
 
         # Getting matched pixel, world coordinate pairs
-        pix_inds = np.array(list(product(list(range(cutout_shape[1])), list(range(cutout_shape[0])))))
+        # We will choose no more than 100 pixels spread evenly throughout the image
+        # Centered on the center pixel.
+        # To do this we the appropriate "step size" between pixel coordinates
+        # (i.e. we take every ith pixel in each row/column) [TODOOOOOO]
+        # For example in a 5x7 cutout with i = 2 we get:
+        #
+        # xxoxoxx
+        # ooooooo
+        # xxoxoxx
+        # ooooooo
+        # xxoxoxx
+        #
+        # Where x denotes the indexes that will be used in the fit.
+        y, x = cutout_shape
+        i = 1
+        while (x/i)*(y/i) > 100:
+            i += 1
+
+        xvals = list(reversed(range(x//2, -1, -i)))[:-1] + list(range(x//2, x, i))
+        if xvals[-1] != x-1:
+            xvals += [x-1]
+        if xvals[0] != 0:
+            xvals = [0] + xvals
+
+        yvals = list(reversed(range(y//2, -1, -i)))[:-1] + list(range(y//2, y, i))
+        if yvals[-1] != y-1:
+            yvals += [y-1]
+        if yvals[0] != 0:
+            yvals = [0] + yvals
+
+        pix_inds = np.array(list(product(xvals, yvals)))
         world_pix = SkyCoord(cutout_wcs.all_pix2world(pix_inds, 0), unit='deg')
 
         # Getting the fit WCS
-        linear_wcs = fit_wcs_from_points([pix_inds[:, 0], pix_inds[:, 1]], world_pix, proj_point=self.center_coord)
+        linear_wcs = fit_wcs_from_points([pix_inds[:, 0], pix_inds[:, 1]], world_pix, proj_point='center')
 
         self.cutout_wcs = linear_wcs
 
-        # Checking the fit
+        # Checking the fit (we want to use all of the pixels for this)
+        pix_inds = np.array(list(product(list(range(cutout_shape[1])), list(range(cutout_shape[0])))))
+        world_pix = SkyCoord(cutout_wcs.all_pix2world(pix_inds, 0), unit='deg')
         world_pix_new = SkyCoord(linear_wcs.all_pix2world(pix_inds, 0), unit='deg')
-    
+
         dists = world_pix.separation(world_pix_new).to('deg')
         sigma = np.sqrt(sum(dists.value**2))
 

astrocut/tests/test_cube_cut.py

@@ -218,11 +218,11 @@ def test_cutout_extras(tmpdir):
     # Test _fit_cutout_wcs #
     ########################
     max_dist, sigma = myfactory._fit_cutout_wcs(cutout_wcs_full, (3, 5))
-    assert round(max_dist.deg, 7) == 7e-07
-    assert round(sigma, 7) == 1.4e-06
+    assert max_dist.deg < 1e-05
+    assert sigma < 1e-05
 
     cry, crx = myfactory.cutout_wcs.wcs.crpix
-    assert round(cry) == 4
+    assert round(cry) == 3
     assert round(crx) == 2
 
 

astrocut/utils/wcs_fitting.py

@@ -166,8 +166,8 @@ def fit_wcs_from_points(xy, world_coords, proj_point='center',
     from astropy.wcs import Sip
     from scipy.optimize import least_squares
 
+    xp, yp = xy
     try:
-        xp, yp = xy
         lon, lat = world_coords.data.lon.deg, world_coords.data.lat.deg
     except AttributeError:
         unit_sph =  world_coords.unit_spherical
@@ -207,7 +207,7 @@ def fit_wcs_from_points(xy, world_coords, proj_point='center',
         raise ValueError("sip_degree must be None, or integer.")
 
     # set pixel_shape to span of input points
-    wcs.pixel_shape = (xp.max()-xp.min(), yp.max()-yp.min())
+    wcs.pixel_shape = (xp.max()+1-xp.min(), yp.max()+1-yp.min())
 
     # determine CRVAL from input
     close = lambda l, p: p[np.argmin(np.abs(l))]
@@ -229,8 +229,15 @@ def fit_wcs_from_points(xy, world_coords, proj_point='center',
     # use (1, 0, 0, 1) as initial guess, in case input wcs was passed in
     # and cd terms are way off.
     p0 = np.concatenate([wcs.wcs.cd.flatten(), wcs.wcs.crpix.flatten()])
+
+    xpmin, xpmax, ypmin, ypmax = xp.min(), xp.max(), yp.min(), yp.max()
+    if xpmin==xpmax: xpmin, xpmax = xpmin-0.5, xpmax+0.5
+    if ypmin==ypmax: ypmin, ypmax = ypmin-0.5, ypmax+0.5
+
     fit = least_squares(_linear_wcs_fit, p0,
-                        args=(lon, lat, xp, yp, wcs))
+                        args=(lon, lat, xp, yp, wcs),
+                        bounds=[[-np.inf,-np.inf,-np.inf,-np.inf, xpmin, ypmin],
+                                [ np.inf, np.inf, np.inf, np.inf, xpmax, ypmax]])
     wcs.wcs.crpix = np.array(fit.x[4:6])
     wcs.wcs.cd = np.array(fit.x[0:4].reshape((2, 2)))