From e9a757d540418324f6b0a5a30982d693d0d3f477 Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Wed, 12 Jun 2019 16:05:24 -0400
Subject: [PATCH 01/15] rough cut adding wcs fitting

---
 astrocut/cube_cut.py          | 171 ++++++++++++++++++++++++++++-
 astrocut/utils/wcs_fitting.py | 197 ++++++++++++++++++++++++++++++++++
 2 files changed, 366 insertions(+), 2 deletions(-)
 create mode 100644 astrocut/utils/wcs_fitting.py

diff --git a/astrocut/cube_cut.py b/astrocut/cube_cut.py
index 86b899f9..981c116d 100644
--- a/astrocut/cube_cut.py
+++ b/astrocut/cube_cut.py
@@ -9,6 +9,13 @@
 from astropy.coordinates import SkyCoord
 from astropy import wcs
 
+from itertools import product
+
+try:
+    from astropy.utils import fit_wcs_from_points
+except: # astropy version does have the function
+    from .utils.wcs_fitting import fit_wcs_from_points
+
 from time import time
 
 import os
@@ -33,6 +40,7 @@ def __init__(self):
         """
 
         self.cube_wcs = None  # WCS information from the image cube
+        self.cutout_wcs = None # WCS information (linear) for the cutout
         self.cutout_lims = np.zeros((2, 2), dtype=int)  # Cutout pixel limits, [[ymin,ymax],[xmin,xmax]]
         self.center_coord = None  # Central skycoord
         
@@ -189,6 +197,156 @@ def _get_cutout_limits(self, cutout_size):
         self.cutout_lims = lims
 
 
+    def _fit_cutout_wcs(self, cutout_wcs, cutout_shape):
+        """
+        Given a full (including SIP coefficients) wcs for the cutout, 
+        calculate the best fit linear wcs, and a measure of the goodness-of-fit.
+        
+        The news WCS is stored in ``self.cutout_wcs``.
+
+        Parameters
+        ----------
+        cutout_wcs : 
+            The full (including SIP coefficients) cutout WCS object 
+        cutout_shape : tuple
+            The shape of the cutout in the form (width, height).
+
+        Returns
+        -------
+        response : tuple
+            Goodness-of-fit statistics. (max dist, sigma)
+        """
+
+        # Getting matched pixel, world coordinate pairs
+        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, 1), unit='deg')
+
+        # Getting the fit WCS
+        linear_wcs = fit_wcs_from_points(pix_inds[:,0], pix_inds[:,1], world_pix, mode='wcs', proj_point="center")
+
+        self.cutout_wcs = linear_wcs
+
+        # Checking the fit
+        #world_pix_nosip = SkyCoord(cutout_wcs.wcs_pix2world(pix_inds,1), unit='deg')
+        world_pix_new = SkyCoord(linear_wcs.all_pix2world(pix_inds,1), unit='deg')
+    
+        dists = world_pix.separation(world_pix_new)
+        #dists_nosip = world_pix.separation(world_pix_nosip)
+    
+        sigma = np.sqrt(sum(dists.value**2))
+        #nosip_sigma = np.sqrt(sum(dists_nosip.value**2))
+
+        # TODO: If just wacking off the SIP coefficients is better than the fit, just do that
+        #if nosip_sigma < sigma:
+        #    pass
+
+        return (dists.max(), sigma)
+    
+
+    def _get_full_cutout_wcs(self, cube_table_header):
+        """
+        TODO: Document
+        """
+        
+        wcs_header = self.cube_wcs.to_header(relax=True)
+
+        # Using table comment rather than the default ones if available
+        for kwd in wcs_header:
+            if cube_table_header.get(kwd, None):
+                wcs_header.comments[kwd] = cube_table_header[kwd]
+
+        # Adjusting the CRPIX values
+        wcs_header["CRPIX1"] -= self.cutout_lims[0, 0]
+        wcs_header["CRPIX2"] -= self.cutout_lims[1, 0]
+
+        # Adding the physical wcs keywords
+        wcs_header.set("WCSNAMEP", "PHYSICAL", "name of world coordinate system alternate P")
+        wcs_header.set("WCSAXESP", 2, "number of WCS physical axes")
+    
+        wcs_header.set("CTYPE1P", "RAWX", "physical WCS axis 1 type CCD col")
+        wcs_header.set("CUNIT1P", "PIXEL", "physical WCS axis 1 unit")
+        wcs_header.set("CRPIX1P", 1, "reference CCD column")
+        wcs_header.set("CRVAL1P", self.cutout_lims[0, 0] + 1, "value at reference CCD column")
+        wcs_header.set("CDELT1P", 1.0, "physical WCS axis 1 step")
+                
+        wcs_header.set("CTYPE2P", "RAWY", "physical WCS axis 2 type CCD col")
+        wcs_header.set("CUNIT2P", "PIXEL", "physical WCS axis 2 unit")
+        wcs_header.set("CRPIX2P", 1, "reference CCD row")
+        wcs_header.set("CRVAL2P", self.cutout_lims[1, 0] + 1, "value at reference CCD row")
+        wcs_header.set("CDELT2P", 1.0, "physical WCS axis 2 step")
+
+        return wcs.WCS(wcs_header)
+
+    
+    def _get_cutout_wcs_dict(self):
+        """
+        Transform the cutout WCS object into the cutout column WCS keywords.
+        Adds the physical keywords for transformation back from cutout to location on FFI.
+        This is a very TESS specific function.
+        
+        Returns
+        -------
+        response: dict
+            Cutout wcs column header keywords as dictionary of 
+            ``{<kwd format string>: [value, desc]} pairs.``
+        """
+        wcs_header = self.cutout_wcs.to_header()
+
+        cutout_wcs_dict = dict()
+
+        
+        ## Cutout array keywords ##
+
+        cutout_wcs_dict["WCAX{}"] = [wcs_header['WCSAXES'], "number of WCS axes"]
+        # TODO: check for 2? this must be two
+
+        cutout_wcs_dict["1CTYP{}"] = [wcs_header["CTYPE1"], "right ascension coordinate type"]
+        cutout_wcs_dict["2CTYP{}"] = [wcs_header["CTYPE2"], "declination coordinate type"]
+        
+        cutout_wcs_dict["1CRPX{}"] = [wcs_header["CRPIX1"], "[pixel] reference pixel along image axis 1"]
+        cutout_wcs_dict["2CRPX{}"] = [wcs_header["CRPIX2"], "[pixel] reference pixel along image axis 2"]
+    
+        cutout_wcs_dict["1CRVL{}"] = [wcs_header["CRVAL1"], "[deg] right ascension at reference pixel"]
+        cutout_wcs_dict["2CRVL{}"] = [wcs_header["CRVAL2"], "[deg] declination at reference pixel"]
+
+        cutout_wcs_dict["1CUNI{}"] = [wcs_header["CUNIT1"], "physical unit in column dimension"]
+        cutout_wcs_dict["2CUNI{}"] = [wcs_header["CUNIT2"], "physical unit in row dimension"]
+
+        cutout_wcs_dict["1CDLT{}"] = [wcs_header["CDELT1"], "[deg] pixel scale in RA dimension"]
+        cutout_wcs_dict["2CDLT{}"] = [wcs_header["CDELT1"], "[deg] pixel scale in DEC dimension"]
+
+        cutout_wcs_dict["11PC{}"] = [wcs_header["PC1_1"], "Coordinate transformation matrix element"]
+        cutout_wcs_dict["12PC{}"] = [wcs_header["PC1_2"], "Coordinate transformation matrix element"]
+        cutout_wcs_dict["21PC{}"] = [wcs_header["PC2_1"], "Coordinate transformation matrix element"]
+        cutout_wcs_dict["22PC{}"] = [wcs_header["PC2_2"], "Coordinate transformation matrix element"]
+
+        
+        ## Physical keywords ##
+        # TODO: Make sure these are correct
+        cutout_wcs_dict["WCSN{}P"] = ["PHYSICAL", "table column WCS name"]
+        cutout_wcs_dict["WCAX{}P"] = [2, "table column physical WCS dimensions"]
+    
+        cutout_wcs_dict["1CTY{}P"] = ["RAWX", "table column physical WCS axis 1 type, CCD col"]
+        cutout_wcs_dict["2CTY{}P"] = ["RAWY", "table column physical WCS axis 2 type, CCD row"]
+    
+        cutout_wcs_dict["1CUN{}P"] = ["PIXEL", "table column physical WCS axis 1 unit"]
+        cutout_wcs_dict["2CUN{}P"] = ["PIXEL", "table column physical WCS axis 2 unit"]
+    
+        cutout_wcs_dict["1CRV{}P"] = [self.cutout_lims[0, 0] + 1,
+                                      "table column physical WCS ax 1 ref value"]
+        cutout_wcs_dict["2CRV{}P"] = [self.cutout_lims[1, 0] + 1,
+                                      "table column physical WCS ax 2 ref value"]
+
+        # TODO: can we calculate these? or are they fixed?
+        cutout_wcs_dict["1CDL{}P"] = [1.0, "table column physical WCS a1 step"]    
+        cutout_wcs_dict["2CDL{}P"] = [1.0, "table column physical WCS a2 step"]
+    
+        cutout_wcs_dict["1CRP{}P"] = [1, "table column physical WCS a1 reference"]
+        cutout_wcs_dict["2CRP{}P"] = [1, "table column physical WCS a2 reference"]
+
+        return cutout_wcs_dict
+
+    
     def _get_cutout_wcs(self):
         """
         Transform the cube WCS object into the cutout column WCS keywords.
@@ -585,7 +743,12 @@ def _build_tpf(self, cube_fits, img_cube, uncert_cube, cutout_wcs_dict, aperture
         # Building the aperture HDU
         aperture_hdu = fits.ImageHDU(data=aperture)
         aperture_hdu.header['EXTNAME'] = 'APERTURE'
-        self._add_aperture_wcs(aperture_hdu.header, cube_fits[2].header)
+        for kwd, val, cmt in self.cutout_wcs.to_header().cards: 
+            aperture_hdu.header.set(kwd, val, cmt)
+        #self._add_aperture_wcs(aperture_hdu.header, cube_fits[2].header)
+        # Adding extra aperture keywords (TESS specific)
+        aperture_hdu.header.set("NPIXMISS", None, "Number of op. aperture pixels not collected")
+        aperture_hdu.header.set("NPIXSAP", None, "Number of pixels in optimal aperture")
         aperture_hdu.header['INHERIT'] = True
     
         cutout_hdu_list = fits.HDUList([primary_hdu, table_hdu, aperture_hdu])
@@ -678,7 +841,11 @@ def cube_cut(self, cube_file, coordinates, cutout_size,
         img_cutout, uncert_cutout, aperture = self._get_cutout(cube[1].data, verbose=verbose)
 
         # Get cutout wcs info
-        cutout_wcs_dict = self._get_cutout_wcs()
+        #cutout_wcs_dict = self._get_cutout_wcs()
+        cutout_wcs_full = self._get_full_cutout_wcs(cube[2].header)
+        max_dist,sigma = self._fit_cutout_wcs(cutout_wcs_full, img_cutout.shape[1:])
+        # TODO: put the max_dist, sigma somewhere
+        cutout_wcs_dict = self._get_cutout_wcs_dict()
     
         # Build the TPF
         tpf_object = self._build_tpf(cube, img_cutout, uncert_cutout, cutout_wcs_dict, aperture)
diff --git a/astrocut/utils/wcs_fitting.py b/astrocut/utils/wcs_fitting.py
new file mode 100644
index 00000000..bbb3a70c
--- /dev/null
+++ b/astrocut/utils/wcs_fitting.py
@@ -0,0 +1,197 @@
+import numpy as np
+
+from astropy import units as u
+from astropy.wcs.utils import celestial_frame_to_wcs
+
+
+def _linear_transformation_fit(params, lon, lat, x, y, w_obj):
+    """ Objective function for fitting linear terms."""
+    pc = params[0:4]
+    crpix = params[4:6]
+
+    w_obj.wcs.pc = ((pc[0],pc[1]),(pc[2],pc[3]))
+    w_obj.wcs.crpix = crpix
+    lon2, lat2 = w_obj.wcs_pix2world(x,y,1)
+    resids = np.concatenate((lon-lon2, lat-lat2))
+
+    return resids
+
+
+def _sip_fit(params, lon, lat, u, v, w_obj, a_order, b_order, a_coeff_names,
+             b_coeff_names):
+        """ Objective function for fitting SIP."""
+        from astropy.modeling.models import SIP # here instead of top to avoid circular import
+
+        crpix = params[0:2]
+        cdx = params[2:6].reshape((2,2))
+        a_params, b_params = params[6:6+len(a_coeff_names)],params[6+len(a_coeff_names):]
+        w_obj.wcs.pc = cdx
+        w_obj.wcs.crpix = crpix
+        x,y = w_obj.wcs_world2pix(lon,lat,1) #'intermediate world coordinates', x & y
+        x,y = np.dot(w_obj.wcs.pc,(x - w_obj.wcs.crpix[0],y - w_obj.wcs.crpix[1]))
+
+        a_params, b_params = params[6:6+len(a_coeff_names)], params[6+len(a_coeff_names):]
+        a_coeff, b_coeff = {}, {}
+        for i in range(len(a_coeff_names)):
+            a_coeff['A_' + a_coeff_names[i]] = a_params[i]
+        for i in range(len(b_coeff_names)):
+            b_coeff['B_' + b_coeff_names[i]] = b_params[i]
+
+        sip = SIP(crpix=crpix, a_order=a_order, b_order=b_order, a_coeff=a_coeff, \
+                  b_coeff=b_coeff)
+        fuv, guv = sip(u,v)
+        xo, yo = np.dot(cdx, np.array([u+fuv-crpix[0], v+guv-crpix[1]]))
+        resids = np.concatenate((x-xo, y-yo))
+        return resids
+
+
+def fit_wcs_from_points(xp, yp, coords, mode, projection='TAN', proj_point=None, order=None, inwcs=None):
+    """Given a set of matched x,y pixel positions and celestial coordinates, solves for
+    WCS parameters and returns a WCS with these best fit values and other keywords based
+    on projection type, frame, and units.
+    Along with the matched coordinate pairs, users must provide the projection type
+    (e.g. 'TAN'), celestial coordinate pair for the projection point (or 'center' to use
+    the center of input sky coordinates), mode ('wcs' to fit only linear terms, or 'all'
+    to fit a SIP to the data). If a coordinate pair is passed to 'proj_point', it is
+    assumed to be in the same units respectivley as the input (lon, lat) coordinates.
+    Additionaly, if mode is set to 'all', the polynomial order must be provided.
+    Optionally, an existing ~astropy.wcs.WCS object with some may be passed in. For
+    example, this is useful if the user wishes to refit an existing WCS with better
+    astrometry, or are using a projection type with non-standard keywords. If any overlap
+    between keyword arguments passed to the function and values in the input WCS exist,
+    the keyword argument will override, including the CD/PC convention.
+    NAXIS1 and NAXIS2 will set as the span of the input points, unless an input WCS is
+    provided.
+    All output will be in degrees.
+    Parameters
+    ----------
+    xp, yp : float or `numpy.ndarray`
+        x & y pixel coordinates.
+    coords : `~astropy.coordinates.SkyCoord`
+        Skycoord object.
+    mode : str
+        Whether SIP distortion should be fit (``'all'``), or only the linear
+        terms of the core WCS transformation (``'wcs'``).
+    inwcs : None or `~astropy.wcs.WCS`
+        Optional input WCS object. Populated keyword values will be used. Keyword
+        arguments passed to the function will override any in 'inwcs' if both are present.
+    projection : None or str
+        Three-letter projection code. Defaults to TAN.
+    proj_point: `~astropy.coordinates.SkyCoord` or str
+        Celestial coordinates of projection point (lat, lon). If 'center', the geometric
+        center of input coordinates will be used for the projection. If None, and input
+        wcs with CRVAL set must be passed in.
+    order : None, int, or tuple of ints
+        A order, B order, respectivley, for SIP polynomial to be fit, or single integer
+        for both. Must be provided if mode is 'sip'.
+    Returns
+    -------
+    wcs : `~astropy.wcs.WCS`
+        WCS object with best fit coefficients given the matched set of X, Y pixel and
+        sky positions of sources.
+    Notes
+    -----
+    Please pay careful attention to the logic that applies when both an input WCS with
+    keywords populated is passed in. For example, if no 'CUNIT' is set in this input
+    WCS, the units of the CRVAL etc. are assumed to be the same as the input Skycoord.
+    Additionally, if 'RADESYS' is not set in the input WCS, this will be taken from the
+    Skycoord as well. """
+
+    from scipy.optimize import least_squares
+    from astropy.wcs import Sip
+    import copy
+
+    lon, lat = coords.data.lon.deg, coords.data.lat.deg
+    mode = mode.lower()
+    if mode not in ("wcs", "all"):
+        raise ValueError("mode must be 'wcs' or 'all'")
+    if inwcs is not None:
+        inwcs = copy.deepcopy(inwcs)
+    if projection is None:
+        if (inwcs is None) or ('' in inwcs.wcs.ctype):
+            raise ValueError("Must provide projection type or input WCS with CTYPE.")
+        projection = inwcs.wcs.ctype[0].replace('-SIP','').split('-')[-1]
+
+    # template wcs
+    wcs = celestial_frame_to_wcs(frame=coords.frame, projection=projection)
+
+    # Determine CRVAL from input, and
+    close = lambda l, p: p[np.where(np.abs(l) == min(np.abs(l)))[0][0]]
+    if str(proj_point) == 'center': #use center of input points
+        wcs.wcs.crval = ((max(lon)+min(lon))/2.,(max(lat)+min(lat))/2.)
+        wcs.wcs.crpix = ((max(xp)+min(xp))/2.,(max(yp)+min(yp))/2.) #initial guess
+    elif (proj_point is None) and (inwcs is None):
+        raise ValueError("Must give proj_point as argument or as CRVAL in input wcs.")
+    elif proj_point is not None: # convert units + rough initial guess for crpix for fit
+        lon_u, lat_u = u.Unit(coords.data.lon.unit), u.Unit(coords.data.lat.unit)
+        wcs.wcs.crval = (proj_point[0]*lon_u.to(u.deg), proj_point[1]*lat_u.to(u.deg))
+        wcs.wcs.crpix = (close(lon-wcs.wcs.crval[0],xp),close(lon-wcs.wcs.crval[1],yp))
+
+    if inwcs is not None:
+        if inwcs.wcs.radesys == '': # no frame specified, use Skycoord's frame (warn???)
+            inwcs.wcs.radesys = coords.frame.name
+        wcs.wcs.radesys = inwcs.wcs.radesys
+        coords.transform_to(inwcs.wcs.radesys.lower()) #work in inwcs units
+        if proj_point is None: # crval, crpix from wcs. crpix used for initial guess.
+            wcs.wcs.crval, wcs.wcs.crpix = inwcs.wcs.crval, inwcs.wcs.crpix
+            if wcs.wcs.crpix[0] == wcs.wcs.crpix[1] == 0: # assume 0 wasn't intentional
+                wcs.wcs.crpix = (close(lon-wcs.wcs.crval[0],xp), \
+                                 close(lon-wcs.wcs.crval[1],yp))
+        if inwcs.wcs.has_pc():
+            wcs.wcs.pc = inwcs.wcs.pc # work internally with pc
+        if inwcs.wcs.has_cd():
+            wcs.wcs.pc = inwcs.wcs.cd
+        # create dictionaries of both wcs (input and kwarg), merge them favoring kwargs
+        wcs_dict = dict(wcs.to_header(relax=False))
+        in_wcs_dict = dict(inwcs.to_header(relax=False))
+        wcs = WCS({**in_wcs_dict,**wcs_dict})
+        wcs.pixel_shape = inwcs.pixel_shape
+    else:
+        wcs.pixel_shape = (max(xp) - min(xp), max(yp) - min(yp))
+
+    # fit
+    p0 = np.concatenate((wcs.wcs.pc.flatten(),wcs.wcs.crpix.flatten()))
+    fit = least_squares(_linear_transformation_fit, p0, args=(lon, lat, xp, yp, wcs))
+
+    # put fit values in wcs
+    wcs.wcs.crpix = np.array(fit.x[4:6])
+    wcs.wcs.pc = np.array(fit.x[0:4].reshape((2,2)))
+
+    if mode == "all":
+        wcs.wcs.ctype = [x + '-SIP' for x in wcs.wcs.ctype]
+        if (order is None) or (type(order) == float):
+            raise ValueError("Must provide integer or tuple (a_order, b_order) for SIP.")
+        if type(order) == int:
+            a_order = b_order = order
+        elif (len(order) == 2):
+            if (type(order[0]) != int) or (type(order[1]) != int):
+                raise ValueError("Must provide integer or tuple (a_order, b_order) for SIP.")
+            a_order, b_order = order
+        else:
+            raise ValueError("Must provide integer or tuple (a_order, b_order) for SIP.")
+        a_coef_names = ['{0}_{1}'.format(i,j) for i in range(a_order+1) \
+                        for j in range(a_order+1) if (i+j) < (a_order+1) and (i+j) > 1]
+        b_coef_names = ['{0}_{1}'.format(i,j) for i in range(b_order+1) \
+                        for j in range(b_order+1) if (i+j) < (b_order + 1) and (i+j) > 1]
+        # fit
+        p0 = np.concatenate((np.array(wcs.wcs.crpix), wcs.wcs.pc.flatten(), \
+                             np.zeros(len(a_coef_names)+len(b_coef_names))))
+        fit = least_squares(_sip_fit, p0, args=(lon, lat, xp, yp, wcs, a_order, b_order,
+                                                a_coef_names, b_coef_names))
+        # put fit values in wcs
+        wcs.wcs.pc = fit.x[2:6].reshape((2,2))
+        wcs.wcs.crpix = fit.x[0:2]
+        coef_fit = (list(fit.x[6:6+len(a_coef_names)]),list(fit.x[6+len(b_coef_names):]))
+        a_vals, b_vals = np.zeros((a_order+1,a_order+1)), np.zeros((b_order+1,b_order+1))
+        for coef_name in a_coef_names:
+            a_vals[int(coef_name[0])][int(coef_name[2])] = coef_fit[0].pop(0)
+        for coef_name in b_coef_names:
+            b_vals[int(coef_name[0])][int(coef_name[2])] = coef_fit[1].pop(0)
+        wcs.sip = Sip(a_vals, b_vals, a_vals * -1., b_vals * -1., wcs.wcs.crpix)
+
+    if (inwcs is not None): # maintain cd matrix if it was in input wcs
+        if inwcs.wcs.has_cd():
+            wcs.wcs.cd = wcs.wcs.pc
+            wcs.wcs.__delattr__('pc')
+
+    return wcs

From 79452c0996989e210529dae6ff87a18e5f9bd54a Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Thu, 13 Jun 2019 08:54:39 -0400
Subject: [PATCH 02/15] adding goodness of fit print statement

---
 astrocut/cube_cut.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/astrocut/cube_cut.py b/astrocut/cube_cut.py
index 981c116d..918e58aa 100644
--- a/astrocut/cube_cut.py
+++ b/astrocut/cube_cut.py
@@ -844,6 +844,10 @@ def cube_cut(self, cube_file, coordinates, cutout_size,
         #cutout_wcs_dict = self._get_cutout_wcs()
         cutout_wcs_full = self._get_full_cutout_wcs(cube[2].header)
         max_dist,sigma = self._fit_cutout_wcs(cutout_wcs_full, img_cutout.shape[1:])
+        if verbose:
+            print("Maximum distance between approximate and true location: {}".format(max_dist))
+            print("Error in approximate WCS (sigma): {}".format(sigma))
+        
         # TODO: put the max_dist, sigma somewhere
         cutout_wcs_dict = self._get_cutout_wcs_dict()
     

From 78f2c2a4bd717a22bee3bdffd01a0f1a2e55e9fc Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Thu, 13 Jun 2019 14:55:40 -0400
Subject: [PATCH 03/15] adding fallback option

---
 astrocut/cube_cut.py | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/astrocut/cube_cut.py b/astrocut/cube_cut.py
index 918e58aa..ebf6416a 100644
--- a/astrocut/cube_cut.py
+++ b/astrocut/cube_cut.py
@@ -227,18 +227,22 @@ def _fit_cutout_wcs(self, cutout_wcs, cutout_shape):
         self.cutout_wcs = linear_wcs
 
         # Checking the fit
-        #world_pix_nosip = SkyCoord(cutout_wcs.wcs_pix2world(pix_inds,1), unit='deg')
+        world_pix_nosip = SkyCoord(cutout_wcs.wcs_pix2world(pix_inds,1), unit='deg')
         world_pix_new = SkyCoord(linear_wcs.all_pix2world(pix_inds,1), unit='deg')
     
         dists = world_pix.separation(world_pix_new)
-        #dists_nosip = world_pix.separation(world_pix_nosip)
+        dists_nosip = world_pix.separation(world_pix_nosip)
     
         sigma = np.sqrt(sum(dists.value**2))
-        #nosip_sigma = np.sqrt(sum(dists_nosip.value**2))
+        nosip_sigma = np.sqrt(sum(dists_nosip.value**2))
 
         # TODO: If just wacking off the SIP coefficients is better than the fit, just do that
-        #if nosip_sigma < sigma:
-        #    pass
+        if nosip_sigma < sigma:
+            if verbose:
+                print("Falling back to original WCS.")
+            self.cutout_wcs = cutout_wcs
+            dists = dists_nosip
+            sigma = nosip_sigma
 
         return (dists.max(), sigma)
     

From a958f0ca794405645cfa86fa561eca46d443dc33 Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Thu, 13 Jun 2019 15:01:59 -0400
Subject: [PATCH 04/15] fixing bug

---
 astrocut/cube_cut.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/astrocut/cube_cut.py b/astrocut/cube_cut.py
index ebf6416a..ca97e5a7 100644
--- a/astrocut/cube_cut.py
+++ b/astrocut/cube_cut.py
@@ -197,7 +197,7 @@ def _get_cutout_limits(self, cutout_size):
         self.cutout_lims = lims
 
 
-    def _fit_cutout_wcs(self, cutout_wcs, cutout_shape):
+    def _fit_cutout_wcs(self, cutout_wcs, cutout_shape, verbose=False):
         """
         Given a full (including SIP coefficients) wcs for the cutout, 
         calculate the best fit linear wcs, and a measure of the goodness-of-fit.
@@ -847,7 +847,7 @@ def cube_cut(self, cube_file, coordinates, cutout_size,
         # Get cutout wcs info
         #cutout_wcs_dict = self._get_cutout_wcs()
         cutout_wcs_full = self._get_full_cutout_wcs(cube[2].header)
-        max_dist,sigma = self._fit_cutout_wcs(cutout_wcs_full, img_cutout.shape[1:])
+        max_dist,sigma = self._fit_cutout_wcs(cutout_wcs_full, img_cutout.shape[1:], verbose=verbose)
         if verbose:
             print("Maximum distance between approximate and true location: {}".format(max_dist))
             print("Error in approximate WCS (sigma): {}".format(sigma))

From bba4d7d28b8c3987afc1f9bff0680692356ddf2e Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Fri, 14 Jun 2019 13:45:25 -0400
Subject: [PATCH 05/15] fixing wcs fit bugs

---
 astrocut/cube_cut.py          | 21 +++++++++++----------
 astrocut/utils/wcs_fitting.py | 19 +++++++++++++++++++
 2 files changed, 30 insertions(+), 10 deletions(-)

diff --git a/astrocut/cube_cut.py b/astrocut/cube_cut.py
index ca97e5a7..a0da3144 100644
--- a/astrocut/cube_cut.py
+++ b/astrocut/cube_cut.py
@@ -222,27 +222,28 @@ def _fit_cutout_wcs(self, cutout_wcs, cutout_shape, verbose=False):
         world_pix = SkyCoord(cutout_wcs.all_pix2world(pix_inds, 1), unit='deg')
 
         # Getting the fit WCS
-        linear_wcs = fit_wcs_from_points(pix_inds[:,0], pix_inds[:,1], world_pix, mode='wcs', proj_point="center")
+        linear_wcs = fit_wcs_from_points(pix_inds[:,0], pix_inds[:,1], world_pix, mode='wcs',
+                                         proj_point=[self.center_coord.data.lon.value,self.center_coord.data.lat.value])
 
         self.cutout_wcs = linear_wcs
 
         # Checking the fit
-        world_pix_nosip = SkyCoord(cutout_wcs.wcs_pix2world(pix_inds,1), unit='deg')
+        #world_pix_nosip = SkyCoord(cutout_wcs.wcs_pix2world(pix_inds,1), unit='deg')
         world_pix_new = SkyCoord(linear_wcs.all_pix2world(pix_inds,1), unit='deg')
     
         dists = world_pix.separation(world_pix_new)
-        dists_nosip = world_pix.separation(world_pix_nosip)
+        #dists_nosip = world_pix.separation(world_pix_nosip)
     
         sigma = np.sqrt(sum(dists.value**2))
-        nosip_sigma = np.sqrt(sum(dists_nosip.value**2))
+        #nosip_sigma = np.sqrt(sum(dists_nosip.value**2))
 
         # TODO: If just wacking off the SIP coefficients is better than the fit, just do that
-        if nosip_sigma < sigma:
-            if verbose:
-                print("Falling back to original WCS.")
-            self.cutout_wcs = cutout_wcs
-            dists = dists_nosip
-            sigma = nosip_sigma
+        #if nosip_sigma < sigma:
+        #    if verbose:
+        #        print("Falling back to original WCS.")
+        #    self.cutout_wcs = cutout_wcs
+        #    dists = dists_nosip
+        #    sigma = nosip_sigma
 
         return (dists.max(), sigma)
     
diff --git a/astrocut/utils/wcs_fitting.py b/astrocut/utils/wcs_fitting.py
index bbb3a70c..2a0985f8 100644
--- a/astrocut/utils/wcs_fitting.py
+++ b/astrocut/utils/wcs_fitting.py
@@ -1,3 +1,19 @@
+#################################################################################
+#
+# Licensed under a 3-clause BSD style license
+#           - see https://github.com/astropy/astropy/blob/master/LICENSE.rst
+#
+# wcs fitting functionality
+#          by Clare Shanahan (shannnnnyyy @github)
+#
+# Will be added to Astropy (PR: https://github.com/astropy/astropy/pull/7884)
+#
+# Astropy version is used preferenetially, this is supplied prior to the
+# addition of this code to Astropy, and for users using older versions of Astropy
+#
+#################################################################################
+
+
 import numpy as np
 
 from astropy import units as u
@@ -12,7 +28,10 @@ def _linear_transformation_fit(params, lon, lat, x, y, w_obj):
     w_obj.wcs.pc = ((pc[0],pc[1]),(pc[2],pc[3]))
     w_obj.wcs.crpix = crpix
     lon2, lat2 = w_obj.wcs_pix2world(x,y,1)
+    
     resids = np.concatenate((lon-lon2, lat-lat2))
+    resids[resids > 180] = 360 - resids[resids > 180]
+    resids[resids < -180] = 360	+ resids[resids < -180]
 
     return resids
 

From b5f0c8ec07e1b337d2cca9d3ee75cca81fb47cbe Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Fri, 14 Jun 2019 15:50:32 -0400
Subject: [PATCH 06/15] removing unused code and cleaning up a bit

---
 astrocut/cube_cut.py | 155 +++++--------------------------------------
 1 file changed, 17 insertions(+), 138 deletions(-)

diff --git a/astrocut/cube_cut.py b/astrocut/cube_cut.py
index a0da3144..bad50442 100644
--- a/astrocut/cube_cut.py
+++ b/astrocut/cube_cut.py
@@ -12,8 +12,8 @@
 from itertools import product
 
 try:
-    from astropy.utils import fit_wcs_from_points
-except: # astropy version does have the function
+    from astropy.wcs.utils import fit_wcs_from_points
+except: # astropy version does not have the function
     from .utils.wcs_fitting import fit_wcs_from_points
 
 from time import time
@@ -41,6 +41,9 @@ def __init__(self):
 
         self.cube_wcs = None  # WCS information from the image cube
         self.cutout_wcs = None # WCS information (linear) for the cutout
+        self.cutout_wcs_fit = {'WCS_MSEP': [None, "[deg] Max offset between cutout WCS and FFI WCS"],
+                               'WCS_SIG': [None, "[deg] Error measurement of cutout WCS fit"]}
+        
         self.cutout_lims = np.zeros((2, 2), dtype=int)  # Cutout pixel limits, [[ymin,ymax],[xmin,xmax]]
         self.center_coord = None  # Central skycoord
         
@@ -197,12 +200,13 @@ def _get_cutout_limits(self, cutout_size):
         self.cutout_lims = lims
 
 
-    def _fit_cutout_wcs(self, cutout_wcs, cutout_shape, verbose=False):
+    def _fit_cutout_wcs(self, cutout_wcs, cutout_shape):
         """
         Given a full (including SIP coefficients) wcs for the cutout, 
         calculate the best fit linear wcs, and a measure of the goodness-of-fit.
         
-        The news WCS is stored in ``self.cutout_wcs``.
+        The new WCS is stored in ``self.cutout_wcs``.
+        Goodness-of-fit measures are returned and stored in ``self.cutout_wcs_fit``.
 
         Parameters
         ----------
@@ -228,22 +232,13 @@ def _fit_cutout_wcs(self, cutout_wcs, cutout_shape, verbose=False):
         self.cutout_wcs = linear_wcs
 
         # Checking the fit
-        #world_pix_nosip = SkyCoord(cutout_wcs.wcs_pix2world(pix_inds,1), unit='deg')
         world_pix_new = SkyCoord(linear_wcs.all_pix2world(pix_inds,1), unit='deg')
     
         dists = world_pix.separation(world_pix_new)
-        #dists_nosip = world_pix.separation(world_pix_nosip)
-    
         sigma = np.sqrt(sum(dists.value**2))
-        #nosip_sigma = np.sqrt(sum(dists_nosip.value**2))
 
-        # TODO: If just wacking off the SIP coefficients is better than the fit, just do that
-        #if nosip_sigma < sigma:
-        #    if verbose:
-        #        print("Falling back to original WCS.")
-        #    self.cutout_wcs = cutout_wcs
-        #    dists = dists_nosip
-        #    sigma = nosip_sigma
+        self.cutout_wcs_fit['WCS_MSEP'][0] = dists.max().to('deg').value
+        self.cutout_wcs_fit['WCS_SIG'][0] = sigma
 
         return (dists.max(), sigma)
     
@@ -351,80 +346,6 @@ def _get_cutout_wcs_dict(self):
 
         return cutout_wcs_dict
 
-    
-    def _get_cutout_wcs(self):
-        """
-        Transform the cube WCS object into the cutout column WCS keywords.
-        Adds the physical keywords for transformation back from cutout to location on FFI.
-        This is a very TESS specific function.
-        
-        Returns
-        -------
-        response: dict
-            Cutout wcs column header keywords as dictionary of 
-            ``{<kwd format string>: [value, desc]} pairs.``
-        """
-        cube_wcs_header = self.cube_wcs.to_header(relax=True)
-
-        cutout_wcs_dict = dict()
-
-        
-        ## Cutout array keywords ##
-
-        cutout_wcs_dict["WCAX{}"] = [cube_wcs_header['WCSAXES'], "number of WCS axes"]
-        # TODO: check for 2? this must be two
-
-        cutout_wcs_dict["1CTYP{}"] = [cube_wcs_header["CTYPE1"], "right ascension coordinate type"]
-        cutout_wcs_dict["2CTYP{}"] = [cube_wcs_header["CTYPE2"], "declination coordinate type"]
-        
-        cutout_wcs_dict["1CRPX{}"] = [cube_wcs_header["CRPIX1"] - self.cutout_lims[0, 0],
-                                      "[pixel] reference pixel along image axis 1"]
-        cutout_wcs_dict["2CRPX{}"] = [cube_wcs_header["CRPIX2"] - self.cutout_lims[1, 0],
-                                      "[pixel] reference pixel along image axis 2"]
-    
-        cutout_wcs_dict["1CRVL{}"] = [cube_wcs_header["CRVAL1"], "[deg] right ascension at reference pixel"]
-        cutout_wcs_dict["2CRVL{}"] = [cube_wcs_header["CRVAL2"], "[deg] declination at reference pixel"]
-
-        cunits = self.cube_wcs.wcs.cunit
-        cutout_wcs_dict["1CUNI{}"] = [str(cunits[0]), "physical unit in column dimension"]
-        cutout_wcs_dict["2CUNI{}"] = [str(cunits[1]), "physical unit in row dimension"]
-
-        px_scales = wcs.utils.proj_plane_pixel_scales(self.cube_wcs)
-        cutout_wcs_dict["1CDLT{}"] = [px_scales[0], "[deg] pixel scale in RA dimension"]
-        cutout_wcs_dict["2CDLT{}"] = [px_scales[1], "[deg] pixel scale in DEC dimension"]
-
-        # TODO: THIS IS FILLER, HAVE TO FIGURE OUT HOW TO DO THE TRANSFORMATION FOR REAL
-        cutout_wcs_dict["11PC{}"] = [1, "linear transformation matrix element - unfilled"]
-        cutout_wcs_dict["12PC{}"] = [1, "linear transformation matrix element - unfilled"]
-        cutout_wcs_dict["21PC{}"] = [1, "linear transformation matrix element - unfilled"]
-        cutout_wcs_dict["22PC{}"] = [1, "linear transformation matrix element - unfilled"]
-
-        
-        ## Physical keywords ##
-    
-        cutout_wcs_dict["WCSN{}P"] = ["PHYSICAL", "table column WCS name"]
-        cutout_wcs_dict["WCAX{}P"] = [2, "table column physical WCS dimensions"]
-    
-        cutout_wcs_dict["1CTY{}P"] = ["RAWX", "table column physical WCS axis 1 type, CCD col"]
-        cutout_wcs_dict["2CTY{}P"] = ["RAWY", "table column physical WCS axis 2 type, CCD row"]
-    
-        cutout_wcs_dict["1CUN{}P"] = ["PIXEL", "table column physical WCS axis 1 unit"]
-        cutout_wcs_dict["2CUN{}P"] = ["PIXEL", "table column physical WCS axis 2 unit"]
-    
-        cutout_wcs_dict["1CRV{}P"] = [self.cutout_lims[0, 0] + 1,
-                                      "table column physical WCS ax 1 ref value"]
-        cutout_wcs_dict["2CRV{}P"] = [self.cutout_lims[1, 0] + 1,
-                                      "table column physical WCS ax 2 ref value"]
-
-        # TODO: can we calculate these? or are they fixed?
-        cutout_wcs_dict["1CDL{}P"] = [1.0, "table column physical WCS a1 step"]    
-        cutout_wcs_dict["2CDL{}P"] = [1.0, "table column physical WCS a2 step"]
-    
-        cutout_wcs_dict["1CRP{}P"] = [1, "table column physical WCS a1 reference"]
-        cutout_wcs_dict["2CRP{}P"] = [1, "table column physical WCS a2 reference"]
-
-        return cutout_wcs_dict
-    
 
     def _get_cutout(self, transposed_cube, verbose=True):
         """
@@ -577,51 +498,6 @@ def _add_column_wcs(self, table_header, wcs_dict):
                 for wcs_key, (val, com) in wcs_dict.items():
                     table_header.insert(kwd, (wcs_key.format(int(kwd[-1])-1), val, com))
 
-
-    def _add_aperture_wcs(self, aperture_header, cube_table_header):
-        """
-        Adds the full cube WCS information with adjustments for the cutout 
-        to the aperture extension header in place.
-
-        Parameters
-        ----------
-        aperture_header : `~astropy.io.fits.Header`
-            The aperture extension header.  It will be modified in place.
-        cube_table_header : `~astropy.io.fits.Header`
-            The header for the table extension header from the cube file.
-        """
-
-        cube_wcs_header = self.cube_wcs.to_header(relax=True)
-
-        # Adding the wcs keywords
-        for kwd, val, cmt in cube_wcs_header.cards: 
-            aperture_header.set(kwd, val, cube_table_header.get(kwd, cmt))
-            # using table comment rather than the default ones if available
-
-        # Adjusting the CRPIX values
-        aperture_header["CRPIX1"] -= self.cutout_lims[0, 0]
-        aperture_header["CRPIX2"] -= self.cutout_lims[1, 0]
-
-        # Adding the physical wcs keywords
-        aperture_header.set("WCSNAMEP", "PHYSICAL", "name of world coordinate system alternate P")
-        aperture_header.set("WCSAXESP", 2, "number of WCS physical axes")
-    
-        aperture_header.set("CTYPE1P", "RAWX", "physical WCS axis 1 type CCD col")
-        aperture_header.set("CUNIT1P", "PIXEL", "physical WCS axis 1 unit")
-        aperture_header.set("CRPIX1P", 1, "reference CCD column")
-        aperture_header.set("CRVAL1P", self.cutout_lims[0, 0] + 1, "value at reference CCD column")
-        aperture_header.set("CDELT1P", 1.0, "physical WCS axis 1 step")
-                
-        aperture_header.set("CTYPE2P", "RAWY", "physical WCS axis 2 type CCD col")
-        aperture_header.set("CUNIT2P", "PIXEL", "physical WCS axis 2 unit")
-        aperture_header.set("CRPIX2P", 1, "reference CCD row")
-        aperture_header.set("CRVAL2P", self.cutout_lims[1, 0] + 1, "value at reference CCD row")
-        aperture_header.set("CDELT2P", 1.0, "physical WCS axis 2 step")
-
-        # Adding extra aperture keywords (TESS specific)
-        aperture_header.set("NPIXMISS", None, "Number of op. aperture pixels not collected")
-        aperture_header.set("NPIXSAP", None, "Number of pixels in optimal aperture")
-
         
     def _add_img_kwds(self, table_header):
         """
@@ -750,10 +626,15 @@ def _build_tpf(self, cube_fits, img_cube, uncert_cube, cutout_wcs_dict, aperture
         aperture_hdu.header['EXTNAME'] = 'APERTURE'
         for kwd, val, cmt in self.cutout_wcs.to_header().cards: 
             aperture_hdu.header.set(kwd, val, cmt)
-        #self._add_aperture_wcs(aperture_hdu.header, cube_fits[2].header)
+
         # Adding extra aperture keywords (TESS specific)
         aperture_hdu.header.set("NPIXMISS", None, "Number of op. aperture pixels not collected")
         aperture_hdu.header.set("NPIXSAP", None, "Number of pixels in optimal aperture")
+
+        # Adding goodness-of-fit keywords
+        for key in self.cutout_wcs_fit:
+            aperture_hdu.header[key] = tuple(self.cutout_wcs_fit[key])
+        
         aperture_hdu.header['INHERIT'] = True
     
         cutout_hdu_list = fits.HDUList([primary_hdu, table_hdu, aperture_hdu])
@@ -846,14 +727,12 @@ def cube_cut(self, cube_file, coordinates, cutout_size,
         img_cutout, uncert_cutout, aperture = self._get_cutout(cube[1].data, verbose=verbose)
 
         # Get cutout wcs info
-        #cutout_wcs_dict = self._get_cutout_wcs()
         cutout_wcs_full = self._get_full_cutout_wcs(cube[2].header)
-        max_dist,sigma = self._fit_cutout_wcs(cutout_wcs_full, img_cutout.shape[1:], verbose=verbose)
+        max_dist, sigma = self._fit_cutout_wcs(cutout_wcs_full, img_cutout.shape[1:])
         if verbose:
             print("Maximum distance between approximate and true location: {}".format(max_dist))
             print("Error in approximate WCS (sigma): {}".format(sigma))
         
-        # TODO: put the max_dist, sigma somewhere
         cutout_wcs_dict = self._get_cutout_wcs_dict()
     
         # Build the TPF

From ecbc79a6d51b86808d8e4c35225452b906463e26 Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Mon, 17 Jun 2019 11:13:35 -0400
Subject: [PATCH 07/15] adding unit tests for new functionality

---
 astrocut/tests/test_cube_cut.py | 30 ++++++++++++++++++++----------
 1 file changed, 20 insertions(+), 10 deletions(-)

diff --git a/astrocut/tests/test_cube_cut.py b/astrocut/tests/test_cube_cut.py
index 06c6b87a..223ffb88 100644
--- a/astrocut/tests/test_cube_cut.py
+++ b/astrocut/tests/test_cube_cut.py
@@ -204,21 +204,31 @@ def test_cutout_extras(tmpdir):
     assert (xmax-xmin) == 1
     assert (ymax-ymin) == 1
 
-    #########################
-    # Test  _get_cutout_wcs #
-    #########################
+    #############################
+    # Test _get_full_cutout_wcs #
+    #############################
     cutout_size = [5, 3]
     out_file = myfactory.cube_cut(cube_file, coord, cutout_size, verbose=False)
 
-    cutout_wcs_dict = myfactory._get_cutout_wcs()
+    cutout_wcs_full = myfactory._get_full_cutout_wcs(fits.getheader(cube_file, 2))
+    assert (cutout_wcs_full.wcs.crpix == [1045 - myfactory.cutout_lims[0, 0],
+                                          1001 - myfactory.cutout_lims[1, 0]]).all()    
 
-    assert cutout_wcs_dict["1CTYP{}"][0] == 'RA---TAN-SIP'
-    assert cutout_wcs_dict["1CRPX{}"][0] == 1043
-    assert round(cutout_wcs_dict["1CRVL{}"][0], 4) == 250.3497
+    ########################
+    # Test _fit_cutout_wcs #
+    ########################
+    max_dist, sigma = myfactory._fit_cutout_wcs(cutout_wcs_full, (3,5))
+    assert round(max_dist.deg,9) == 6.97e-07
+    assert round(sigma,7) == 1.4e-06
 
-    ###########################
-    # Test  target pixel file #
-    ###########################
+    cry, crx =  myfactory.cutout_wcs.wcs.crpix
+    assert round(cry) == 4
+    assert round(crx) == 2
+    
+
+    ##########################
+    # Test target pixel file #
+    ##########################
 
     # Testing the cutout content is in test_cube_cutout
     # this tests that the format of the tpf is what it should be

From c70564456154a9d7a9ddd294a74073ec219e04e2 Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Tue, 18 Jun 2019 14:26:21 -0400
Subject: [PATCH 08/15] documenting the new wcs functionality

---
 astrocut/cube_cut.py    | 93 ++++++++++++++++++++++++-----------------
 docs/astrocut/index.rst |  9 ++++
 2 files changed, 63 insertions(+), 39 deletions(-)

diff --git a/astrocut/cube_cut.py b/astrocut/cube_cut.py
index bad50442..05eb49c5 100644
--- a/astrocut/cube_cut.py
+++ b/astrocut/cube_cut.py
@@ -145,6 +145,9 @@ def _parse_table_info(self, table_header, table_data, verbose=False):
         # Filling the img_kwds dictionary while we are here
         for kwd in self.img_kwds:
             self.img_kwds[kwd][0] = wcs_header.get(kwd)
+        # Adding the info about which FFI we got the 
+        self.img_kwds["WCS_FFI"] = [table_data[data_ind]["FFI_FILE"],
+                                    "FFI filname used in cutout WCS calculation."]
 
             
     def _get_cutout_limits(self, cutout_size):
@@ -200,6 +203,53 @@ def _get_cutout_limits(self, cutout_size):
         self.cutout_lims = lims
 
 
+    def _get_full_cutout_wcs(self, cube_table_header):
+        """
+        Starting with the full FFI WCS and adjusting it for the cutout WCS.
+        Adjusts CRPIX values and adds physical WCS keywords.
+
+        Parameters
+        ----------
+        cube_table_header :  `~astropy.io.fits.Header`
+           The FFI cube header for the data table extension. This allows the cutout WCS information
+           to more closely match the mission TPF format.
+
+        Resturns
+        --------
+        response :  `~astropy.wcs.WCS`
+            The cutout WCS object including SIP distortions.
+        """
+        
+        wcs_header = self.cube_wcs.to_header(relax=True)
+
+        # Using table comment rather than the default ones if available
+        for kwd in wcs_header:
+            if cube_table_header.get(kwd, None):
+                wcs_header.comments[kwd] = cube_table_header[kwd]
+
+        # Adjusting the CRPIX values
+        wcs_header["CRPIX1"] -= self.cutout_lims[0, 0]
+        wcs_header["CRPIX2"] -= self.cutout_lims[1, 0]
+
+        # Adding the physical wcs keywords
+        wcs_header.set("WCSNAMEP", "PHYSICAL", "name of world coordinate system alternate P")
+        wcs_header.set("WCSAXESP", 2, "number of WCS physical axes")
+    
+        wcs_header.set("CTYPE1P", "RAWX", "physical WCS axis 1 type CCD col")
+        wcs_header.set("CUNIT1P", "PIXEL", "physical WCS axis 1 unit")
+        wcs_header.set("CRPIX1P", 1, "reference CCD column")
+        wcs_header.set("CRVAL1P", self.cutout_lims[0, 0] + 1, "value at reference CCD column")
+        wcs_header.set("CDELT1P", 1.0, "physical WCS axis 1 step")
+                
+        wcs_header.set("CTYPE2P", "RAWY", "physical WCS axis 2 type CCD col")
+        wcs_header.set("CUNIT2P", "PIXEL", "physical WCS axis 2 unit")
+        wcs_header.set("CRPIX2P", 1, "reference CCD row")
+        wcs_header.set("CRVAL2P", self.cutout_lims[1, 0] + 1, "value at reference CCD row")
+        wcs_header.set("CDELT2P", 1.0, "physical WCS axis 2 step")
+
+        return wcs.WCS(wcs_header)
+
+    
     def _fit_cutout_wcs(self, cutout_wcs, cutout_shape):
         """
         Given a full (including SIP coefficients) wcs for the cutout, 
@@ -210,7 +260,7 @@ def _fit_cutout_wcs(self, cutout_wcs, cutout_shape):
 
         Parameters
         ----------
-        cutout_wcs : 
+        cutout_wcs :  `~astropy.wcs.WCS`
             The full (including SIP coefficients) cutout WCS object 
         cutout_shape : tuple
             The shape of the cutout in the form (width, height).
@@ -234,50 +284,15 @@ def _fit_cutout_wcs(self, cutout_wcs, cutout_shape):
         # Checking the fit
         world_pix_new = SkyCoord(linear_wcs.all_pix2world(pix_inds,1), unit='deg')
     
-        dists = world_pix.separation(world_pix_new)
+        dists = world_pix.separation(world_pix_new).to('deg')
         sigma = np.sqrt(sum(dists.value**2))
 
-        self.cutout_wcs_fit['WCS_MSEP'][0] = dists.max().to('deg').value
+        self.cutout_wcs_fit['WCS_MSEP'][0] = dists.max().value
         self.cutout_wcs_fit['WCS_SIG'][0] = sigma
 
         return (dists.max(), sigma)
     
-
-    def _get_full_cutout_wcs(self, cube_table_header):
-        """
-        TODO: Document
-        """
-        
-        wcs_header = self.cube_wcs.to_header(relax=True)
-
-        # Using table comment rather than the default ones if available
-        for kwd in wcs_header:
-            if cube_table_header.get(kwd, None):
-                wcs_header.comments[kwd] = cube_table_header[kwd]
-
-        # Adjusting the CRPIX values
-        wcs_header["CRPIX1"] -= self.cutout_lims[0, 0]
-        wcs_header["CRPIX2"] -= self.cutout_lims[1, 0]
-
-        # Adding the physical wcs keywords
-        wcs_header.set("WCSNAMEP", "PHYSICAL", "name of world coordinate system alternate P")
-        wcs_header.set("WCSAXESP", 2, "number of WCS physical axes")
-    
-        wcs_header.set("CTYPE1P", "RAWX", "physical WCS axis 1 type CCD col")
-        wcs_header.set("CUNIT1P", "PIXEL", "physical WCS axis 1 unit")
-        wcs_header.set("CRPIX1P", 1, "reference CCD column")
-        wcs_header.set("CRVAL1P", self.cutout_lims[0, 0] + 1, "value at reference CCD column")
-        wcs_header.set("CDELT1P", 1.0, "physical WCS axis 1 step")
-                
-        wcs_header.set("CTYPE2P", "RAWY", "physical WCS axis 2 type CCD col")
-        wcs_header.set("CUNIT2P", "PIXEL", "physical WCS axis 2 unit")
-        wcs_header.set("CRPIX2P", 1, "reference CCD row")
-        wcs_header.set("CRVAL2P", self.cutout_lims[1, 0] + 1, "value at reference CCD row")
-        wcs_header.set("CDELT2P", 1.0, "physical WCS axis 2 step")
-
-        return wcs.WCS(wcs_header)
-
-    
+  
     def _get_cutout_wcs_dict(self):
         """
         Transform the cutout WCS object into the cutout column WCS keywords.
diff --git a/docs/astrocut/index.rst b/docs/astrocut/index.rst
index bee57681..c9981123 100644
--- a/docs/astrocut/index.rst
+++ b/docs/astrocut/index.rst
@@ -119,6 +119,15 @@ be saved, if unspecified it defaults to the current directory.
                 2  APERTURE      1 ImageHDU        45   (5, 5)   float64  
 
 
+A note about the cutout WCS object
+""""""""""""""""""""""""""""""""""
+
+TESS FFIs are large and therefore are described by WCS objects that have many non-linear terms. Astrocut creates a new simpler (linear) WCS object from the matched set of cutout pixel coordinates and sky coordinates (from the FFI WCS). This linear WCS object will generally work very well, however at larger cutout sizes (100-200 pixels per side and above) the linear WCS fit will start to be noticeably incorrect at the edges of the cutout.  Three header keywords have been added to the PIXELS extension to give additional information about the cutout WCS:
+
+* **WCS_FFI:** The name of the FFI file used to build the original WCS from which the cutout and cutout WCS were calculated.
+* **WCS_MSEP:** The maximum separation in degrees between the cutout's linear WCS and the FFI's full WCS.
+* **WCS_SIG:** The error in the cutout's linear WCS, calculated as sqrt((dist(Po_ij, Pl_ij)^2) where dist(Po_ij, Pl_ij) is the angular distance in degrees between the sky position of of pixel i,j in the original full WCS and the new linear WCS.
+
 
 .. automodapi:: astrocut
     :skip: test

From a3f6b7291a3fe8fac034370dd8d7a89465684345 Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Tue, 18 Jun 2019 15:30:16 -0400
Subject: [PATCH 09/15] adding code to make output directories if needed

---
 astrocut/cube_cut.py  |  4 ++++
 astrocut/make_cube.py | 15 +++++++++------
 2 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/astrocut/cube_cut.py b/astrocut/cube_cut.py
index 05eb49c5..51c1c731 100644
--- a/astrocut/cube_cut.py
+++ b/astrocut/cube_cut.py
@@ -768,6 +768,10 @@ def cube_cut(self, cube_file, coordinates, cutout_size,
         
         if verbose:
             print("Target pixel file: {}".format(target_pixel_file))
+
+        # Make sure the output directory exists
+        if not os.path.exists(output_path):
+            os.makedirs(output_path)
         
         # Write the TPF
         tpf_object.writeto(target_pixel_file, overwrite=True, checksum=True)
diff --git a/astrocut/make_cube.py b/astrocut/make_cube.py
index a0cbe936..7f988c47 100644
--- a/astrocut/make_cube.py
+++ b/astrocut/make_cube.py
@@ -6,13 +6,13 @@
 """
 
 import numpy as np
+import os
 
 from astropy.io import fits
 from astropy.table import Table, Column
 
 from time import time
 from datetime import date
-from os import path
 from copy import deepcopy
 
 
@@ -122,7 +122,7 @@ def make_cube(self, file_list, cube_file="img-cube.fits", sector=None, verbose=T
         fle = file_list[good_header_ind]
         
         if verbose:
-            print("Using {} to initialize the image header table.".format(path.basename(fle)))
+            print("Using {} to initialize the image header table.".format(os.path.basename(fle)))
             
         ffi_data = fits.open(fle)
 
@@ -141,7 +141,7 @@ def make_cube(self, file_list, cube_file="img-cube.fits", sector=None, verbose=T
                     
             cols.append(Column(name=kwd, dtype=tpe, length=len(file_list), meta={"comment": cmt}))
                     
-        cols.append(Column(name="FFI_FILE", dtype="S" + str(len(path.basename(fle))), length=len(file_list)))
+        cols.append(Column(name="FFI_FILE", dtype="S" + str(len(os.path.basename(fle))), length=len(file_list)))
             
         img_info_table = Table(cols)
 
@@ -172,7 +172,7 @@ def make_cube(self, file_list, cube_file="img-cube.fits", sector=None, verbose=T
             img_cube[:, :, i, 1] = ffi_data[2].data
             for kwd in img_info_table.columns:
                 if kwd == "FFI_FILE":
-                    img_info_table[kwd][i] = path.basename(fle)
+                    img_info_table[kwd][i] = os.path.basename(fle)
                 else:
                     nulval = None
                     if img_info_table[kwd].dtype.name == "int32":
@@ -221,8 +221,11 @@ def make_cube(self, file_list, cube_file="img-cube.fits", sector=None, verbose=T
         if verbose:
             writeTime = time()
 
-        # TODO: check for path in cube and make sure directory exists
-        #       and create it if it doesn't
+        # Making sure the output directory exists
+        direc, _ = os.path.split(cube_file)
+        if direc and not os.path.exists(direc):
+            os.makedirs(direc)
+            
         hdu_list.writeto(cube_file, overwrite=True) 
 
         if verbose:

From cdb7c9e43621bb99bb172a333f2568b4a168694c Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Tue, 18 Jun 2019 15:34:29 -0400
Subject: [PATCH 10/15] updating required packages

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index 0ad402bd..f8217cc5 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -39,7 +39,7 @@ edit_on_github = False
 github_project = spacetelescope/astrocut
 # install_requires should be formatted as a comma-separated list, e.g.:
 # install_requires = astropy, scipy, matplotlib
-install_requires = astropy
+install_requires = astropy, scipy
 # version should be PEP440 compatible (https://www.python.org/dev/peps/pep-0440/)
 version = 0.4.dev
 # Note: you will also need to change this in your package's __init__.py

From facf9f5abd4cadee61f2a1d5499c2c6c232dcbc1 Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Wed, 19 Jun 2019 20:44:44 -0400
Subject: [PATCH 11/15] fixing overlong comment and updating wcs check per mit

---
 astrocut/cube_cut.py  | 6 +++---
 astrocut/make_cube.py | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/astrocut/cube_cut.py b/astrocut/cube_cut.py
index 51c1c731..57c048e6 100644
--- a/astrocut/cube_cut.py
+++ b/astrocut/cube_cut.py
@@ -105,8 +105,8 @@ def _parse_table_info(self, table_header, table_data, verbose=False):
         # Making sure we have a row with wcs info
         while table_row is None:
             table_row = table_data[data_ind]
-            ra_col = int([x for x in table_header.cards if x[1] == "CTYPE1"][0][0].replace("TTYPE", "")) - 1
-            if "RA" not in table_row[ra_col]:
+            ra_col = int([x for x in table_header.cards if x[1] == "WCSAXES"][0][0].replace("TTYPE", "")) - 1
+            if table_row[ra_col] == 2:
                 table_row is None
                 data_ind += 1
                 if data_ind == len(table_data):
@@ -147,7 +147,7 @@ def _parse_table_info(self, table_header, table_data, verbose=False):
             self.img_kwds[kwd][0] = wcs_header.get(kwd)
         # Adding the info about which FFI we got the 
         self.img_kwds["WCS_FFI"] = [table_data[data_ind]["FFI_FILE"],
-                                    "FFI filname used in cutout WCS calculation."]
+                                    "FFI used for cutout WCS"]
 
             
     def _get_cutout_limits(self, cutout_size):
diff --git a/astrocut/make_cube.py b/astrocut/make_cube.py
index 7f988c47..0df343f9 100644
--- a/astrocut/make_cube.py
+++ b/astrocut/make_cube.py
@@ -104,7 +104,7 @@ def make_cube(self, file_list, cube_file="img-cube.fits", sector=None, verbose=T
             start_times[i] = ffi_data[1].header.get("TSTART")  # TODO: optionally pass this in?
             
             if good_header_ind is None:  # Only check this if we don't already have it
-                if ffi_data[1].header.get("CTYPE1"):  # Checking for WCS info
+                if ffi_data[1].header.get("WCSAXES", 0) == 2:  # Checking for WCS info
                     good_header_ind = i
             
             ffi_data.close()

From 16fabc19c271a375027a778bb5c15cef6156d45b Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Wed, 19 Jun 2019 20:52:34 -0400
Subject: [PATCH 12/15] chengelog

---
 CHANGES.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGES.rst b/CHANGES.rst
index 828696e5..5a6f25d9 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -3,6 +3,7 @@
 
 - Adding more unit tests and coveralls setup [#11]
 - Adding workaround for FFIs with bad WCS info [#12]
+- Adding linear WCS approximation for cutouts [#14]
 
 
 0.3 (2019-05-03)

From 1048a01304c6fceff94078e993045ccf0c36c9e8 Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Thu, 20 Jun 2019 14:45:26 -0400
Subject: [PATCH 13/15] setting min astropy version

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index f8217cc5..1751f746 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -39,7 +39,7 @@ edit_on_github = False
 github_project = spacetelescope/astrocut
 # install_requires should be formatted as a comma-separated list, e.g.:
 # install_requires = astropy, scipy, matplotlib
-install_requires = astropy, scipy
+install_requires = astropy>=3.0, scipy
 # version should be PEP440 compatible (https://www.python.org/dev/peps/pep-0440/)
 version = 0.4.dev
 # Note: you will also need to change this in your package's __init__.py

From eafcd008e148cf587a88680cad2537d2f8295782 Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Thu, 20 Jun 2019 15:26:35 -0400
Subject: [PATCH 14/15] removing lts test from travis

---
 .travis.yml | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index b7accf91..11aae72d 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -89,8 +89,10 @@ matrix:
         - os: linux
           env: ASTROPY_VERSION=development
                EVENT_TYPE='pull_request push cron'
-        - os: linux
-          env: ASTROPY_VERSION=lts
+
+        # Min required version of astropy is 3.0, so skipping this test.
+        #- os: linux
+        #  env: ASTROPY_VERSION=lts
 
         # Try all python versions and Numpy versions. Since we can assume that
         # the Numpy developers have taken care of testing Numpy with different

From 911337fb975a9242e9400e5dfc0b70e8f40bfb3f Mon Sep 17 00:00:00 2001
From: Clara Brasseur <cbrasseur@stsci.edu>
Date: Fri, 21 Jun 2019 11:07:34 -0400
Subject: [PATCH 15/15] excluting utils file and pep8 fixes

---
 astrocut/cube_cut.py            | 12 ++++++------
 astrocut/tests/test_cube_cut.py |  8 ++++----
 astrocut/utils/wcs_fitting.py   |  8 ++++----
 3 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/astrocut/cube_cut.py b/astrocut/cube_cut.py
index 57c048e6..eac070a1 100644
--- a/astrocut/cube_cut.py
+++ b/astrocut/cube_cut.py
@@ -13,7 +13,7 @@
 
 try:
     from astropy.wcs.utils import fit_wcs_from_points
-except: # astropy version does not have the function
+except ImportError:  # astropy version does not have the function
     from .utils.wcs_fitting import fit_wcs_from_points
 
 from time import time
@@ -40,7 +40,7 @@ def __init__(self):
         """
 
         self.cube_wcs = None  # WCS information from the image cube
-        self.cutout_wcs = None # WCS information (linear) for the cutout
+        self.cutout_wcs = None  # WCS information (linear) for the cutout
         self.cutout_wcs_fit = {'WCS_MSEP': [None, "[deg] Max offset between cutout WCS and FFI WCS"],
                                'WCS_SIG': [None, "[deg] Error measurement of cutout WCS fit"]}
         
@@ -272,17 +272,17 @@ 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])))))
+        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, 1), unit='deg')
 
         # Getting the fit WCS
-        linear_wcs = fit_wcs_from_points(pix_inds[:,0], pix_inds[:,1], world_pix, mode='wcs',
-                                         proj_point=[self.center_coord.data.lon.value,self.center_coord.data.lat.value])
+        linear_wcs = fit_wcs_from_points(pix_inds[:, 0], pix_inds[:, 1], world_pix, mode='wcs',
+                                         proj_point=[self.center_coord.data.lon.value, self.center_coord.data.lat.value])
 
         self.cutout_wcs = linear_wcs
 
         # Checking the fit
-        world_pix_new = SkyCoord(linear_wcs.all_pix2world(pix_inds,1), unit='deg')
+        world_pix_new = SkyCoord(linear_wcs.all_pix2world(pix_inds, 1), unit='deg')
     
         dists = world_pix.separation(world_pix_new).to('deg')
         sigma = np.sqrt(sum(dists.value**2))
diff --git a/astrocut/tests/test_cube_cut.py b/astrocut/tests/test_cube_cut.py
index 223ffb88..1e489a29 100644
--- a/astrocut/tests/test_cube_cut.py
+++ b/astrocut/tests/test_cube_cut.py
@@ -217,11 +217,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,9) == 6.97e-07
-    assert round(sigma,7) == 1.4e-06
+    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
 
-    cry, crx =  myfactory.cutout_wcs.wcs.crpix
+    cry, crx = myfactory.cutout_wcs.wcs.crpix
     assert round(cry) == 4
     assert round(crx) == 2
     
diff --git a/astrocut/utils/wcs_fitting.py b/astrocut/utils/wcs_fitting.py
index 2a0985f8..1b8e7d76 100644
--- a/astrocut/utils/wcs_fitting.py
+++ b/astrocut/utils/wcs_fitting.py
@@ -13,14 +13,14 @@
 #
 #################################################################################
 
+# flake8: noqa
 
 import numpy as np
 
 from astropy import units as u
 from astropy.wcs.utils import celestial_frame_to_wcs
 
-
-def _linear_transformation_fit(params, lon, lat, x, y, w_obj):
+def _linear_transformation_fit(params, lon, lat, x, y, w_obj):  # pragma: no cover
     """ Objective function for fitting linear terms."""
     pc = params[0:4]
     crpix = params[4:6]
@@ -37,7 +37,7 @@ def _linear_transformation_fit(params, lon, lat, x, y, w_obj):
 
 
 def _sip_fit(params, lon, lat, u, v, w_obj, a_order, b_order, a_coeff_names,
-             b_coeff_names):
+             b_coeff_names):  # pragma: no cover
         """ Objective function for fitting SIP."""
         from astropy.modeling.models import SIP # here instead of top to avoid circular import
 
@@ -64,7 +64,7 @@ def _sip_fit(params, lon, lat, u, v, w_obj, a_order, b_order, a_coeff_names,
         return resids
 
 
-def fit_wcs_from_points(xp, yp, coords, mode, projection='TAN', proj_point=None, order=None, inwcs=None):
+def fit_wcs_from_points(xp, yp, coords, mode, projection='TAN', proj_point=None, order=None, inwcs=None):  # pragma: no cover
     """Given a set of matched x,y pixel positions and celestial coordinates, solves for
     WCS parameters and returns a WCS with these best fit values and other keywords based
     on projection type, frame, and units.