diff --git a/development/update_tess_sectors.ipynb b/development/update_tess_sectors.ipynb
new file mode 100644
index 0000000..d6b237e
--- /dev/null
+++ b/development/update_tess_sectors.ipynb
@@ -0,0 +1,238 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "f699122a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import pandas as pd \n",
+ "from astropy.time import Time"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "07bebdbf",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/var/folders/db/hdghk6ts5g11hr10jq0ss625nf1ny2/T/ipykernel_10556/1019568273.py:1: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support skipfooter; you can avoid this warning by specifying engine='python'.\n",
+ " data = pd.read_csv('https://tess.mit.edu/public/files/TESS_orbit_times.csv',skipfooter=1)\n"
+ ]
+ }
+ ],
+ "source": [
+ "data = pd.read_csv('https://tess.mit.edu/public/files/TESS_orbit_times.csv',skipfooter=1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "3ebad035",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sectors = data['Sector'].unique()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "id": "6d49e9f4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "df = pd.DataFrame()\n",
+ "for sector in sectors:\n",
+ " s = data.loc[data['Sector'] == sector]\n",
+ " start = Time(s['Start of Orbit'].values[0])\n",
+ " end = Time(s['End of Orbit'].values[1])\n",
+ " entry = pd.DataFrame()\n",
+ " entry['Sector'] = [sector]\n",
+ " entry['mjd_start'] = [start.mjd]\n",
+ " entry['mjd_end'] = [end.mjd]\n",
+ " df = pd.concat([df,entry])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "id": "dfc55f06",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " Sector | \n",
+ " mjd_start | \n",
+ " mjd_end | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " 1 | \n",
+ " 58324.815972 | \n",
+ " 58352.684028 | \n",
+ "
\n",
+ " \n",
+ " | 0 | \n",
+ " 2 | \n",
+ " 58353.607639 | \n",
+ " 58381.020833 | \n",
+ "
\n",
+ " \n",
+ " | 0 | \n",
+ " 3 | \n",
+ " 58382.225694 | \n",
+ " 58408.875000 | \n",
+ "
\n",
+ " \n",
+ " | 0 | \n",
+ " 4 | \n",
+ " 58410.406250 | \n",
+ " 58436.357639 | \n",
+ "
\n",
+ " \n",
+ " | 0 | \n",
+ " 5 | \n",
+ " 58437.482639 | \n",
+ " 58463.795139 | \n",
+ "
\n",
+ " \n",
+ " | ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ "
\n",
+ " \n",
+ " | 0 | \n",
+ " 93 | \n",
+ " 60829.357639 | \n",
+ " 60842.638889 | \n",
+ "
\n",
+ " \n",
+ " | 0 | \n",
+ " 94 | \n",
+ " 60855.763889 | \n",
+ " 60868.236111 | \n",
+ "
\n",
+ " \n",
+ " | 0 | \n",
+ " 95 | \n",
+ " 60881.833333 | \n",
+ " 60894.173611 | \n",
+ "
\n",
+ " \n",
+ " | 0 | \n",
+ " 96 | \n",
+ " 60907.274306 | \n",
+ " 60919.989583 | \n",
+ "
\n",
+ " \n",
+ " | 0 | \n",
+ " 97 | \n",
+ " 60933.159722 | \n",
+ " 60946.371528 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
97 rows × 3 columns
\n",
+ "
')[-1].split('<')[0])
max_t = deepcopy(disc_t)
- c = SkyCoord(ra,dec, unit=(u.hourangle, u.deg))
- ra = c.ra.deg
- dec = c.dec.deg
+ # c = SkyCoord(ra,dec, unit=(u.hourangle, u.deg))
+ # ra = c.ra.deg
+ # dec = c.dec.deg
- c = SkyCoord(ra,dec, unit=(u.hourangle, u.deg))
+ c = SkyCoord(ra,dec, unit=(u.deg, u.deg))
ra = c.ra.deg
dec = c.dec.deg
@@ -713,17 +741,17 @@ def par_psf_initialise(flux,camera,ccd,sector,column,row,cutoutSize,loc,time_ind
prf = create_psf(prf,cutoutSize)
return prf, cutout
-def par_psf_flux(image,prf,shift=[0,0],bkg_poly_order=3,kernel=None):
+def par_psf_flux(image,error,prf,shift=[0,0],bkg_poly_order=3,kernel=None):
if np.isnan(shift)[0]:
shift = np.array([0,0])
- prf.psf_flux(image,ext_shift=shift,poly_order=bkg_poly_order,kernel=kernel)
+ prf.psf_flux(image,error,ext_shift=shift,poly_order=bkg_poly_order,kernel=kernel)
return prf.flux, prf.eflux
-def par_psf_full(cutout,prf,shift=[0,0],xlim=0.5,ylim=0.5):
+def par_psf_full(cutout,error,prf,shift=[0,0],xlim=0.5,ylim=0.5):
if np.isnan(shift)[0]:
shift = np.array([0,0])
- prf.psf_position(cutout,ext_shift=shift,limx=xlim,limy=ylim)
- prf.psf_flux(cutout)
+ prf.psf_position(cutout,error,ext_shift=shift,limx=xlim,limy=ylim)
+ prf.psf_flux(cutout,error)
pos = [prf.source_x, prf.source_y]
return prf.flux, prf.eflux, pos
@@ -1207,3 +1235,130 @@ def grad_clip_fill_bkg(bkg,sigma=3,max_size=1000):
estimate = inpaint.inpaint_biharmonic(data,points)
return estimate
+def grad_clip(data,box_size=100):
+ """
+ Perform a local sigma clip of points based on the gradient of the points.
+ Pixels with large gradients are contaminated by stars/galaxies.
+
+ Inputs
+ ------
+ data : array
+ 1d array of the data to clip
+ box_size : int
+ integer defining the box size to clip over
+ Output
+ ------
+ gradind : bool
+
+ """
+ gradind = np.zeros_like(data)
+
+ for i in range(len(data)):
+ if i < box_size//2:
+ d = data[:i+box_size//2]
+ elif len(data) - i < box_size//2:
+ d = data[i-box_size//2:]
+ else:
+ d = data[i-box_size//2:i+box_size//2]
+
+ ind = np.isfinite(d)
+ d = d[ind]
+ if len(d) > 5:
+ gind = ~sigma_clip(np.gradient(abs(d))+d,sigma_lower=5,).mask
+ if i < box_size//2:
+ gradind[:i+box_size//2][ind] = gind
+ elif len(data) - i < box_size//2:
+ gradind[i-box_size//2:][ind] = gind
+ else:
+ gradind[i-box_size//2:i+box_size//2][ind] = gind
+
+ gradind = gradind > 0
+ return gradind
+
+def fit_strap(data,mask):
+ """
+ interpolate over missing data
+
+ """
+ x = np.arange(0,len(data))
+ y = data.copy()
+ p = np.ones_like(x)
+ if len(y[mask]) > 5:
+ p = interp1d(x[mask], y[mask],bounds_error=False,fill_value=np.nan,kind='linear')
+ p = p(x)
+ if np.isnan(p).any():
+ p2 = interp1d(x[mask], y[mask],bounds_error=False,fill_value='extrapolate',kind='nearest')
+ p2 = p2(x)
+ p[np.isnan(p)] = p2[np.isnan(p)]
+ return p
+
+def parallel_strap_fit(frame,frame_bkg,frame_err,mask,repeats=3,tol=3):
+ norm = frame / frame_bkg
+ sind = np.where(np.nansum(mask,axis=0)>0)[0]
+ qe = np.ones_like(frame)
+ for i in sind:
+ y = frame[:,i]
+ d = abs(np.gradient(y))
+ m, med, std = sigma_clipped_stats(d,maxiters=10)
+ nm = (d > med + std) * 1
+ nm = np.convolve(nm,np.ones(3),mode='same')
+ nm = (nm == 0)
+ #for r in range(repeats):
+ q = fit_strap(norm[:,i],nm)
+ #q /= frame_bkg[:,i]
+ if len(y) > 110:
+ q = savgol_filter(q,101,1)
+ else:
+ mq = np.nanmedian(q)
+ q[:] = mq
+
+ qe[:,i] = q
+
+ #qe[:,np.nanmean((qe-1),axis=) < 5e-3] = 1
+ return qe
+
+
+def simulate_epsf(camera,ccd,sector,column,row,size=13,realizations=2000,oversample=3):
+ from photutils.psf import EPSFBuilder, EPSFStar, EPSFStars
+
+ psfs = []
+ randpos = np.random.uniform(-1,1,size=(2,realizations)) + 6
+ prf_mod = TESS_PRF(cam=camera, ccd=ccd,sector=sector,
+ colnum=column,rownum=row)
+
+ for i in range(len(randpos)):
+ psfs += [prf_mod.locate(randpos[0,i],randpos[1,i],(13,13))]
+ psfs = np.array(psfs)
+
+ stars = []
+ for i in range(len(psfs)):
+ stars += [EPSFStar(psfs[i],cutout_center=randpos[:,i])]
+ stars = EPSFStars(stars)
+
+ epsf_builder = EPSFBuilder(oversampling=oversample, maxiters=20, progress_bar=False)
+ epsf, fitted_stars = epsf_builder(stars)
+
+ return epsf
+
+def parallel_photutils(cutout,e_cutout,psf_phot,init_params=None,return_pos=False):
+ if np.nansum(abs(cutout)) > 0:
+ phot = psf_phot(cutout, error=e_cutout,init_params=init_params)
+ phot = phot.to_pandas()
+ f = phot[['flux_fit']].values[0]; ef = phot[['flux_err']].values[0]
+ pos = phot[['x_fit','y_fit']].values
+ epos = phot[['x_err','y_err']].values
+ if return_pos:
+ return f, ef, pos, epos
+ else:
+ return f, ef
+ #return float(phot[['flux_fit']].values), float(phot[['flux_err']].values)
+ else:
+ if return_pos:
+ pos = np.array([np.nan,np.nan])
+ return np.array([np.nan]), np.array([np.nan]), pos, pos
+ else:
+ return np.array([np.nan]), np.array([np.nan])
+
+
+
+
diff --git a/tessreduce/psf_photom.py b/tessreduce/psf_photom.py
index f344513..d2e8073 100644
--- a/tessreduce/psf_photom.py
+++ b/tessreduce/psf_photom.py
@@ -91,7 +91,7 @@ def source(self,shiftx=0,shifty=0,ext_shift=[0,0]):
psf = shift(psf,ext_shift)
self.psf = psf/np.nansum(psf)
- def minimize_position(self,coeff,image,ext_shift):
+ def minimize_position(self,coeff,image,error,ext_shift):#,surface=True,order=2):
"""
Applies an exponential function using psf residuals to optimise the psf fit.
@@ -114,7 +114,15 @@ def minimize_position(self,coeff,image,ext_shift):
optimization model used for psf fitting
"""
-
+ # if surface:
+ # x = np.arange(image.shape[1])
+ # y = np.arange(image.shape[0])
+ # yy,xx = np.meshgrid(y,x)
+ # plane_coeff = coeff[2:]
+ # s = polynomial_surface(xx,yy,plane_coeff,order)
+ # else:
+ # s = 0
+
self.source_x = coeff[0]
self.source_y = coeff[1]
@@ -122,11 +130,11 @@ def minimize_position(self,coeff,image,ext_shift):
self.source(shiftx = self.source_x, shifty = self.source_y,ext_shift=ext_shift)
# -- calculate residuals -- #
- diff = abs(image - self.psf)
- residual = np.nansum(diff**2)
+ diff = abs(image - self.psf)**2
+ residual = np.nansum(diff/error)
return residual#np.exp(residual)
- def psf_position(self,image,limx=0.5,limy=0.5,ext_shift=[0,0]):
+ def psf_position(self,image,error,limx=0.8,limy=0.8,ext_shift=[0,0]):#,surface=False,order=2):
"""
Finds the optimal psf fit
@@ -150,17 +158,33 @@ def psf_position(self,image,limx=0.5,limy=0.5,ext_shift=[0,0]):
Optimal psf fit to input image
"""
+ if error is None:
+ error = np.ones_like(image)
+ #brightloc =
+ if np.nansum(image) > 0:
+ if np.isfinite(ext_shift).all():
+ ext_shift[0] = 0; ext_shift[1] = 0
+ normimage = image / np.nansum(image) # normalise the image
+
+ # if surface:
+ # num_coeffs = (poly_order + 1) * (poly_order + 2) // 2
+ # coeff = np.zeros(num_coeffs + 2)
+ # coeff[0] = self.source_x; coeff[1] = self.source_y
+ # lims = [[-limx,limx],[-limy,limy]]
+ # for i in range(num_coeffs):
+ # lims += [-np.inf,np.inf]
+ # else:
+ coeff = [self.source_x,self.source_y]
+ lims = [[-limx,limx],[-limy,limy]]
+
+ # -- Optimize -- #
+ res = minimize(self.minimize_position, coeff, args=(normimage,error,ext_shift), method='Powell',bounds=lims)
+ print('Optimal PSF shift: ', res.x)
+ self.psf_fit = res
+ else:
+ self.psf_fit = None
- normimage = image / np.nansum(image) # normalise the image
- coeff = [self.source_x,self.source_y]
- lims = [[-limx,limx],[-limy,limy]]
-
- # -- Optimize -- #
- res = minimize(self.minimize_position, coeff, args=(normimage,ext_shift), method='Powell',bounds=lims)
- print(res.x)
- self.psf_fit = res
-
- def minimize_psf_flux(self,coeff,image,surface=True,order=2,kernel=None):
+ def minimize_psf_flux(self,coeff,image,error=None,surface=True,order=2,kernel=None):
"""
@@ -194,10 +218,10 @@ def minimize_psf_flux(self,coeff,image,surface=True,order=2,kernel=None):
if kernel is not None:
self.psf = fftconvolve(self.psf, kernel, mode='same')
- res = np.nansum((image - self.psf*coeff[0] - s)**2)
+ res = np.nansum((image - self.psf*coeff[0] - s)**2/error)
return res
- def psf_flux(self,image,ext_shift=None,surface=True,poly_order=3,kernel=None):
+ def psf_flux(self,image,error=None,ext_shift=None,surface=True,poly_order=3,kernel=None):
"""
@@ -222,10 +246,11 @@ def psf_flux(self,image,ext_shift=None,surface=True,poly_order=3,kernel=None):
residual of optimal psf flux fit
"""
-
+ if error is None:
+ error = np.ones_like(image)
if self.psf is None:
self.source(shiftx=self.source_x,shifty=self.source_y)
- if ext_shift is not None:
+ if (ext_shift is not None) & np.isfinite(ext_shift).all():
self.source(ext_shift=ext_shift)
mask = np.zeros_like(self.psf)
mask[self.psf > np.nanpercentile(self.psf,90)] = 1
@@ -240,8 +265,10 @@ def psf_flux(self,image,ext_shift=None,surface=True,poly_order=3,kernel=None):
else:
initial = f0
- res = minimize(self.minimize_psf_flux,initial,args=(image,surface,poly_order,kernel),method='BFGS')
+ res = minimize(self.minimize_psf_flux,initial,args=(image,error,surface,poly_order,kernel),method='BFGS')
+ #res = minimize(self.minimize_psf_flux,initial,args=(image,error,surface,poly_order,kernel),method='Powell')
error = np.sqrt(np.diag(res['hess_inv']))
+
self.res = res
self.flux = res.x[0]
self.eflux = error[0]
diff --git a/tessreduce/sector_mjd.csv b/tessreduce/sector_mjd.csv
index 3367f7f..0effd67 100755
--- a/tessreduce/sector_mjd.csv
+++ b/tessreduce/sector_mjd.csv
@@ -1,84 +1,98 @@
-Sector,Dates,mjd_start,mjd_end
-1,07/25/18-08/22/18,58324,58352
-2,08/22/18-09/20/18,58352,58381
-3,09/20/18-10/18/18,58381,58409
-4,10/18/18-11/15/18,58409,58437
-5,11/15/18-12/11/18,58437,58463
-6,12/11/18-01/07/19,58463,58490
-7,01/07/19-02/02/19,58490,58516
-8,02/02/19-02/28/19,58516,58542
-9,02/28/19-03/26/19,58542,58568
-10,03/26/19-04/22/19,58568,58595
-11,04/22/19-05/21/19,58595,58624
-12,05/21/19-06/19/19,58624,58653
-13,06/19/19-07/18/19,58653,58682
-14,07/18/19-08/15/19,58682,58710
-15,08/15/19-09/11/19,58710,58737
-16,09/11/19-10/07/19,58737,58763
-17,10/07/19-11/02/19,58763,58789
-18,11/02/19-11/27/19,58789,58814
-19,11/27/19-12/24/19,58814,58841
-20,12/24/19-01/21/20,58841,58869
-21,01/21/20-02/18/20,58869,58897
-22,02/18/20-03/18/20,58897,58926
-23,03/18/20-04/16/20,58926,58955
-24,04/16/20-05/13/20,58955,58982
-25,05/13/20-06/08/20,58982,59008
-26,06/08/20-07/04/20,59008,59034
-27,07/04/20-07/30/20,59034,59060
-28,07/30/20-08/26/20,59060,59087
-29,08/26/20-09/22/20,59087,59114
-30,09/22/20-10/21/20,59114,59143
-31,10/21/20-11/19/20,59143,59172
-32,11/19/20-12/17/20,59172,59200
-33,12/17/20-01/13/21,59200,59227
-34,01/13/21-02/09/21,59227,59254
-35,02/09/21-03/07/21,59254,59280
-36,03/07/21-04/02/21,59280,59306
-37,04/02/21-04/28/21,59306,59332
-38,04/28/21-05/26/21,59332,59360
-39,05/26/21-06/24/21,59360,59389
-40,06/24/21-07/23/21,59389,59418
-41,07/23/21-08/20/21,59418,59446
-42,08/20/21-09/16/21,59446,59473
-43,09/16/21-10/12/21,59473,59499
-44,10/12/21-11/06/21,59499,59524
-45,11/06/21-12/02/21,59524,59550
-46,12/02/21-12/30/21,59550,59578
-47,12/30/21-01/28/22,59578,59607
-48,01/28/22-02/26/22,59607,59636
-49,02/26/22-03/26/22,59636,59664
-50,03/26/22-04/22/22,59664,59691
-51,04/22/22-05/18/22,59691,59717
-52,05/18/22-06/13/22,59717,59743
-53,06/13/22-07/09/22,59743,59769
-54,07/09/22-08/05/22,59769,59796
-55,08/05/22-09/01/22,59796,59823
-56,09/01/22-09/30/22,59823,59852
-57,09/30/22-10/29/22,59852,59881
-58,10/29/22-11/26/22,59881,59909
-59,11/26/22-12/23/22,59909,59936
-60,12/23/22-01/18/23,59936,59962
-61,01/18/23-02/12/23,59962,59987
-62,02/12/23-03/10/23,59987,60013
-63,03/10/23-04/06/23,60013,60040
-64,04/06/23-05/04/23,60040,60068
-65,05/04/23-06/02/23,60068,60097
-66,06/02/23-07/01/23,60097,60126
-67,07/01/23-07/29/23,60126,60154
-68,07/29/23-08/25/23,60154,60181
-69,08/25/23-09/20/23,60181,60207
-70,20/09/23-16/10/23,60207,60233
-71,16/10/23-11/11/23,60233,60259
-72,11/11/23-07/12/23,60259,60285
-73,07/12/23-03/01/24,60285,60312
-74,03/01/24-30/01/24,60312,60339
-75,30/01/24-26/02/24,60339,60366
-76,26/02/24-26/03/24,60366,60395
-77,26/03/24-23/04/24,60395,60423
-78,23/04/24-21/05/24,60423,60451
-79,21/05/24-18/06/24,60451,60479
-80,18/06/24-15/07/24,60479,60506
-81,15/07/24-10/08/24,60506,60532
-82,10/08/24-05/09/24,60532,60558
-83,05/09/24-01/10/24,60558,60584
\ No newline at end of file
+Sector,mjd_start,mjd_end
+1,58324.81597222222,58352.68402777778
+2,58353.60763888889,58381.020833333336
+3,58382.225694444445,58408.875
+4,58410.40625,58436.35763888889
+5,58437.48263888889,58463.79513888889
+6,58464.711805555555,58489.552083333336
+7,58491.131944444445,58515.59375
+8,58516.84722222222,58541.506944444445
+9,58542.72222222222,58567.98263888889
+10,58568.9375,58595.1875
+11,58596.27777777778,58623.399305555555
+12,58624.45486111111,58652.399305555555
+13,58653.42013888889,58681.864583333336
+14,58682.854166666664,58709.711805555555
+15,58710.864583333336,58736.916666666664
+16,58738.15277777778,58762.82638888889
+17,58764.18402777778,58789.20138888889
+18,58790.15625,58814.538194444445
+19,58815.583333333336,58840.65625
+20,58842.00347222222,58868.32986111111
+21,58869.93402777778,58897.288194444445
+22,58898.805555555555,58926.0
+23,58927.604166666664,58954.381944444445
+24,58955.29513888889,58981.788194444445
+25,58983.131944444445,59008.8125
+26,59009.76736111111,59034.64236111111
+27,59035.77777777778,59060.149305555555
+28,59061.350694444445,59086.604166666664
+29,59087.739583333336,59113.94097222222
+30,59115.385416666664,59142.729166666664
+31,59144.01388888889,59171.53472222222
+32,59172.572916666664,59199.739583333336
+33,59201.23263888889,59227.07986111111
+34,59228.25,59253.572916666664
+35,59254.489583333336,59279.48611111111
+36,59280.40277777778,59305.49652777778
+37,59306.739583333336,59332.086805555555
+38,59333.354166666664,59360.05902777778
+39,59361.270833333336,59389.225694444445
+40,59390.15277777778,59418.36111111111
+41,59419.489583333336,59446.086805555555
+42,59447.19097222222,59472.666666666664
+43,59473.666666666664,59498.395833333336
+44,59499.6875,59523.947916666664
+45,59525.006944444445,59550.131944444445
+46,59551.06597222222,59578.27777777778
+47,59579.302083333336,59606.447916666664
+48,59607.43402777778,59635.493055555555
+49,59636.97222222222,59663.822916666664
+50,59664.770833333336,59691.01736111111
+51,59692.447916666664,59717.041666666664
+52,59718.135416666664,59742.583333333336
+53,59743.49652777778,59768.48611111111
+54,59769.399305555555,59795.635416666664
+55,59796.600694444445,59823.770833333336
+56,59824.756944444445,59838.020833333336
+57,59852.854166666664,59866.75347222222
+58,59881.82986111111,59895.5625
+59,59909.76388888889,59922.854166666664
+60,59936.40277777778,59949.274305555555
+61,59962.29861111111,59974.70486111111
+62,59987.94097222222,60000.430555555555
+63,60013.868055555555,60026.899305555555
+64,60040.614583333336,60054.586805555555
+65,60068.239583333336,60082.98263888889
+66,60097.17361111111,60111.48263888889
+67,60126.14236111111,60139.756944444445
+68,60154.11111111111,60167.77777777778
+69,60181.854166666664,60194.5
+70,60207.854166666664,60220.631944444445
+71,60233.53472222222,60246.51388888889
+72,60259.68402777778,60272.506944444445
+73,60285.29861111111,60298.447916666664
+74,60312.364583333336,60325.354166666664
+75,60339.28472222222,60353.0625
+76,60367.197916666664,60380.89236111111
+77,60394.98611111111,60408.774305555555
+78,60423.26388888889,60437.166666666664
+79,60452.038194444445,60465.10763888889
+80,60479.38888888889,60492.850694444445
+81,60506.052083333336,60519.07986111111
+82,60532.89236111111,60544.96875
+83,60558.927083333336,60571.21527777778
+84,60584.09027777778,60596.805555555555
+85,60610.055555555555,60622.739583333336
+86,60635.760416666664,60649.114583333336
+87,60662.541666666664,60675.947916666664
+88,60689.65277777778,60703.395833333336
+89,60717.63888888889,60732.45138888889
+90,60746.65972222222,60760.239583333336
+91,60774.791666666664,60788.76388888889
+92,60802.47222222222,60815.98611111111
+93,60829.35763888889,60842.63888888889
+94,60855.76388888889,60868.23611111111
+95,60881.833333333336,60894.17361111111
+96,60907.274305555555,60919.989583333336
+97,60933.15972222222,60946.37152777778
diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py
index 4d40336..d91b451 100644
--- a/tessreduce/tessreduce.py
+++ b/tessreduce/tessreduce.py
@@ -56,7 +56,7 @@
package_directory = os.path.dirname(os.path.abspath(__file__)) + '/'
fig_width_pt = 240.0 # Get this from LaTeX using \showthe\columnwidth
-inches_per_pt = 1.0/72.27 # Convert pt to inches
+inches_per_pt = 1.0/72.27 # Convert pt to inches
golden_mean = (np.sqrt(5)-1.0)/2.0 # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt # width in inches
@@ -65,8 +65,8 @@ class tessreduce():
def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sector=None,
reduce=True,align=True,diff=True,corr_correction=True,kernel_match=False,calibrate=True,sourcehunt=True,
phot_method='aperture',imaging=False,parallel=True,num_cores=-1,diagnostic_plot=False,plot=True,
- savename=None,quality_bitmask='default',cache_dir=None,catalogue_path=False,
- shift_method='difference',prf_path=None,verbose=1):
+ savename=None,quality_bitmask='default',cache_dir=None,cache=True,catalogue_path=False,
+ shift_method='difference',use_error_image=False,prf_path=None,verbose=1,col_offset=0):
"""
Class for extracting reduced TESS photometry around a target coordinate or event.
@@ -145,6 +145,7 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect
self.kernel_match = kernel_match
self.imaging = imaging
self.parallel = parallel
+ self._col_offset = col_offset
if type(num_cores) == str:
self.num_cores = multiprocessing.cpu_count()
else:
@@ -153,6 +154,7 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect
self._sourcehunt = sourcehunt
self.verbose = verbose
self._shift_method = shift_method
+ self._use_error_image = use_error_image
# Offline Paths
if catalogue_path is None:
@@ -171,6 +173,7 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect
# Calculated
self.mask = None
+ #self._over_sub = None
self.shift = None
self.bkg = None
self.flux = None
@@ -186,6 +189,7 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect
self.zp_e = None
self.sn_name = None
self.ebv = 0
+ self.epsf = None
# repeat for backup
self.tzp = None
self.tzp_e = None
@@ -223,19 +227,24 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect
if type(tpf) == str:
self.tpf = lk.TessTargetPixelFile(tpf)
self.flux = strip_units(self.tpf.flux)
+ if self._use_error_image:
+ self.eflux = strip_units(self.tpf.flux_err)
+ else:
+ self.eflux = None
+ self.flux[np.isnan(self.flux)] = 0
self.wcs = self.tpf.wcs
self.ra = self.tpf.ra
self.dec = self.tpf.dec
self.size = self.tpf.flux.shape[1]
self.sector = self.tpf.sector
if self.sector is None:
- self.sector = 30
+ self.sector = 999
# Retrieve TPF
elif self.check_coord():
if self.verbose>0:
- print('getting TPF from TESScut')
- self.get_TESS(quality_bitmask=quality_bitmask,cache_dir=cache_dir)
+ print('Downloading TPF from TESScut')
+ self.get_TESS(quality_bitmask=quality_bitmask,cache_dir=cache_dir,cache=cache)
self._get_gaia()
self.ground = ground(ra = self.ra, dec = self.dec)
@@ -278,9 +287,9 @@ def _get_gaia(self,maglim=21):
result = Get_Catalogue(self.tpf, Catalog = 'gaia')
result = result[result.Gmag < maglim]
result = result.rename(columns={'RA_ICRS': 'ra',
- 'DE_ICRS': 'dec',
- 'e_RA_ICRS': 'e_ra',
- 'e_DE_ICRS': 'e_dec',})
+ 'DE_ICRS': 'dec',
+ 'e_RA_ICRS': 'e_ra',
+ 'e_DE_ICRS': 'e_dec',})
# Convert star RA/DEC to pixel values and input into dataframe
x,y = self.wcs.all_world2pix(result['ra'].values,result['dec'].values,0)
@@ -325,7 +334,13 @@ def _assign_phot_method(self,phot_method):
m = 'phot_mehtod must be a string equal to either "psf", or "aperture".'
raise ValueError(m)
- def get_TESS(self,ra=None,dec=None,name=None,size=None,sector=None,quality_bitmask='default',cache_dir=None):
+ def __clean_lk_cache(self,cache_dir=None):
+ if cache_dir is None:
+ cache = lk.config.get_cache_dir()
+
+
+ def get_TESS(self,ra=None,dec=None,name=None,size=None,sector=None,
+ quality_bitmask='default',cache_dir=None,cache=True):
"""
Use the lightcurve interface with TESScut to get an FFI cutout
of a region around the given coords.
@@ -382,6 +397,14 @@ def get_TESS(self,ra=None,dec=None,name=None,size=None,sector=None,quality_bitma
# Download
tpf = tess.download(quality_bitmask=quality_bitmask,cutout_size=size,download_dir=cache_dir)
+ if not cache:
+ try:
+ call = f'rm {tpf.path}'
+ os.system(call)
+ if self.verbose > 0:
+ print('Cache removed')
+ except:
+ print(f'Failed to remove: {tpf.path}')
# Check to ensure it succeeded
if tpf is None:
@@ -390,6 +413,11 @@ def get_TESS(self,ra=None,dec=None,name=None,size=None,sector=None,quality_bitma
self.tpf = tpf
self.flux = strip_units(tpf.flux) # Stripping astropy units so only numbers are returned
+ self.flux[np.isnan(self.flux)] = 0
+ if self._use_error_image:
+ self.eflux = strip_units(tpf.flux_err)
+ else:
+ self.eflux = None
self.wcs = tpf.wcs
def make_mask(self,catalogue_path=None,maglim=19,scale=1,strapsize=6,useref=False):
@@ -429,9 +457,9 @@ def make_mask(self,catalogue_path=None,maglim=19,scale=1,strapsize=6,useref=Fals
# Generate mask from source catalogue
if useref:
- mask, cat = Cat_mask(self.tpf,catalogue_path,maglim,scale,strapsize,ref=self.ref)
+ mask, cat = Cat_mask(self.tpf,catalogue_path,maglim,scale,strapsize,ref=self.ref,col_offset=self._col_offset)
else:
- mask, cat = Cat_mask(self.tpf,catalogue_path,maglim,scale,strapsize)
+ mask, cat = Cat_mask(self.tpf,catalogue_path,maglim,scale,strapsize,col_offset=self._col_offset)
# Generate sky background as the inverse of mask
sky = ((mask & 1)+1 ==1) * 1.
@@ -487,20 +515,30 @@ def psf_source_mask(self,sigma=5):
"""
+ col = self.tpf.column + int(self.size//2) # find column and row, when specifying location on a *say* 90x90 px cutout
+ row = self.tpf.row + int(self.size//2)
+
+ col += 45 # add on the non-science columns
+ row += 1 # add on the non-science row
+ if col > 2090:
+ col = 2090
+ if row > 2040:
+ row = 2040
+
# Find PRF for cutout (depends on Sector, Camera, CCD, Pixel Row, Pixel Column)
if self._catalogue_path is not None:
if self.sector < 4:
prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.sector,
- self.tpf.column+self.flux.shape[2]/2,self.tpf.row+self.flux.shape[1]/2,
+ col,row,
localdatadir=f'{self._prf_path}/Sectors1_2_3')
else:
prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.sector,
- self.tpf.column+self.flux.shape[2]/2,self.tpf.row+self.flux.shape[1]/2,
+ col,row,
localdatadir=f'{self._prf_path}/Sectors4+')
else:
prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.sector,
- self.tpf.column+self.flux.shape[2]/2,self.tpf.row+self.flux.shape[1]/2)
+ col,row)
self.prf = prf.locate(5,5,(11,11))
@@ -522,7 +560,62 @@ def psf_source_mask(self,sigma=5):
m[i] = par_psf_source_mask(data[i],self.prf,sigma)
return m * 1.0
- def background(self,gauss_smooth=2,calc_qe=True, strap_iso=True,source_hunt=False,interpolate=True):
+ def _calc_qe(self,flux,bkg_smth):#,flux_e):
+ '''
+ Calculate the effective quantum efficiency enhancement of the detector from scattered light.
+ '''
+ norm = flux / bkg_smth
+ straps = norm * ((self.mask & 4)>0)
+
+ straps[straps==0] = np.nan
+ m,med,std = sigma_clipped_stats(straps,axis=(1,2),maxiters=10,mask_value=np.nan)
+ masks = straps < (med + 2*std)[:,np.newaxis,np.newaxis]
+ m = (np.nansum(masks,axis=0) > 0) * 1.
+ m[m==0] = np.nan
+ ratio = flux / bkg_smth * m
+ ratio[ratio < 1] = np.nan
+ #m, med, std = sigma_clipped_stats(ratio,axis=1,sigma_upper=2)
+ qe_1d = np.nanpercentile(ratio,10,axis=1)
+ qe = np.ones_like(norm)
+ qe[:,:,:] = qe_1d[:,np.newaxis,:]
+ #qe[:,:,:] = med[:,np.newaxis,:]
+ qe[np.isnan(qe)] = 1
+ qe[qe<1] = 1
+ m, med, std = sigma_clipped_stats(qe,axis=0)
+ qe[:] = med
+ return qe
+
+ # straps[straps==0] = np.nan
+ # m,med,std = sigma_clipped_stats(straps,axis=(1,2),maxiters=10,mask_value=np.nan)
+ # masks = straps < (med + 2*std)[:,np.newaxis,np.newaxis]
+ # m = (np.nansum(masks,axis=0) > 0) * 1.
+ # m[m==0] = np.nan
+ # qe = np.ones_like(flux) * 1.
+ # if flux.shape[1] < 10000:
+ # ratio = (flux / bkg_smth) * m
+
+ # m, med, std = sigma_clipped_stats(ratio,axis=1,sigma_upper=2,sigma_lower=3)
+ # #qe_1d = np.nanpercentile(ratio,10,axis=1)
+
+ # #qe[:,:,:] = qe_1d[:,np.newaxis,:]
+ # qe[:,:,:] = med[:,np.newaxis,:]
+
+ # else:
+ # index = np.arange(len(flux),dtype=int)
+ # if self.parallel:
+ # qe = Parallel(n_jobs=self.num_cores)(delayed(parallel_strap_fit)(flux[i],bkg_smth[i],flux_e[i],m) for i in index)
+ # else:
+ # qe = []
+ # for i in index:
+ # qe += [parallel_strap_fit(flux[i],bkg_smth[i],flux_e[i],m)]
+ # qe = np.array(qe)
+ # #filler = np.nanmedian(qe,axis=1)
+
+ # qe[np.isnan(qe)] = 1
+ # qe[qe<1] = 1
+ # return qe
+
+ def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False,interpolate=True,rerun_negative=False):
"""
Calculate the temporal and spatial variation in the background.
@@ -569,6 +662,29 @@ def background(self,gauss_smooth=2,calc_qe=True, strap_iso=True,source_hunt=Fals
bkg_smth = np.zeros_like(flux) * np.nan
if self.parallel:
bkg_smth = Parallel(n_jobs=self.num_cores)(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m)
+ if rerun_negative:
+ if self._use_error_image:
+ over_sub = (deepcopy(self.flux) - bkg_smth) < -self.eflux # -0.5
+ else:
+ over_sub = (deepcopy(self.flux) - bkg_smth) < -0.5
+ over_sub = np.nansum(over_sub,axis=0) > 0
+ #self._over_sub = over_sub
+ #print('overshape ',over_sub.shape)
+ #print('m ',m.shape)
+ strap_mask = (self.mask & 4) > 0
+ if len(strap_mask.shape) == 3:
+ strap_mask = strap_mask[0]
+ if strap_iso:
+ over_sub[strap_mask] = 0
+ if source_hunt | (len(self.mask.shape) == 3):
+ m[:,over_sub[:,:]] = 1
+ else:
+ m[over_sub] = 1
+ self._bkgmask = m
+ bkg_smth = Parallel(n_jobs=self.num_cores)(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m)
+
+
+
else:
for i in range(flux.shape[0]):
bkg_smth[i] = Smooth_bkg((flux*m)[i],gauss_smooth,interpolate)
@@ -579,25 +695,10 @@ def background(self,gauss_smooth=2,calc_qe=True, strap_iso=True,source_hunt=Fals
# Calculate quantum efficiency
if calc_qe:
- strap = (self.mask == 4) * 1.0
- strap[strap==0] = np.nan
- # check if its a time varying mask
- if len(strap.shape) == 3:
- strap = strap[self.ref_ind]
- mask = ((self.mask & 1) == 0) * 1.0
- mask[mask==0] = np.nan
-
- data = strip_units(self.flux) * mask
- norm = self.flux / bkg_smth
- straps = norm * ((self.mask & 4)>0)
- limit = np.nanpercentile(straps,60,axis=1)
- straps[limit[:,np.newaxis,:] < straps] = np.nan
- straps[straps==0] = 1
-
- value = np.nanmedian(straps,axis=1)
- qe = np.ones_like(bkg_smth) * value[:,np.newaxis,:]
- bkg = bkg_smth * qe
+ self.bkg = bkg_smth
+ qe = self._calc_qe(flux,bkg_smth)#,self.eflux)
self.qe = qe
+ bkg = bkg_smth * qe
else:
bkg = np.array(bkg_smth)
self.bkg = bkg
@@ -657,7 +758,7 @@ def _bkg_round_3(self,iters=5):
dist_mask = convolve(dist_mask,kern) > 0
if self.parallel:
bkg_3 = Parallel(n_jobs=self.num_cores)(delayed(parallel_bkg3)(self.bkg[i],dist_mask[i])
- for i in np.arange(len(dist_mask)))
+ for i in np.arange(len(dist_mask)))
else:
bkg_3 = []
bkg_smth = np.zeros_like(dist_mask)
@@ -682,7 +783,7 @@ def _clip_background(self,sigma=5,ideal_size=90):
if self.parallel:
bkg_clip = Parallel(n_jobs=self.num_cores)(delayed(clip_background)(self.bkg[i],self.mask,sigma,ideal_size)
- for i in np.arange(len(self.bkg)))
+ for i in np.arange(len(self.bkg)))
else:
bkg_clip = []
for i in range(len(self.bkg)):
@@ -708,7 +809,7 @@ def _grad_bkg_clip(self,sigma=3,max_size=1000):
if self.parallel:
bkg_clip = Parallel(n_jobs=self.num_cores)(delayed(grad_clip_fill_bkg)(self.bkg[i],sigma,max_size)
- for i in np.arange(len(self.bkg)))
+ for i in np.arange(len(self.bkg)))
else:
bkg_clip = []
for i in range(len(self.bkg)):
@@ -798,7 +899,7 @@ def centroids_shifts_starfind(self,plot=None,savename=None):
mean, med, std = sigma_clipped_stats(m, sigma=3.0)
prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.tpf.sector,
- self.tpf.column+self.flux.shape[2]/2,self.tpf.row+self.flux.shape[1]/2)
+ self.tpf.column+self.flux.shape[2]/2,self.tpf.row+self.flux.shape[1]/2)
self.prf = prf.locate(5,5,(11,11))
finder = StarFinder(2*std,kernel=self.prf,exclude_border=True)
@@ -876,15 +977,21 @@ def fit_shift(self,smooth=True,plot=None,savename=None):
f = self.flux
m = self.ref.copy() * sources
m[m==0] = np.nan
+ #eref = self.eflux[self.ref_ind]
+
if self.parallel:
+ ind = np.arange(len(f))
shifts = Parallel(n_jobs=self.num_cores)(
- delayed(difference_shifts)(frame,m) for frame in f)
+ #delayed(difference_shifts)(f[i],m,self.eflux[i],eref) for i in ind)
+ delayed(difference_shifts)(f[i],m) for i in ind)
shifts = np.array(shifts)
else:
shifts = np.zeros((len(f),2)) * np.nan
for i in range(len(f)):
+ #shifts[i,:] = difference_shifts(f[i],m,self.eflux[i],eref)
shifts[i,:] = difference_shifts(f[i],m)
+ sraw = deepcopy(shifts)
shifts = Smooth_motion(shifts,self.tpf)
@@ -898,8 +1005,10 @@ def fit_shift(self,smooth=True,plot=None,savename=None):
ind = np.where(np.diff(t) > .5)[0]
shifts[ind,:] = np.nan
plt.figure(figsize=(1.5*fig_width,1*fig_width))
- plt.plot(t,shifts[:,0],'.',label='Row shift',alpha =0.5)
- plt.plot(t,shifts[:,1],'.',label='Col shift',alpha =0.5)
+ plt.plot(t,sraw[:,0],'.',label='Row shift',alpha = 0.5)
+ plt.plot(t,sraw[:,1],'.',label='Col shift',alpha = 0.5)
+ plt.plot(t,shifts[:,0],'-',label='Smoothed row shift')
+ plt.plot(t,shifts[:,1],'-',label='Smoothed col shift')
plt.ylabel('Shift (pixels)',fontsize=15)
plt.xlabel('Time (MJD)',fontsize=15)
plt.legend()
@@ -945,7 +1054,7 @@ def shift_images(self,median=False):
if median:
for i in range(len(shifted)):
if np.nansum(abs(shifted[i])) > 0:
- shifted[i] = shift(self.ref,[-self.shift[i,1],-self.shift[i,0]])
+ shifted[i] = shift(self.ref,[-self.shift[i,1],-self.shift[i,0]], mode='nearest',order=5)
self.flux -= shifted
else:
@@ -1071,7 +1180,7 @@ def check_trend(self,lc=None,limit=0.6):
p1 = pearsonr(lc[finite],self.shift[finite,0])[0]
p2 = pearsonr(lc[finite],self.shift[finite,1])[0]
if (abs(p1) > limit) | (abs(p2) > limit):
- if p1 > p2:
+ if abs(p1) > abs(p2):
m = 'x shift'
p = p1
else:
@@ -1151,12 +1260,12 @@ def diff_lc(self,time=None,x=None,y=None,ra=None,dec=None,tar_ap=3,
if (ra is not None) & (dec is not None) & (self.tpf is not None):
x,y = self.wcs.all_world2pix(ra,dec,0)
- x = int(x + 0.5)
- y = int(y + 0.5)
+ x = int(np.round(x,0))
+ y = int(np.round(y,0))
elif (x is None) & (y is None):
x,y = self.wcs.all_world2pix(self.ra,self.dec,0)
- x = int(x + 0.5)
- y = int(y + 0.5)
+ x = int(np.round(x,0))
+ y = int(np.round(y,0))
ap_tar = np.zeros_like(data[0])
ap_sky = np.zeros_like(data[0])
@@ -1186,7 +1295,10 @@ def diff_lc(self,time=None,x=None,y=None,ra=None,dec=None,tar_ap=3,
else:
tar = np.nansum((data+self.ref)*ap_tar,axis=(1,2))
tar -= sky_med * tar_ap**2
- tar_err = sky_std * tar_ap**2
+ if self._use_error_image:
+ tar_err = np.nansum((self.eflux)*ap_tar,axis=(1,2))#sky_std * tar_ap**2
+ else:
+ tar_err = sky_std * tar_ap**2
if phot_method == 'psf':
if psf_snap is None:
psf_snap = 'brightest'
@@ -1210,7 +1322,7 @@ def diff_lc(self,time=None,x=None,y=None,ra=None,dec=None,tar_ap=3,
if savename is not None:
plt.savefig(savename + '_diff_diag.pdf', bbox_inches = "tight")
- p = self.check_trend(lc=lc[1])
+ #p = self.check_trend(lc=lc[1])
return lc, sky
def dif_diag_plot(self,ap_tar,ap_sky,lc=None,sky=None,data=None):
@@ -1267,9 +1379,9 @@ def dif_diag_plot(self,ap_tar,ap_sky,lc=None,sky=None,data=None):
nonan1 = np.isfinite(d)
nonan2 = np.isfinite(d*ap)
plt.imshow(data[maxind],origin='lower',
- vmin=np.nanpercentile(d,16),
- vmax=np.nanpercentile(d[nonan2],80),
- aspect='auto')
+ vmin=np.nanpercentile(d,16),
+ vmax=np.nanpercentile(d[nonan2],80),
+ aspect='auto')
cbar = plt.colorbar()
cbar.set_label('$e^-/s$',fontsize=15)
plt.xlabel('Column',fontsize=15)
@@ -1547,20 +1659,38 @@ def _psf_initialise(self,cutoutSize,loc,time_ind=None,ref=False):
if time_ind is None:
time_ind = np.arange(0,len(self.flux))
+
+ col = self.tpf.column - int(self.size//2) + loc[0] # find column and row, when specifying location on a *say* 90x90 px cutout
+ row = self.tpf.row - int(self.size//2) + loc[1]
+
if isinstance(loc[0], (float, np.floating, np.float32, np.float64)):
- loc[0] = int(loc[0] + 0.5)
+ loc[0] = int(np.round(loc[0],0))
if isinstance(loc[1], (float, np.floating, np.float32, np.float64)):
- loc[1] = int(loc[1] + 0.5)
-
- col = self.tpf.column - int(self.size/2-1) + loc[0] # find column and row, when specifying location on a *say* 90x90 px cutout
- row = self.tpf.row - int(self.size/2-1) + loc[1]
-
- prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row) # initialise psf kernel
+ loc[1] = int(np.round(loc[1]))
+
+ col += 45 # add on the non-science columns
+ row += 1 # add on the non-science row
+ if col > 2090:
+ col = 2090
+ if row > 2040:
+ row = 2040
+ row = np.max([row,10])
+ col = np.max([col,45])
+ try:
+ prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row) # initialise psf kernel
+ except:
+ print(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row)
+ raise ValueError
if ref:
cutout = (self.flux+self.ref)[time_ind,loc[1]-cutoutSize//2:loc[1]+1+cutoutSize//2,loc[0]-cutoutSize//2:loc[0]+1+cutoutSize//2] # gather cutouts
else:
cutout = self.flux[time_ind,loc[1]-cutoutSize//2:loc[1]+1+cutoutSize//2,loc[0]-cutoutSize//2:loc[0]+1+cutoutSize//2] # gather cutouts
- return prf, cutout
+ if self._use_error_image:
+ ecutout = self.eflux[time_ind,loc[1]-cutoutSize//2:loc[1]+1+cutoutSize//2,loc[0]-cutoutSize//2:loc[0]+1+cutoutSize//2] # gather cutouts
+ else:
+ ecutout = np.ones_like(cutout) * 0.1
+ #else:
+ return prf, cutout, ecutout
def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2):
"""
@@ -1597,14 +1727,14 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2):
raise ValueError(m)
inds = np.arange(0,len(xpos))
if self.parallel:
- prfs, cutouts = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd,
- self.tpf.sector,self.tpf.column,self.tpf.row,
- size,[xpos[i],ypos[i]],time_ind) for i in inds))
+ prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd,
+ self.tpf.sector,self.tpf.column,self.tpf.row,
+ size,[xpos[i],ypos[i]],time_ind) for i in inds))
else:
prfs = []
cutouts = []
for i in range(len(time_ind)):
- prf, cutout = self._psf_initialise(size,[xpos[i],ypos[i]],time_ind=time_ind[i])
+ prf, cutout, ecutouts = self._psf_initialise(size,[xpos[i],ypos[i]],time_ind=time_ind[i])
prfs += [prf]
cutouts += [cutout]
cutouts = np.array(cutouts)
@@ -1623,6 +1753,96 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2):
pos[0,:] += xpos; pos[1,:] += ypos
return flux, pos
+ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None,
+ flux=None,eflux=None,ref=False,return_pos=False,
+ snap='brightest'):
+ from photutils.psf import PSFPhotometry
+ from astropy.table import Table
+ rad = size // 2
+
+ if flux is None:
+ flux = self.flux
+ if (flux.shape[1] < size) | (flux.shape[2] < size):
+ e = 'Image dimensions must be larger than the cutout size'
+ raise ValueError(e)
+ if eflux is None:
+ eflux = self.eflux
+
+ if (xPix is None) | (yPix is None):
+ xPix = flux.shape[2]//2
+ yPix = flux.shape[1]//2
+
+ if epsf is None:
+ if self.epsf is None:
+ col = self.tpf.column - int(self.size//2) + yPix # find column and row, when specifying location on a *say* 90x90 px cutout
+ row = self.tpf.row - int(self.size//2) + xPix
+ col += 45 # add on the non-science columns
+ row += 1 # add on the non-science row
+ if col > 2090:
+ col = 2090
+ if row > 2040:
+ row = 2040
+
+ self.epsf = simulate_epsf(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row)
+ epsf = self.epsf
+
+ if local_bkg:
+ localbkg_estimator = LocalBackground(1.5, 7, bkgstat)
+ else:
+ localbkg_estimator = None
+
+ if ref:
+ flux = flux + self.ref
+
+ cutouts = flux[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1]
+ ecutouts = eflux[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1]
+ fit_shape = (size, size)
+ psfphot = PSFPhotometry(epsf, fit_shape, finder=None,aperture_radius=1.5,
+ localbkg_estimator=localbkg_estimator,xy_bounds=(0.05))
+ init = Table()
+ init['x_init'] = [size//2]
+ init['y_init'] = [size//2]
+ if snap.lower() != 'all':
+ if snap.lower() == 'brightest':
+ weight = np.abs(np.nansum((cutouts/ecutouts)[:,1:4,1:4],axis=(1,2)))
+ weight[np.isnan(weight)] = 0
+ ind = np.argmax(weight)
+ rcut = cutouts[ind]
+ ecut = ecutouts[ind]
+ elif snap.lower() == 'ref':
+ rcut = self.ref[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1]
+ ecut = ecutouts[self.ref_ind]
+
+
+ phot = psfphot(rcut, error=ecut,init_params=init)
+ # fit with the best position
+ init2 = Table()
+ init['x_init'] = phot['x_fit']
+ init['y_init'] = phot['y_fit']
+ flux = np.zeros(len(self.flux)) * np.nan
+ eflux = np.zeros(len(self.flux)) * np.nan
+ psfphot2 = PSFPhotometry(epsf, fit_shape, finder=None,aperture_radius=1.5,
+ xy_bounds=(0.05),localbkg_estimator=localbkg_estimator)
+ f,ef = zip(*Parallel(n_jobs=self.num_cores)(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot2,init) for i in np.arange(len(cutouts))))
+ f = np.array(f).flatten()
+ ef = np.array(ef).flatten()
+ phot = phot.to_pandas()
+ pos = phot[['x_fit','y_fit']].values + np.array([xPix,yPix]) - size//2
+ epos = phot[['x_err','y_err']].values
+ else:
+ f,ef,pos,epos = zip(*Parallel(n_jobs=self.num_cores)(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot,init,True) for i in np.arange(len(cutouts))))
+ pos['x_fit'] += xPix - size//2
+ pos['y_fit'] += xPix - size//2
+
+ if return_pos:
+ return f, ef, pos, epos
+ else:
+ return f, ef
+
+
+
+
+
def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=False,diff=None,bkg_poly_order=3):
"""
Main PSF Photometry function
@@ -1669,12 +1889,12 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa
if type(snap) == str:
if snap == 'all':
- prf, cutouts = self._psf_initialise(size,(xPix,yPix),ref=(not diff)) # gather base PRF and the array of cutouts data
+ prf, cutouts, ecutouts = self._psf_initialise(size,(xPix,yPix),ref=(not diff)) # gather base PRF and the array of cutouts data
inds = np.arange(len(cutouts))
base = create_psf(prf,size)
flux, eflux, pos = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_full)(cutouts[i],base,self.shift[i]) for i in inds))
- #prf, cutouts = self._psf_initialise(size,(xPix,yPix)) # gather base PRF and the array of cutouts data
+ #prf, cutouts = self._psf_initialise(size,(xPix,yPix)) # gather base PRF and the array of cutouts data
#xShifts = []
#yShifts = []
#for cutout in tqdm(cutouts):
@@ -1694,53 +1914,54 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa
# ax[2].set_ylabel('yShift')
else:
if snap == 'brightest': # each cutout has position snapped to brightest frame fit position
- prf, cutouts = self._psf_initialise(size,(xPix,yPix),ref=(not diff)) # gather base PRF and the array of cutouts data
+ prf, cutouts, ecutouts = self._psf_initialise(size,(xPix,yPix),ref=(not diff)) # gather base PRF and the array of cutouts data
bkg = self.bkg[:,int(yPix),int(xPix)]
#lowbkg = bkg < np.nanpercentile(bkg,16)
- weight = np.abs(np.nansum(cutouts[:,int(yPix)-1:int(yPix)+2,int(xPix)-1:int(xPix)+2],axis=(1,2))) / bkg
+ #weight = np.abs(np.nansum(cutouts[:,int(yPix)-1:int(yPix)+2,int(xPix)-1:int(xPix)+2],axis=(1,2))) / bkg
+ weight = np.abs(np.nansum((cutouts / ecutouts)[:,int(yPix)-1:int(yPix)+2,int(xPix)-1:int(xPix)+2],axis=(1,2)))
weight[np.isnan(weight)] = 0
- ind = np.argmax(abs(weight))
+ ind = np.argmax(weight)
#ind = np.where(cutouts==np.nanmax(cutouts))[0][0]
ref = cutouts[ind]
base = create_psf(prf,size)
- base.psf_position(ref,ext_shift=self.shift[ind])
+ base.psf_position(ref,ecutouts[ind],ext_shift=self.shift[ind])
elif snap == 'ref':
- prf, cutouts = self._psf_initialise(size,(xPix,yPix),ref=True) # gather base PRF and the array of cutouts data
+ prf, cutouts,ecutouts = self._psf_initialise(size,(xPix,yPix),ref=True) # gather base PRF and the array of cutouts data
ref = cutouts[self.ref_ind]
base = create_psf(prf,size)
- base.psf_position(ref)
+ base.psf_position(ref,ecutouts[self.ref_ind])
if diff:
- _, cutouts = self._psf_initialise(size,(xPix,yPix),ref=False)
+ _, cutouts,ecutouts = self._psf_initialise(size,(xPix,yPix),ref=False)
elif snap == 'fixed':
- prf, cutouts = self._psf_initialise(size,(xPix,yPix),ref=(not diff)) # gather base PRF and the array of cutouts data
+ prf, cutouts, ecutouts = self._psf_initialise(size,(xPix,yPix),ref=(not diff)) # gather base PRF and the array of cutouts data
base = create_psf(prf,size)
if self.parallel:
inds = np.arange(len(cutouts))
if self.delta_kernel is not None:
- flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds))
+ flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds))
else:
- flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],base,self.shift[i],bkg_poly_order) for i in inds))
+ flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds))
else:
for i in range(len(cutouts)):
- flux += [par_psf_flux(cutouts[i],base,self.shift[i])]
+ flux += [par_psf_flux(cutouts[i],ecutouts[i],base,self.shift[i])]
if plot:
plt.figure()
plt.plot(flux)
plt.ylabel('Flux')
- elif type(snap) == int: # each cutout has position snapped to 'snap' frame fit position (snap is integer)
+ elif type(snap) == int: # each cutout has position snapped to 'snap' frame fit position (snap is integer)
base = create_psf(prf,size)
base.psf_position(cutouts[snap])
if self.parallel:
inds = np.arange(len(cutouts))
if self.delta_kernel is not None:
- flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds))
+ flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds))
else:
- flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],base,self.shift[i],bkg_poly_order) for i in inds))
+ flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds))
else:
for i in range(len(cutouts)):
- flux += [par_psf_flux(cutouts[i],base,self.shift[i])]
+ flux += [par_psf_flux(cutouts[i],ecutouts[i],base,self.shift[i])]
if plot:
fig,ax = plt.subplots(ncols=1,figsize=(12,4))
ax.plot(flux)
@@ -1855,19 +2076,24 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None,
#print('mask frac ',frac)
if frac < 0.05:
print('!!!WARNING!!! mask is too dense, lowering mask_scale to 0.5, and raising maglim to 15. Background quality will be reduced.')
- self.make_mask(catalogue_path=self._catalogue_path,maglim=15,strapsize=7,scale=0.5)
+ self.make_mask(catalogue_path=self._catalogue_path,maglim=15,strapsize=7,scale=0.2)
if self.verbose > 0:
print('made source mask')
else:
self.mask = mask
if self.verbose > 0:
print('assigned source mask')
+ if moving_mask is not None:
+ moving_mask = moving_mask > 0
+ temp = np.zeros_like(self.flux,dtype=int)
+ temp[:,:,:] = self.mask
+ self.mask = temp | moving_mask
# calculate background for each frame
if self.verbose > 0:
print('calculating background')
# calculate the background
- self.background()
+ self.background(rerun_negative=True)
if np.isnan(self.bkg).all():
@@ -1927,7 +2153,7 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None,
if self.verbose > 0:
print('!!Re-running for difference image!!')
# reseting to do diffim
- self._flux_aligned = deepcopy(self.flux)
+
self.flux = strip_units(self.tpf.flux)
self.flux = self.flux / self.qe
@@ -1936,7 +2162,7 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None,
if self.verbose > 0:
print('shifting images')
-
+ self._flux_aligned = deepcopy(self.flux)
if test_seed is not None:
self.flux += test_seed
self.flux[np.nansum(self.tpf.flux.value,axis=(1,2))==0] = np.nan
@@ -1945,13 +2171,14 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None,
self.flux -= self.ref
self.ref -= self.bkg[self.ref_ind]
+ self._ref_bkg = self.bkg[self.ref_ind]
# remake mask
- self.make_mask(catalogue_path=self._catalogue_path,maglim=18,strapsize=7,scale=mask_scale*.8,useref=False)#Source_mask(ref,grid=0)
+ self.make_mask(catalogue_path=self._catalogue_path,maglim=17,strapsize=7,scale=mask_scale*.5,useref=False)#Source_mask(ref,grid=0)
frac = np.nansum((self.mask== 0) * 1.) / (self.mask.shape[0] * self.mask.shape[1])
#print('mask frac ',frac)
if frac < 0.05:
print('!!!WARNING!!! mask is too dense, lowering mask_scale to 0.5, and raising maglim to 15. Background quality will be reduced.')
- self.make_mask(catalogue_path=self._catalogue_path,maglim=15,strapsize=7,scale=0.5)
+ self.make_mask(catalogue_path=self._catalogue_path,maglim=15,strapsize=7,scale=mask_scale*.5*.5)
# assuming that the target is in the centre, so masking it out
#m_tar = np.zeros_like(self.mask,dtype=int)
#m_tar[self.ref.shape[0]//2,self.ref.shape[1]//2]= 1
@@ -2052,7 +2279,7 @@ def make_lc(self,aperture = None,bin_size=0,zeropoint=None,scale='counts',clip =
elif type(aperture) == np.ndarray:
aper = aperture * 1.
- lc = Lightcurve(flux,aper) #,scale = scale)
+ lc = Lightcurve(flux,aper) #,scale = scale)
if clip:
mask = ~sigma_mask(lc)
lc[mask] = np.nan
@@ -2420,7 +2647,7 @@ def bin_interp(self,lc=None,time_bin=6/24):
smooth = f1(lc[0])
return smooth
- def detrend_star(self,lc=None):
+ def detrend_star(self,lc=None,normalize=False):
"""
Removes trends, e.g. background or stellar variability from the lightcurve data.
@@ -2467,7 +2694,10 @@ def detrend_star(self,lc=None):
smooth = f1(temp[0])
detrended = deepcopy(lc)
- detrended[1] -= smooth
+ if normalize:
+ detrended[1] /= smooth
+ else:
+ detrended[1] -= smooth
return detrended
@@ -2532,12 +2762,12 @@ def isolated_star_lcs(self):
ind = dist < 1.5
close = tab.iloc[ind]
- d['gmag'].iloc[i] = -2.5*np.log10(np.nansum(mag2flux(close.gmag.values,25))) + 25
- d['rmag'].iloc[i] = -2.5*np.log10(np.nansum(mag2flux(close.rmag.values,25))) + 25
- d['imag'].iloc[i] = -2.5*np.log10(np.nansum(mag2flux(close.imag.values,25))) + 25
- d['zmag'].iloc[i] = -2.5*np.log10(np.nansum(mag2flux(close.zmag.values,25))) + 25
+ d.loc[i,'gmag'] = -2.5*np.log10(np.nansum(maselfflux(close.gmag.values,25))) + 25
+ d.loc[i,'rmag'] = -2.5*np.log10(np.nansum(maselfflux(close.rmag.values,25))) + 25
+ d.loc[i,'imag'] = -2.5*np.log10(np.nansum(maselfflux(close.imag.values,25))) + 25
+ d.loc[i,'zmag'] = -2.5*np.log10(np.nansum(maselfflux(close.zmag.values,25))) + 25
if system == 'ps1':
- d['ymag'].iloc[i] = -2.5*np.log10(np.nansum(mag2flux(close.ymag.values,25))) + 25
+ d.loc[i,'ymag'] = -2.5*np.log10(np.nansum(maselfflux(close.ymag.values,25))) + 25
# convert to tess mags
if len(d) < 10:
print('!!!WARNING!!! field calibration is unreliable, using the default zp = 20.44')
@@ -2559,11 +2789,11 @@ def isolated_star_lcs(self):
for i in range(len(d)):
#if self.phot_method == 'aperture':
mask = np.zeros_like(self.ref)
- mask[int(d.row.values[i] + .5),int(d.col.values[i] + .5)] = 1
+ mask[int(np.round(d.row.values[i],0)),int(np.round(d.col.values[i],0))] = 1
mask = convolve(mask,np.ones((3,3)))
flux += [np.nansum(tflux*mask,axis=(1,2))]
m2 = np.zeros_like(self.ref)
- m2[int(d.row.values[i] + .5),int(d.col.values[i] + .5)] = 1
+ m2[int(np.round(d.row.values[i],0)),int(np.round(d.col.values[i]))] = 1
m2 = convolve(m2,np.ones((7,7))) - convolve(m2,np.ones((5,5)))
eflux += [np.nansum(tflux*m2,axis=(1,2))]
mag = -2.5*np.log10(np.nansum((self.ref*m2))) + 20.44
@@ -2618,7 +2848,6 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None):
Error in the photometric zeropoint
"""
-
if plot is None:
plot = self.diagnostic_plot
if savename is None:
@@ -2651,31 +2880,32 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None):
tflux = self.flux
- ind = (table.imag.values < 19) & (table.imag.values > 14)
+ ind = (table.imag.values < 19) & (table.imag.values > 13.5)
tab = table.iloc[ind]
e, dat = Tonry_reduce(tab,plot=plot,savename=savename,system=system)
self.ebv = e[0]
- gr = (dat.gmag - dat.rmag).values
- ind = (gr < 1) & (dat.imag.values < 15)
+ gr = (dat.gmag - dat['rmag']).values
+ ind = (gr < 1) & (dat['imag'].values < 16)
d = dat.iloc[ind]
x,y = self.wcs.all_world2pix(d.RAJ2000.values,d.DEJ2000.values,0)
d['col'] = x
d['row'] = y
- pos_ind = (5 < x) & (x < self.ref.shape[1]-5) & (5 < y) & (y < self.ref.shape[0]-5)
+ pos_ind = (-5 < x) & (x < self.ref.shape[1]+5) & (-5 < y) & (y < self.ref.shape[0]+5)
d = d.iloc[pos_ind]
+ x = x[pos_ind]; y = y[pos_ind]
# account for crowding
for i in range(len(d)):
- x = d.col.values[i]
- y = d.row.values[i]
+ xx = d.col.values[i]
+ yy = d.row.values[i]
- dist = np.sqrt((tab.col.values-x)**2 + (tab.row.values-y)**2)
+ dist = np.sqrt((tab['col'].values-xx)**2 + (tab['row'].values-yy)**2)
- ind = dist < 1.5
+ ind = dist < 1
close = tab.iloc[ind]
d['gmag'].iloc[i] = -2.5*np.log10(np.nansum(mag2flux(close.gmag.values,25))) + 25
@@ -2684,6 +2914,7 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None):
d['zmag'].iloc[i] = -2.5*np.log10(np.nansum(mag2flux(close.zmag.values,25))) + 25
if system == 'ps1':
d['ymag'].iloc[i] = -2.5*np.log10(np.nansum(mag2flux(close.ymag.values,25))) + 25
+
# convert to tess mags
if len(d) < 10:
print('!!!WARNING!!! field calibration is unreliable, using the default zp = 20.44')
@@ -2699,20 +2930,56 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None):
d = SM_to_TESS_mag(d,ebv=self.ebv)
+ maglim = 16.5
+ maglim_bright = 10
+ dist = np.sqrt((x[:,np.newaxis] - x[np.newaxis,:])**2 + (y[:,np.newaxis] - y[np.newaxis,:])**2)
+ mag_diff = d['tmag'].values[:,np.newaxis] - d['tmag'].values[np.newaxis,:]
+ mag_diff[mag_diff == 0] = np.nan
+ mag_diff[mag_diff < -2] = np.nan
+ mag_diff = np.isnan(mag_diff)
+ dist[dist==0] = np.nan
+ dist[mag_diff] = np.nan
+ min_dist = np.nanmin(dist,axis=1)
+ ind = min_dist > 5
+ if sum(ind) < 10:
+ ind = min_dist > 3
+ d = d.iloc[ind]
+ ind2 = (d['tmag'].values < maglim) & (d['tmag'].values > maglim_bright)
+ d = d.iloc[ind2]
+ xx = x[ind][ind2]; yy = y[ind][ind2]
+ ind = (xx > 5) & (xx < self.ref.shape[1]-5) & (yy > 5) & (yy < self.ref.shape[0]-5)
+ d = d.iloc[ind]
+ xx = xx[ind]; yy = yy[ind]
+ d['col'] = xx; d['row'] = yy
+
+ #self.cat['cal'] = 0
+ #self.cat['cal'].iloc[ind] = 1
+
+ if len(d) == 0:
+ print('!!! No suitable calibration sources !!!\nSetting to the default value of zp=20.44')
+ self.zp = 20.44
+ self.zp_e = 0.05
+ # backup for when messing around with flux later
+ self.tzp = 20.44
+ self.tzp_e = 0.05
+ return
+
flux = []
eflux = []
eind = np.zeros(len(d))
for i in range(len(d)):
mask = np.zeros_like(self.ref)
- xx = int(d.col.values[i] + .5); yy = int(d.row.values[i] + .5)
+ xx = int(np.round(d['col'].values[i],0)); yy = int(np.round(d['row'].values[i]))
mask[yy,xx] = 1
mask = convolve(mask,np.ones((3,3)))
if self.phot_method == 'aperture':
flux += [np.nansum(tflux*mask,axis=(1,2))]
elif self.phot_method == 'psf':
- flux += [self.psf_photometry(xPix=xx,yPix=yy,snap='ref',diff=False)]
+ f, e = self.psf_photometry(xPix=xx,yPix=yy,snap='ref',diff=False)
+ flux += [f]
+ eflux += [e]
m2 = np.zeros_like(self.ref)
- m2[int(d.row.values[i] + .5),int(d.col.values[i] + .5)] = 1
+ m2[int(np.round(d['row'].values[i],0)),int(np.round(d['col'].values[i]))] = 1
m2 = convolve(m2,np.ones((7,7))) - convolve(m2,np.ones((5,5)))
eflux += [np.nansum(tflux*m2,axis=(1,2))]
mag = -2.5*np.log10(np.nansum((ref*m2))) + 20.44
@@ -2727,10 +2994,9 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None):
#eind = abs(eflux) > 20
if self.phot_method == 'aperture':
flux[~eind] = np.nan
-
#calculate the zeropoint
- zp = d.tmag.values[:,np.newaxis] + 2.5*np.log10(flux)
+ zp = d['tmag'].values[:,np.newaxis] + 2.5*np.log10(flux)
if len(zp) == 0:
zp = np.array([20.44])
@@ -2750,12 +3016,14 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None):
if plot:
plt.figure()
nonan = np.isfinite(self.ref)
- plt.imshow(ref,origin='lower',vmax = np.percentile(ref[nonan],80),vmin=np.percentile(ref[nonan],10))
- plt.scatter(d.col.iloc[eind],d.row.iloc[eind],color='r')
+ plt.imshow(ref,origin='lower',vmax = np.percentile(ref[nonan],90),vmin=np.percentile(ref[nonan],10))
+ cbar = plt.colorbar()
+ cbar.ax.set_ylabel('Counts',fontsize=12)
+ plt.plot(d.col.iloc[eind],d.row.iloc[eind],'rx')
plt.title('Calibration sources')
plt.ylabel('Row',fontsize=15)
plt.xlabel('Column',fontsize=15)
- plt.colorbar()
+
plt.show()
if savename is not None:
plt.savefig(savename + 'cal_sources.pdf', bbox_inches = "tight")
@@ -2811,10 +3079,10 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None):
if compare:
print('!!!WARNING!!! field calibration is unreliable, using the default zp = 20.44')
self.zp = 20.44
- self.zp_e = 0.5
+ self.zp_e = 0.05
# backup for when messing around with flux later
self.tzp = 20.44
- self.tzp_e = 0.5
+ self.tzp_e = 0.05
else:
self.zp = mzp
self.zp_e = stdzp
diff --git a/tests/test_tessreduce_main.py b/tests/test_tessreduce_main.py
index 515679e..6225542 100755
--- a/tests/test_tessreduce_main.py
+++ b/tests/test_tessreduce_main.py
@@ -27,8 +27,8 @@ def test_align(self):
#def test_Shift_images(self):
# self.tess.shift_images()
- def test_field_calibrate(self):
- self.tess.field_calibrate()
+ #def test_field_calibrate(self):
+ # self.tess.field_calibrate()
def test_Diff_lc(self):
self.tess.diff_lc()