From bb67362ed2e6f659fd087f1f3c3b06d88b694dd5 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Sat, 31 May 2025 16:59:10 +1200 Subject: [PATCH 01/31] current working --- tessreduce/cat_mask.py | 13 ++-- tessreduce/delta_function_fitting.py | 14 ++-- tessreduce/helpers.py | 34 +++++++-- tessreduce/psf_photom.py | 4 +- tessreduce/tessreduce.py | 103 +++++++++++++++++---------- 5 files changed, 115 insertions(+), 53 deletions(-) diff --git a/tessreduce/cat_mask.py b/tessreduce/cat_mask.py index e392658..2e56531 100755 --- a/tessreduce/cat_mask.py +++ b/tessreduce/cat_mask.py @@ -105,10 +105,15 @@ def gaia_auto_mask(table,Image,scale=1): mags = [[18,17],[17,16],[16,15],[15,14],[14,13.5],[13.5,12],[12,10],[10,9],[9,8],[8,7]] size = (np.array([3,4,5,6,7,8,10,14,16,18])*scale).astype(int) for i, mag in enumerate(mags): - m = ((magim > mag[1]) & (magim <= mag[0])) * 1. - k = np.ones((size[i],size[i])) - conv = fftconvolve(m, k,mode='same')#.astype(int) - masks[str(mag[0])] = (conv >.1) * 1. + ind = (m > mag[1]) & (m <= mag[0]) + magim = np.zeros_like(image) + magim[y[ind],x[ind]] = 1. + if size[i] >0: + k = np.ones((size[i],size[i])) + conv = fftconvolve(magim, k,mode='same')#.astype(int) + masks[str(mag[0])] = (conv >.1) * 1. + else: + conv = magim masks['all'] = np.zeros_like(image,dtype=float) for key in masks: masks['all'] += masks[key] diff --git a/tessreduce/delta_function_fitting.py b/tessreduce/delta_function_fitting.py index fa16a79..c5518d9 100644 --- a/tessreduce/delta_function_fitting.py +++ b/tessreduce/delta_function_fitting.py @@ -5,8 +5,8 @@ from scipy.optimize import minimize -def Delta_basis(Size = 11): - kernel = np.zeros((Size,Size)) +def Delta_basis(size = 11): + kernel = np.zeros((size,size)) x,y = np.where(kernel==0) middle = int(len(x)/2) basis = [] @@ -22,17 +22,17 @@ def Delta_basis(Size = 11): coeff = np.ones(len(basis)) return basis, coeff -def Delta_kernel(reference,image,Size=11,mask=None): +def Delta_kernel(reference,image,size=11,mask=None): if mask is None: mask = np.ones_like(image) mask[mask == 0] = np.nan - Basis, coeff_0 = Delta_basis(Size) + Basis, coeff_0 = Delta_basis(size) bds = [] for i in range(len(coeff_0)): bds += [(0,1)] coeff_0 *= 0.01 - coeff_0[Size//2+1] = 0.95 - res = minimize(optimize_delta, coeff_0, args=(Basis,reference,image,Size,mask), + coeff_0[size//2+1] = 0.95 + res = minimize(optimize_delta, coeff_0, args=(Basis,reference,image,size,mask), bounds=bds,method='Powell') k = np.nansum(res.x[:,np.newaxis,np.newaxis]*Basis,axis=0) return k @@ -50,7 +50,7 @@ def optimize_delta(Coeff, Basis, reference, image,size,mask): def parallel_delta_diff(image,reference,mask=None,size=11): - kernel = Delta_kernel(reference,image,Size=11,mask=mask) + kernel = Delta_kernel(reference,image,size=size,mask=mask) template = signal.fftconvolve(reference, kernel, mode='same') diff = image - template return diff, kernel diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 8c99555..f6289bf 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -221,7 +221,7 @@ def Smooth_bkg(data, gauss_smooth=2, interpolate=False, extrapolate=True): # end inpaint estimate = inpaint.inpaint_biharmonic(data,mask) #estimate = signal.fftconvolve(estimate,self.prf,mode='same') - if np.nanmedian(estimate) < 20: + if (np.nanmedian(estimate) < 100) & (np.nanstd(estimate) < 3): # magic numbers to define a well behaved background gauss_smooth = gauss_smooth * 3 estimate = gaussian_filter(estimate,gauss_smooth) else: @@ -282,7 +282,7 @@ def Calculate_shifts(data,mx,my,finder): def image_sub(theta, image, ref): dx, dy = theta - s = shift(image,([dx,dy]),order=5) + s = shift(image,([dx,dy]),order=5, mode='nearest') #translation = np.float64([[1,0,dx],[0,1, dy]]) #s = cv2.warpAffine(image, translation, image.shape[::-1], flags=cv2.INTER_CUBIC,borderValue=0) diff = (ref-s)**2 @@ -291,6 +291,31 @@ def image_sub(theta, image, ref): else: return np.nansum(diff[5:-6,5:-6]) + +### TESTING CODE +#def c_shift_image(img, shift): +# return scipy.ndimage.shift(img, shift=shift, order=5, mode='nearest') # cubic interpolation +# +#def cost_function(shift, ref_img, moving_img): +# shifted = shift_image(moving_img, shift) +# diff = ref_img - shifted +# return np.sum(diff[2:-2,2:-2]**2) # Sum of Squared Differences (SSD) +# +#def align_subpixel(ref_img, moving_img, initial_shift=(0,0)): +# ref_img = ref_img.astype(np.float32) +# moving_img = moving_img.astype(np.float32) +# +# # Minimize SSD by shifting moving_img +# result = minimize(cost_function, initial_shift, args=(ref_img, moving_img), method='Powell') +# +# best_shift = result.x +# aligned_img = c_shift_image(moving_img, best_shift) +# +# print(f"Optimal shift (y, x): {best_shift}") +# return aligned_img, best_shift +### + + def difference_shifts(image,ref): """ Calculate the offsets of sources identified by photutils from a reference @@ -317,9 +342,10 @@ def difference_shifts(image,ref): """ if np.nansum(abs(image)) > 0: x0= [0,0] - bds = [(-1,1),(-1,1)] + bds = [(-1.5,1.5),(-1.5,1.5)] res = minimize(image_sub,x0,args=(image,ref),method = 'Powell',bounds= bds) s = res.x + #a,s = align_subpixel(ref,image) else: s = np.zeros((2)) * np.nan if (s == np.ones((2))).any(): @@ -346,7 +372,7 @@ def Smooth_motion(Centroids,tpf): """ smoothed = np.zeros_like(Centroids) * np.nan - skernel = int(len(tpf.flux) * 0.2) #simple way of making the smoothing window 10% of the duration + skernel = int(len(tpf.flux) * 0.01) #simple way of making the smoothing window 10% of the duration skernel = skernel // 2 +1 print('!!! skernel '+ str(skernel)) #skernel = 25 diff --git a/tessreduce/psf_photom.py b/tessreduce/psf_photom.py index f344513..fe03a02 100644 --- a/tessreduce/psf_photom.py +++ b/tessreduce/psf_photom.py @@ -126,7 +126,7 @@ def minimize_position(self,coeff,image,ext_shift): residual = np.nansum(diff**2) 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,limx=0.8,limy=0.8,ext_shift=[0,0]): """ Finds the optimal psf fit @@ -157,7 +157,7 @@ def psf_position(self,image,limx=0.5,limy=0.5,ext_shift=[0,0]): # -- Optimize -- # res = minimize(self.minimize_position, coeff, args=(normimage,ext_shift), method='Powell',bounds=lims) - print(res.x) + print('Optimal PSF shift: ', res.x) self.psf_fit = res def minimize_psf_flux(self,coeff,image,surface=True,order=2,kernel=None): diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 4d40336..1a7dd0d 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -223,6 +223,7 @@ 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) + self.flux[np.isnan(self.flux)] = 0 self.wcs = self.tpf.wcs self.ra = self.tpf.ra self.dec = self.tpf.dec @@ -390,6 +391,7 @@ 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 self.wcs = tpf.wcs def make_mask(self,catalogue_path=None,maglim=19,scale=1,strapsize=6,useref=False): @@ -522,7 +524,28 @@ 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): + ''' + Calculate the effective quantum efficiency enhancement of the detector from scattered light. + ''' + norm = (flux + bkg_smth) / 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 + qe_1d = np.nanpercentile(ratio,10,axis=1) + qe = np.ones_like(norm) + qe[:,:,:] = qe_1d[:,np.newaxis,:] + 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 +592,24 @@ 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: + 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 strap_iso: + over_sub[strap_mask] = 0 + if source_hunt: + 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 +620,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.qe = qe + bkg = bkg_smth * qe else: bkg = np.array(bkg_smth) self.bkg = bkg @@ -885,6 +911,7 @@ def fit_shift(self,smooth=True,plot=None,savename=None): shifts = np.zeros((len(f),2)) * np.nan for i in range(len(f)): shifts[i,:] = difference_shifts(f[i],m) + sraw = deepcopy(shifts) shifts = Smooth_motion(shifts,self.tpf) @@ -898,8 +925,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 +974,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 +1100,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: @@ -1547,13 +1576,14 @@ 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-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] + if isinstance(loc[0], (float, np.floating, np.float32, np.float64)): loc[0] = int(loc[0] + 0.5) 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 if ref: @@ -1855,7 +1885,7 @@ 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: @@ -1867,7 +1897,7 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, print('calculating background') # calculate the background - self.background() + self.background(rerun_negative=True) if np.isnan(self.bkg).all(): @@ -1927,7 +1957,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 +1966,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,8 +1975,9 @@ 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=18,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: @@ -2532,12 +2563,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['gmag'].iloc[i] = -2.5*np.log10(np.nansum(maselfflux(close.gmag.values,25))) + 25 + d['rmag'].iloc[i] = -2.5*np.log10(np.nansum(maselfflux(close.rmag.values,25))) + 25 + d['imag'].iloc[i] = -2.5*np.log10(np.nansum(maselfflux(close.imag.values,25))) + 25 + d['zmag'].iloc[i] = -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['ymag'].iloc[i] = -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') From ef8139fb4ede684515d0cfc29eee4687e94483ea Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Wed, 25 Jun 2025 09:17:05 +1200 Subject: [PATCH 02/31] trying to fix errors --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index fa868c2..bce1216 100755 --- a/setup.py +++ b/setup.py @@ -24,7 +24,7 @@ # What packages are required for this module to be executed? REQUIRED = ['lightkurve>=2.0.0', - 'numpy', + 'numpy<2.0.0', 'photutils>=1.4', 'pandas', 'scipy!=1.4.0,!=1.4.1,>=0.19.0', From ae6130e51abdec723ea6379994a070f3a1f72621 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Wed, 2 Jul 2025 10:52:46 +1200 Subject: [PATCH 03/31] added option to remove lightkurve cache files from tesscut --- tessreduce/calibration_tools.py | 4 ++-- tessreduce/tessreduce.py | 32 +++++++++++++++++++++++--------- 2 files changed, 25 insertions(+), 11 deletions(-) diff --git a/tessreduce/calibration_tools.py b/tessreduce/calibration_tools.py index 0beebc3..0e14f12 100755 --- a/tessreduce/calibration_tools.py +++ b/tessreduce/calibration_tools.py @@ -192,8 +192,8 @@ def Tonry_reduce(Data,plot=False,savename=None,system='ps1'): colours = Make_colours(dat.iloc[ind],tonry,compare,Extinction = res.x, Tonry = True,system=system) clip = Tonry_clip(colours,tonry) #clips += [clip] - dat['locus'].iloc[ind] = clip - dat2 = dat.iloc[dat['locus'].values > 0] + dat.loc[ind,'locus'] = clip * 1 + dat2 = dat.loc[dat['locus'] > 0] #print('Pass ' + str(i+1) + ': ' + str(res.x[0])) #clips[0][clips[0]] = clips[1] if plot: diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 1a7dd0d..e740954 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -65,7 +65,7 @@ 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, + savename=None,quality_bitmask='default',cache_dir=None,cache=True,catalogue_path=False, shift_method='difference',prf_path=None,verbose=1): """ @@ -235,8 +235,8 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect # 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) @@ -326,7 +326,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. @@ -383,6 +389,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: @@ -2563,12 +2577,12 @@ def isolated_star_lcs(self): ind = dist < 1.5 close = tab.iloc[ind] - d['gmag'].iloc[i] = -2.5*np.log10(np.nansum(maselfflux(close.gmag.values,25))) + 25 - d['rmag'].iloc[i] = -2.5*np.log10(np.nansum(maselfflux(close.rmag.values,25))) + 25 - d['imag'].iloc[i] = -2.5*np.log10(np.nansum(maselfflux(close.imag.values,25))) + 25 - d['zmag'].iloc[i] = -2.5*np.log10(np.nansum(maselfflux(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(maselfflux(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') From 20c8f87336fd624ed8b541732b3a3a8b91a01e62 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Wed, 2 Jul 2025 10:59:04 +1200 Subject: [PATCH 04/31] disabling a broken test --- tests/test_tessreduce_main.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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() From c226db22414de4c20f16bf1f04694fae3d0d08ff Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 7 Aug 2025 08:06:47 +1200 Subject: [PATCH 05/31] messing with errors --- tessreduce/psf_photom.py | 35 +++++++++++++++++++++++++++-------- 1 file changed, 27 insertions(+), 8 deletions(-) diff --git a/tessreduce/psf_photom.py b/tessreduce/psf_photom.py index fe03a02..79b9512 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,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] @@ -126,7 +134,7 @@ def minimize_position(self,coeff,image,ext_shift): residual = np.nansum(diff**2) return residual#np.exp(residual) - def psf_position(self,image,limx=0.8,limy=0.8,ext_shift=[0,0]): + def psf_position(self,image,limx=0.8,limy=0.8,ext_shift=[0,0],surface=False,order=2): """ Finds the optimal psf fit @@ -150,10 +158,18 @@ def psf_position(self,image,limx=0.8,limy=0.8,ext_shift=[0,0]): Optimal psf fit to input image """ - + #brightloc = normimage = image / np.nansum(image) # normalise the image - coeff = [self.source_x,self.source_y] - lims = [[-limx,limx],[-limy,limy]] + 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,ext_shift), method='Powell',bounds=lims) @@ -240,8 +256,11 @@ 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') - error = np.sqrt(np.diag(res['hess_inv'])) + #res = minimize(self.minimize_psf_flux,initial,args=(image,surface,poly_order,kernel),method='BFGS') + res = minimize(self.minimize_psf_flux,initial,args=(image,surface,poly_order,kernel),method='Powell') + error = 1 + #error = np.sqrt(np.diag(res['hess_inv'])) + self.res = res self.flux = res.x[0] self.eflux = error[0] From 90bea126e8545ac2b6e7225017cc289651ccc458 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 7 Aug 2025 08:11:02 +1200 Subject: [PATCH 06/31] messing with errors --- tessreduce/psf_photom.py | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/tessreduce/psf_photom.py b/tessreduce/psf_photom.py index 79b9512..ccbc0b9 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,surface=True,order=2): + def minimize_position(self,coeff,image,ext_shift):#,surface=True,order=2): """ Applies an exponential function using psf residuals to optimise the psf fit. @@ -114,14 +114,14 @@ def minimize_position(self,coeff,image,ext_shift,surface=True,order=2): 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 + # 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] @@ -134,7 +134,7 @@ def minimize_position(self,coeff,image,ext_shift,surface=True,order=2): residual = np.nansum(diff**2) return residual#np.exp(residual) - def psf_position(self,image,limx=0.8,limy=0.8,ext_shift=[0,0],surface=False,order=2): + def psf_position(self,image,limx=0.8,limy=0.8,ext_shift=[0,0]):#,surface=False,order=2): """ Finds the optimal psf fit @@ -160,16 +160,16 @@ def psf_position(self,image,limx=0.8,limy=0.8,ext_shift=[0,0],surface=False,orde """ #brightloc = 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]] + # 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,ext_shift), method='Powell',bounds=lims) From 71edda36d047b0457847b447255913bc5ad308b9 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 7 Aug 2025 08:12:08 +1200 Subject: [PATCH 07/31] messing with errors --- tessreduce/psf_photom.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tessreduce/psf_photom.py b/tessreduce/psf_photom.py index ccbc0b9..0c6f213 100644 --- a/tessreduce/psf_photom.py +++ b/tessreduce/psf_photom.py @@ -258,12 +258,11 @@ def psf_flux(self,image,ext_shift=None,surface=True,poly_order=3,kernel=None): #res = minimize(self.minimize_psf_flux,initial,args=(image,surface,poly_order,kernel),method='BFGS') res = minimize(self.minimize_psf_flux,initial,args=(image,surface,poly_order,kernel),method='Powell') - error = 1 #error = np.sqrt(np.diag(res['hess_inv'])) self.res = res self.flux = res.x[0] - self.eflux = error[0] + self.eflux = 1#error[0] if surface: x = np.arange(image.shape[1]) From c8a426cd224929599ec90d40dc775f7223cdae6b Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 7 Aug 2025 18:39:25 +1200 Subject: [PATCH 08/31] added in errors from the error array --- tessreduce/helpers.py | 19 +++--- tessreduce/psf_photom.py | 67 ++++++++++-------- tessreduce/tessreduce.py | 144 +++++++++++++++++++++++++-------------- 3 files changed, 139 insertions(+), 91 deletions(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index f6289bf..d863a52 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -280,12 +280,12 @@ def Calculate_shifts(data,mx,my,finder): shifts[1,:] = np.nan return shifts -def image_sub(theta, image, ref): +def image_sub(theta, image, ref, eimage, eref): dx, dy = theta s = shift(image,([dx,dy]),order=5, mode='nearest') #translation = np.float64([[1,0,dx],[0,1, dy]]) #s = cv2.warpAffine(image, translation, image.shape[::-1], flags=cv2.INTER_CUBIC,borderValue=0) - diff = (ref-s)**2 + diff = (ref-s)**2#/(eimage + eref) if image.shape[0] > 50: return np.nansum(diff[10:-11,10:-11]) else: @@ -316,7 +316,7 @@ def image_sub(theta, image, ref): ### -def difference_shifts(image,ref): +def difference_shifts(image,ref,eimage,eref): """ Calculate the offsets of sources identified by photutils from a reference @@ -343,7 +343,7 @@ def difference_shifts(image,ref): if np.nansum(abs(image)) > 0: x0= [0,0] bds = [(-1.5,1.5),(-1.5,1.5)] - res = minimize(image_sub,x0,args=(image,ref),method = 'Powell',bounds= bds) + res = minimize(image_sub,x0,args=(image,ref,eimage,eref),method = 'Powell',bounds= bds) s = res.x #a,s = align_subpixel(ref,image) else: @@ -374,7 +374,6 @@ def Smooth_motion(Centroids,tpf): smoothed = np.zeros_like(Centroids) * np.nan skernel = int(len(tpf.flux) * 0.01) #simple way of making the smoothing window 10% of the duration skernel = skernel // 2 +1 - print('!!! skernel '+ str(skernel)) #skernel = 25 if skernel < 25: skernel = 25 @@ -739,17 +738,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 diff --git a/tessreduce/psf_photom.py b/tessreduce/psf_photom.py index 0c6f213..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):#,surface=True,order=2): + 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. @@ -130,11 +130,11 @@ def minimize_position(self,coeff,image,ext_shift):#,surface=True,order=2): 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.8,limy=0.8,ext_shift=[0,0]):#,surface=False,order=2): + 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 @@ -158,25 +158,33 @@ def psf_position(self,image,limx=0.8,limy=0.8,ext_shift=[0,0]):#,surface=False,o Optimal psf fit to input image """ + if error is None: + error = np.ones_like(image) #brightloc = - 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,ext_shift), method='Powell',bounds=lims) - print('Optimal PSF shift: ', res.x) - self.psf_fit = res + 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 - 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): """ @@ -210,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): """ @@ -238,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 @@ -256,13 +265,13 @@ 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,surface,poly_order,kernel),method='Powell') - #error = np.sqrt(np.diag(res['hess_inv'])) + 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 = 1#error[0] + self.eflux = error[0] if surface: x = np.arange(image.shape[1]) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index e740954..6b119c5 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -223,6 +223,7 @@ 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) + self.eflux = strip_units(self.tpf.flux_err) self.flux[np.isnan(self.flux)] = 0 self.wcs = self.tpf.wcs self.ra = self.tpf.ra @@ -230,7 +231,7 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect 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(): @@ -406,6 +407,7 @@ def get_TESS(self,ra=None,dec=None,name=None,size=None,sector=None, self.tpf = tpf self.flux = strip_units(tpf.flux) # Stripping astropy units so only numbers are returned self.flux[np.isnan(self.flux)] = 0 + self.eflux = strip_units(tpf.flux_err) self.wcs = tpf.wcs def make_mask(self,catalogue_path=None,maglim=19,scale=1,strapsize=6,useref=False): @@ -551,12 +553,13 @@ def _calc_qe(self,flux,bkg_smth): m = (np.nansum(masks,axis=0) > 0) * 1. m[m==0] = np.nan ratio = flux / bkg_smth * m - qe_1d = np.nanpercentile(ratio,10,axis=1) + 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[:,:,:] = qe_1d[:,np.newaxis,:] + qe[:,:,:] = med[:,np.newaxis,:] 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): @@ -607,11 +610,11 @@ def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False 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: - over_sub = (deepcopy(self.flux) - bkg_smth) < -0.5 + over_sub = (deepcopy(self.flux) - bkg_smth) < -3 * self.eflux over_sub = np.nansum(over_sub,axis=0) > 0 self.over_sub = over_sub - print('overshape ',over_sub.shape) - print('m ',m.shape) + #print('overshape ',over_sub.shape) + #print('m ',m.shape) strap_mask = (self.mask & 4) > 0 if strap_iso: over_sub[strap_mask] = 0 @@ -916,15 +919,18 @@ 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) 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) + shifts[i,:] = difference_shifts(f[i],m,self.eflux[i],eref) sraw = deepcopy(shifts) shifts = Smooth_motion(shifts,self.tpf) @@ -1194,12 +1200,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]) @@ -1229,7 +1235,7 @@ 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 + tar_err = np.nansum((self.eflux)*ap_tar,axis=(1,2))#sky_std * tar_ap**2 if phot_method == 'psf': if psf_snap is None: psf_snap = 'brightest' @@ -1595,16 +1601,17 @@ def _psf_initialise(self,cutoutSize,loc,time_ind=None,ref=False): row = self.tpf.row - int(self.size/2-1) + 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) + loc[1] = int(np.round(loc[1])) prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row) # initialise psf kernel 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 + 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 + return prf, cutout, ecutout def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): """ @@ -1641,14 +1648,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) @@ -1713,7 +1720,7 @@ 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)) @@ -1738,7 +1745,7 @@ 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 @@ -1747,26 +1754,26 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa #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) @@ -1779,12 +1786,12 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa 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) @@ -2604,11 +2611,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 @@ -2663,7 +2670,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: @@ -2696,31 +2702,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 @@ -2729,6 +2736,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') @@ -2744,20 +2752,53 @@ 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 + + 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 @@ -2772,10 +2813,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]) @@ -2856,10 +2896,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 From c236bebc9cbf8dea8e4b772ad60150ac322de242 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Wed, 13 Aug 2025 08:30:21 +1200 Subject: [PATCH 09/31] added in new functions for the qe correction. Now will capture variation in the QE along the strap --- tessreduce/helpers.py | 87 ++++++++++++++++++++++++++++++++++++++++ tessreduce/tessreduce.py | 75 ++++++++++++++++++++-------------- 2 files changed, 132 insertions(+), 30 deletions(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index d863a52..7f5d4dd 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -1232,3 +1232,90 @@ 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] + q = savgol_filter(q,61,2) + #sub = ((frame - (q * frame_bkg)))[:,i] + #ind = (sub < tol * frame_err)[:,i] + #nm[ind] = True + qe[:,i] = q + #qe[(qe-1) < 5e-3] = 1 + return qe + + + + + + + + + diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 6b119c5..793adf6 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 @@ -280,9 +280,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) @@ -518,7 +518,7 @@ def psf_source_mask(self,sigma=5): 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) + self.tpf.column+self.flux.shape[2]/2,self.tpf.row+self.flux.shape[1]/2) self.prf = prf.locate(5,5,(11,11)) @@ -540,24 +540,39 @@ def psf_source_mask(self,sigma=5): m[i] = par_psf_source_mask(data[i],self.prf,sigma) return m * 1.0 - def _calc_qe(self,flux,bkg_smth): + 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) / bkg_smth - straps = norm * ((self.mask & 4)>0) + 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 - 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.ones_like(flux) * 1. + if flux.shape[1] < 30: + 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 @@ -638,7 +653,7 @@ def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False # Calculate quantum efficiency if calc_qe: self.bkg = bkg_smth - qe = self._calc_qe(flux,bkg_smth) + qe = self._calc_qe(flux,bkg_smth,self.eflux) self.qe = qe bkg = bkg_smth * qe else: @@ -700,7 +715,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) @@ -725,7 +740,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)): @@ -751,7 +766,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)): @@ -841,7 +856,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) @@ -1316,9 +1331,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) @@ -1649,7 +1664,7 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): inds = np.arange(0,len(xpos)) if self.parallel: 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, + self.tpf.sector,self.tpf.column,self.tpf.row, size,[xpos[i],ypos[i]],time_ind) for i in inds)) else: prfs = [] @@ -1720,12 +1735,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, ecutouts = 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): @@ -1745,7 +1760,7 @@ 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, ecutouts = 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 @@ -1756,14 +1771,14 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa base = create_psf(prf,size) base.psf_position(ref,ecutouts[ind],ext_shift=self.shift[ind]) elif snap == 'ref': - prf, cutouts,ecutouts = 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,ecutouts[self.ref_ind]) if diff: _, cutouts,ecutouts = self._psf_initialise(size,(xPix,yPix),ref=False) elif snap == 'fixed': - prf, cutouts, ecutouts = 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)) @@ -1780,7 +1795,7 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa 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: @@ -2104,7 +2119,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 From 3ffef8892cd7148960a8170893673356698b5297 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Wed, 13 Aug 2025 08:31:03 +1200 Subject: [PATCH 10/31] small change --- tessreduce/helpers.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 7f5d4dd..64b540b 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -1303,10 +1303,12 @@ def parallel_strap_fit(frame,frame_bkg,frame_err,mask,repeats=3,tol=3): #for r in range(repeats): q = fit_strap(norm[:,i],nm) #q /= frame_bkg[:,i] - q = savgol_filter(q,61,2) - #sub = ((frame - (q * frame_bkg)))[:,i] - #ind = (sub < tol * frame_err)[:,i] - #nm[ind] = True + if y.shape > 70: + q = savgol_filter(q,61,2) + else: + mq = np.nanmed(q) + q[:] = mq + qe[:,i] = q #qe[(qe-1) < 5e-3] = 1 return qe From f546db0168f117b1958146d17be2f618f6463b99 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Wed, 13 Aug 2025 16:58:23 +1200 Subject: [PATCH 11/31] fix --- tessreduce/helpers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 64b540b..8972b0c 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -1303,7 +1303,7 @@ def parallel_strap_fit(frame,frame_bkg,frame_err,mask,repeats=3,tol=3): #for r in range(repeats): q = fit_strap(norm[:,i],nm) #q /= frame_bkg[:,i] - if y.shape > 70: + if len(y) > 70: q = savgol_filter(q,61,2) else: mq = np.nanmed(q) From 0a9cb42ae5f2c452e4db6d62b1c91077906397e6 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 14 Aug 2025 20:56:26 +1200 Subject: [PATCH 12/31] some changes to straps and added in a new parameter for column offset --- tessreduce/cat_mask.py | 18 +++++----- tessreduce/helpers.py | 36 ++++++++++++++++--- tessreduce/tessreduce.py | 76 +++++++++++++++++++++++++++++++++++----- 3 files changed, 108 insertions(+), 22 deletions(-) diff --git a/tessreduce/cat_mask.py b/tessreduce/cat_mask.py index 2e56531..42e7748 100755 --- a/tessreduce/cat_mask.py +++ b/tessreduce/cat_mask.py @@ -25,7 +25,7 @@ def circle_app(rad): """ Makes a kinda circular aperture, probably not worth using. """ - mask = np.zeros((int(rad*2+.5)+1,int(rad*2+.5)+1)) + mask = np.zeros((int(np.round(rad*2,0))+1,int(np.round(rad*2))+1)) c = rad x,y =np.where(mask==0) dist = np.sqrt((x-c)**2 + (y-c)**2) @@ -42,8 +42,8 @@ def ps1_auto_mask(table,Image,scale=1): image = np.zeros_like(Image) x = table.x.values y = table.y.values - x = (x+.5).astype(int) - y = (y+.5).astype(int) + x = (np.round(x,0)).astype(int) + y = (np.round(y,0)).astype(int) m = table.mag.values ind = size_limit(x,y,image) x = x[ind]; y = y[ind]; m = m[ind] @@ -91,8 +91,8 @@ def gaia_auto_mask(table,Image,scale=1): image = np.zeros_like(Image) x = table.x.values y = table.y.values - x = (x+.5).astype(int) - y = (y+.5).astype(int) + x = (np.round(x,0)).astype(int) + y = (np.round(y,0)).astype(int) m = table.mag.values ind = size_limit(x,y,image) x = x[ind]; y = y[ind]; m = m[ind] @@ -148,8 +148,8 @@ def Big_sat(table,Image,scale=1): sat = table.iloc[i] x = sat.x.values y = sat.y.values - x = (x+.5).astype(int) - y = (y+.5).astype(int) + x = (np.round(x,0)).astype(int) + y = (np.round(y,0)).astype(int) m = sat.mag.values ind = size_limit(x,y,image) @@ -222,7 +222,7 @@ def Strap_mask(Image,col,size=4): big_strap = fftconvolve(strap_mask,np.ones((size,size)),mode='same') > .5 return big_strap -def Cat_mask(tpf,catalogue_path=None,maglim=19,scale=1,strapsize=3,ref=None,sigma=3): +def Cat_mask(tpf,catalogue_path=None,maglim=19,scale=1,strapsize=3,ref=None,sigma=3,col_offset=0): """ Make a source mask from the PS1 and Gaia catalogs. @@ -282,7 +282,7 @@ def Cat_mask(tpf,catalogue_path=None,maglim=19,scale=1,strapsize=3,ref=None,sigm sat = (np.nansum(sat,axis=0) > 0).astype(int) * 2 # assign 2 bit if strapsize > 0: - strap = Strap_mask(image,tpf.column,strapsize).astype(int) * 4 # assign 4 bit + strap = Strap_mask(image,tpf.column,strapsize+col_offset).astype(int) * 4 # assign 4 bit else: strap = np.zeros_like(image,dtype=int) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 8972b0c..67a84d5 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -1303,21 +1303,49 @@ def parallel_strap_fit(frame,frame_bkg,frame_err,mask,repeats=3,tol=3): #for r in range(repeats): q = fit_strap(norm[:,i],nm) #q /= frame_bkg[:,i] - if len(y) > 70: - q = savgol_filter(q,61,2) + if len(y) > 110: + q = savgol_filter(q,101,1) else: - mq = np.nanmed(q) + mq = np.nanmedian(q) q[:] = mq qe[:,i] = q - #qe[(qe-1) < 5e-3] = 1 + + #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): + if np.nansum(abs(cutout)) > 0: + phot = psf_phot(cutout, error=e_cutout,init_params=init_params) + phot = phot.to_pandas() + return phot[['flux_fit']].values[0], phot[['flux_err']].values[0] + #return float(phot[['flux_fit']].values), float(phot[['flux_err']].values) + else: + return np.array([np.nan]), np.array([np.nan]) + diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 793adf6..4021efa 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -66,7 +66,7 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect 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,cache=True,catalogue_path=False, - shift_method='difference',prf_path=None,verbose=1): + shift_method='difference',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: @@ -186,6 +187,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 @@ -447,9 +449,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. @@ -553,10 +555,10 @@ def _calc_qe(self,flux,bkg_smth,flux_e): m = (np.nansum(masks,axis=0) > 0) * 1. m[m==0] = np.nan qe = np.ones_like(flux) * 1. - if flux.shape[1] < 30: + 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) + m, med, std = sigma_clipped_stats(ratio,axis=1,sigma_upper=2,sigma_lower=3,maxiters=10) #qe_1d = np.nanpercentile(ratio,10,axis=1) #qe[:,:,:] = qe_1d[:,np.newaxis,:] @@ -1612,8 +1614,8 @@ def _psf_initialise(self,cutoutSize,loc,time_ind=None,ref=False): time_ind = np.arange(0,len(self.flux)) - 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] + 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(np.round(loc[0],0)) @@ -1689,6 +1691,61 @@ 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,yPix,size=5,local_bkg=False,ref=False,return_pos=False): + from photutils.psf import PSFPhotometry + from astropy.table import Table + rad = size // 2 + 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 + self.epsf = simulate_epsf(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row) + if local_bkg: + localbkg_estimator = LocalBackground(1.5, 7, bkgstat) + else: + localbkg_estimator = None + if ref: + flux = self.flux + self.ref + else: + flux = self.flux + cutouts = flux[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1] + ecutouts = self.eflux[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1] + + weight = np.abs(np.nansum((cutouts/ecutouts)[:,1:4,1:4],axis=(1,2))) + weight[np.isnan(weight)] = 0 + ind = np.argmax(weight) + fit_shape = (size, size) + # initial position fit on brightest frame + psfphot = PSFPhotometry(self.epsf, fit_shape, finder=None,aperture_radius=1.5, + localbkg_estimator=localbkg_estimator) + init = Table() + init['x_init'] = [size//2] + init['y_init'] = [size//2] + phot = psfphot(cutouts[ind], error=ecutouts[ind],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(self.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(flux)))) + 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 + 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 @@ -1763,9 +1820,10 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa 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 / ecuts)[:,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) From 5a127b822e79f87a3eb9113816dc2d1a7dfeda0a Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Fri, 15 Aug 2025 09:47:42 +1200 Subject: [PATCH 13/31] reverting qe correction to previous stable state --- tessreduce/tessreduce.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 4021efa..f0380d9 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -558,7 +558,7 @@ def _calc_qe(self,flux,bkg_smth,flux_e): 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,maxiters=10) + 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,:] From 44ab58dd5c0c2e9e65b022553040f98244b119d2 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Fri, 15 Aug 2025 13:14:09 +1200 Subject: [PATCH 14/31] fix --- tessreduce/cat_mask.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tessreduce/cat_mask.py b/tessreduce/cat_mask.py index 42e7748..7cd80d1 100755 --- a/tessreduce/cat_mask.py +++ b/tessreduce/cat_mask.py @@ -282,7 +282,7 @@ def Cat_mask(tpf,catalogue_path=None,maglim=19,scale=1,strapsize=3,ref=None,sigm sat = (np.nansum(sat,axis=0) > 0).astype(int) * 2 # assign 2 bit if strapsize > 0: - strap = Strap_mask(image,tpf.column,strapsize+col_offset).astype(int) * 4 # assign 4 bit + strap = Strap_mask(image,tpf.column+col_offset,strapsize).astype(int) * 4 # assign 4 bit else: strap = np.zeros_like(image,dtype=int) From c283a32e4622d77bf57504d32226cca281d498ef Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Fri, 15 Aug 2025 17:43:55 +1200 Subject: [PATCH 15/31] winding back qe changes... --- tessreduce/helpers.py | 136 ++++++++++++++++++----------------- tessreduce/tessreduce.py | 149 ++++++++++++++++++++++++--------------- 2 files changed, 167 insertions(+), 118 deletions(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 67a84d5..c8a73e8 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -53,18 +53,18 @@ def strip_units(data): """ - Removes the units off of data that was not in a NDarray, such as an astropy table. Returns an NDarray that has no units - - Parameters: - ---------- - data: ArrayLike - ArrayLike set of data that may have associated units that want to be removed. Should be able to return something sensible when .values is called. - - Returns: - ------- - data: ArrayLike - Same shape as input data, but will not have any units - """ + Removes the units off of data that was not in a NDarray, such as an astropy table. Returns an NDarray that has no units + + Parameters: + ---------- + data: ArrayLike + ArrayLike set of data that may have associated units that want to be removed. Should be able to return something sensible when .values is called. + + Returns: + ------- + data: ArrayLike + Same shape as input data, but will not have any units + """ if type(data) != np.ndarray: data = data.value @@ -294,25 +294,25 @@ def image_sub(theta, image, ref, eimage, eref): ### TESTING CODE #def c_shift_image(img, shift): -# return scipy.ndimage.shift(img, shift=shift, order=5, mode='nearest') # cubic interpolation +# return scipy.ndimage.shift(img, shift=shift, order=5, mode='nearest') # cubic interpolation # #def cost_function(shift, ref_img, moving_img): -# shifted = shift_image(moving_img, shift) -# diff = ref_img - shifted -# return np.sum(diff[2:-2,2:-2]**2) # Sum of Squared Differences (SSD) +# shifted = shift_image(moving_img, shift) +# diff = ref_img - shifted +# return np.sum(diff[2:-2,2:-2]**2) # Sum of Squared Differences (SSD) # #def align_subpixel(ref_img, moving_img, initial_shift=(0,0)): -# ref_img = ref_img.astype(np.float32) -# moving_img = moving_img.astype(np.float32) +# ref_img = ref_img.astype(np.float32) +# moving_img = moving_img.astype(np.float32) # -# # Minimize SSD by shifting moving_img -# result = minimize(cost_function, initial_shift, args=(ref_img, moving_img), method='Powell') +# # Minimize SSD by shifting moving_img +# result = minimize(cost_function, initial_shift, args=(ref_img, moving_img), method='Powell') # -# best_shift = result.x -# aligned_img = c_shift_image(moving_img, best_shift) +# best_shift = result.x +# aligned_img = c_shift_image(moving_img, best_shift) # -# print(f"Optimal shift (y, x): {best_shift}") -# return aligned_img, best_shift +# print(f"Optimal shift (y, x): {best_shift}") +# return aligned_img, best_shift ### @@ -450,36 +450,36 @@ def smooth_zp(zp,time): return smoothed, err def grads_rad(flux): - """ - Calculates the radius of the flux from the gradient of the flux, and the double gradient of the flux. - - Parameters: - ---------- - flux: ArrayLike - An array of flux values - - Returns: - ------- - rad: ArrayLike - The radius of the fluxes - """ - rad = np.sqrt(np.gradient(flux)**2+np.gradient(np.gradient(flux))**2) - return rad + """ + Calculates the radius of the flux from the gradient of the flux, and the double gradient of the flux. + + Parameters: + ---------- + flux: ArrayLike + An array of flux values + + Returns: + ------- + rad: ArrayLike + The radius of the fluxes + """ + rad = np.sqrt(np.gradient(flux)**2+np.gradient(np.gradient(flux))**2) + return rad def grad_flux_rad(flux): """ - Calculates the radius of the flux from the gradient of the flux. - - Parameters: - ---------- - flux: ArrayLike - An array of flux values - - Returns: - ------- - rad: ArrayLike - The radius of the fluxes - """ + Calculates the radius of the flux from the gradient of the flux. + + Parameters: + ---------- + flux: ArrayLike + An array of flux values + + Returns: + ------- + rad: ArrayLike + The radius of the fluxes + """ rad = np.sqrt(flux**2+np.gradient(flux)**2) return rad @@ -1249,7 +1249,7 @@ def grad_clip(data,box_size=100): """ gradind = np.zeros_like(data) - + for i in range(len(data)): if i < box_size//2: d = data[:i+box_size//2] @@ -1257,7 +1257,7 @@ def grad_clip(data,box_size=100): 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: @@ -1268,7 +1268,7 @@ def grad_clip(data,box_size=100): gradind[i-box_size//2:][ind] = gind else: gradind[i-box_size//2:i+box_size//2][ind] = gind - + gradind = gradind > 0 return gradind @@ -1337,15 +1337,25 @@ def simulate_epsf(camera,ccd,sector,column,row,size=13,realizations=2000,oversam return epsf -def parallel_photutils(cutout,e_cutout,psf_phot,init_params=None): - if np.nansum(abs(cutout)) > 0: - phot = psf_phot(cutout, error=e_cutout,init_params=init_params) - phot = phot.to_pandas() - return phot[['flux_fit']].values[0], phot[['flux_err']].values[0] - #return float(phot[['flux_fit']].values), float(phot[['flux_err']].values) - else: - return np.array([np.nan]), np.array([np.nan]) - +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/tessreduce.py b/tessreduce/tessreduce.py index f0380d9..87b961e 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -547,38 +547,53 @@ 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 = 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 - 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) - + ratio = flux / bkg_smth * m + #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 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. @@ -1691,52 +1706,76 @@ 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,yPix,size=5,local_bkg=False,ref=False,return_pos=False): + 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 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 - self.epsf = simulate_epsf(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row) + 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 + 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 = self.flux + self.ref - else: + if flux is None: flux = self.flux - cutouts = flux[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1] - ecutouts = self.eflux[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1] + if (flux.shape[1] < size) | (flux.shape[2] < size): + e = 'Image dimensions must be larger than the cutout size' + raise ValueError(e) + if (xPix is None) | (yPix is None): + xPix = flux.shape[2]//2 + yPix = flux.shape[1]//2 + if eflux is None: + eflux = self.eflux + if ref: + flux = flux + self.ref - weight = np.abs(np.nansum((cutouts/ecutouts)[:,1:4,1:4],axis=(1,2))) - weight[np.isnan(weight)] = 0 - ind = np.argmax(weight) + 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) - # initial position fit on brightest frame - psfphot = PSFPhotometry(self.epsf, fit_shape, finder=None,aperture_radius=1.5, + psfphot = PSFPhotometry(epsf, fit_shape, finder=None,aperture_radius=1.5, localbkg_estimator=localbkg_estimator) init = Table() init['x_init'] = [size//2] init['y_init'] = [size//2] - phot = psfphot(cutouts[ind], error=ecutouts[ind],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(self.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(flux)))) - 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 + 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: @@ -1821,7 +1860,7 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa 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 / ecuts)[:,int(yPix)-1:int(yPix)+2,int(xPix)-1:int(xPix)+2],axis=(1,2))) + 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(weight) #ind = np.where(cutouts==np.nanmax(cutouts))[0][0] From a6b9ca5a7a386068f90983f24a5bc004d77161cc Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Fri, 15 Aug 2025 17:46:13 +1200 Subject: [PATCH 16/31] fix --- tessreduce/tessreduce.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 87b961e..6fbadc2 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -1712,6 +1712,10 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, from photutils.psf import PSFPhotometry from astropy.table import Table rad = size // 2 + 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 @@ -1728,9 +1732,6 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, if (flux.shape[1] < size) | (flux.shape[2] < size): e = 'Image dimensions must be larger than the cutout size' raise ValueError(e) - if (xPix is None) | (yPix is None): - xPix = flux.shape[2]//2 - yPix = flux.shape[1]//2 if eflux is None: eflux = self.eflux if ref: From 1f50ae909eaaeeae20ec033a99122e05c84faf66 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Sat, 16 Aug 2025 11:16:17 +1200 Subject: [PATCH 17/31] attempt to make things work again --- tessreduce/tessreduce.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 6fbadc2..60b2d7d 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -555,6 +555,7 @@ def _calc_qe(self,flux,bkg_smth,flux_e): 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) @@ -642,7 +643,7 @@ def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False 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: - over_sub = (deepcopy(self.flux) - bkg_smth) < -3 * self.eflux + over_sub = (deepcopy(self.flux) - bkg_smth) < -self.eflux # -0.5 over_sub = np.nansum(over_sub,axis=0) > 0 self.over_sub = over_sub #print('overshape ',over_sub.shape) @@ -1712,10 +1713,19 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, 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 @@ -1727,13 +1737,7 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, localbkg_estimator = LocalBackground(1.5, 7, bkgstat) else: localbkg_estimator = None - 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 ref: flux = flux + self.ref From 22d18baf8c25137614ba4145c2cec671151b2109 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Sat, 16 Aug 2025 11:45:35 +1200 Subject: [PATCH 18/31] fix to moving mask --- tessreduce/tessreduce.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 60b2d7d..f53c405 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -172,6 +172,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 @@ -645,7 +646,7 @@ def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False if rerun_negative: over_sub = (deepcopy(self.flux) - bkg_smth) < -self.eflux # -0.5 over_sub = np.nansum(over_sub,axis=0) > 0 - self.over_sub = over_sub + #self._over_sub = over_sub #print('overshape ',over_sub.shape) #print('m ',m.shape) strap_mask = (self.mask & 4) > 0 @@ -2030,6 +2031,11 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, 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') From 09195f1da78ec45d7fa372f69da2c65a4560fd33 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Sun, 17 Aug 2025 15:32:29 +1200 Subject: [PATCH 19/31] small fix --- tessreduce/tessreduce.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index f53c405..6f56aec 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -564,6 +564,8 @@ def _calc_qe(self,flux,bkg_smth,flux_e): #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 @@ -650,6 +652,8 @@ def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False #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: From 88ef0213e747df0d4166c1d0207d3132403ee7f5 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Sun, 17 Aug 2025 16:01:44 +1200 Subject: [PATCH 20/31] small fix --- tessreduce/tessreduce.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 6f56aec..d3d4a04 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -656,7 +656,7 @@ def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False strap_mask = strap_mask[0] if strap_iso: over_sub[strap_mask] = 0 - if source_hunt: + if source_hunt | (len(m.shape == 3)): m[:,over_sub[:,:]] = 1 else: m[over_sub] = 1 From d5c0b1a6b3bdc67d1977f6c326f070fcc8631c64 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Sun, 17 Aug 2025 16:46:43 +1200 Subject: [PATCH 21/31] update before freeze --- tessreduce/tessreduce.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index d3d4a04..b15862b 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -656,7 +656,7 @@ def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False strap_mask = strap_mask[0] if strap_iso: over_sub[strap_mask] = 0 - if source_hunt | (len(m.shape == 3)): + if source_hunt | (len(self.mask.shape) == 3): m[:,over_sub[:,:]] = 1 else: m[over_sub] = 1 From 9881e6f34f7cbc41b890904d8b091875b1140a3b Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Wed, 20 Aug 2025 12:06:53 +1200 Subject: [PATCH 22/31] dissabled check_trend --- tessreduce/tessreduce.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index b15862b..63f1d20 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -1297,7 +1297,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): From 931f3cd2c8e55b42c7d59e6c3eb21b94ce8a0ded Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Wed, 20 Aug 2025 14:01:22 +1200 Subject: [PATCH 23/31] updated sector times --- development/update_tess_sectors.ipynb | 238 ++++++++++++++++++++++++++ tessreduce/sector_mjd.csv | 182 +++++++++++--------- 2 files changed, 336 insertions(+), 84 deletions(-) create mode 100644 development/update_tess_sectors.ipynb 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", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Sectormjd_startmjd_end
0158324.81597258352.684028
0258353.60763958381.020833
0358382.22569458408.875000
0458410.40625058436.357639
0558437.48263958463.795139
............
09360829.35763960842.638889
09460855.76388960868.236111
09560881.83333360894.173611
09660907.27430660919.989583
09760933.15972260946.371528
\n", + "

97 rows × 3 columns

\n", + "
" + ], + "text/plain": [ + " Sector mjd_start mjd_end\n", + "0 1 58324.815972 58352.684028\n", + "0 2 58353.607639 58381.020833\n", + "0 3 58382.225694 58408.875000\n", + "0 4 58410.406250 58436.357639\n", + "0 5 58437.482639 58463.795139\n", + ".. ... ... ...\n", + "0 93 60829.357639 60842.638889\n", + "0 94 60855.763889 60868.236111\n", + "0 95 60881.833333 60894.173611\n", + "0 96 60907.274306 60919.989583\n", + "0 97 60933.159722 60946.371528\n", + "\n", + "[97 rows x 3 columns]" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8efcc1f6", + "metadata": {}, + "outputs": [], + "source": [ + "def update_sector_times():\n", + " import pandas as pd \n", + " from astropy.time import Time\n", + " data = pd.read_csv('https://tess.mit.edu/public/files/TESS_orbit_times.csv',skipfooter=1)\n", + " sectors = data['Sector'].unique()\n", + " 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])" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} 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 From 275d80cec550f97bd51af509453e4ff6095d9ca1 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Wed, 20 Aug 2025 15:05:29 +1200 Subject: [PATCH 24/31] added an option to include the error image for memory concerns --- tessreduce/helpers.py | 7 ++++--- tessreduce/tessreduce.py | 37 +++++++++++++++++++++++++++---------- 2 files changed, 31 insertions(+), 13 deletions(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index c8a73e8..774ca57 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -280,7 +280,7 @@ def Calculate_shifts(data,mx,my,finder): shifts[1,:] = np.nan return shifts -def image_sub(theta, image, ref, eimage, eref): +def image_sub(theta, image, ref):#, eimage, eref): dx, dy = theta s = shift(image,([dx,dy]),order=5, mode='nearest') #translation = np.float64([[1,0,dx],[0,1, dy]]) @@ -316,7 +316,7 @@ def image_sub(theta, image, ref, eimage, eref): ### -def difference_shifts(image,ref,eimage,eref): +def difference_shifts(image,ref):#,eimage,eref): """ Calculate the offsets of sources identified by photutils from a reference @@ -343,7 +343,8 @@ def difference_shifts(image,ref,eimage,eref): if np.nansum(abs(image)) > 0: x0= [0,0] bds = [(-1.5,1.5),(-1.5,1.5)] - res = minimize(image_sub,x0,args=(image,ref,eimage,eref),method = 'Powell',bounds= bds) + #res = minimize(image_sub,x0,args=(image,ref,eimage,eref),method = 'Powell',bounds= bds) + res = minimize(image_sub,x0,args=(image,ref),method = 'Powell',bounds= bds) s = res.x #a,s = align_subpixel(ref,image) else: diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 63f1d20..0677af0 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -66,7 +66,7 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect 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,cache=True,catalogue_path=False, - shift_method='difference',prf_path=None,verbose=1,col_offset=0): + 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. @@ -154,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: @@ -226,7 +227,10 @@ 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) - self.eflux = strip_units(self.tpf.flux_err) + 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 @@ -410,7 +414,10 @@ def get_TESS(self,ra=None,dec=None,name=None,size=None,sector=None, self.tpf = tpf self.flux = strip_units(tpf.flux) # Stripping astropy units so only numbers are returned self.flux[np.isnan(self.flux)] = 0 - self.eflux = strip_units(tpf.flux_err) + 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): @@ -543,7 +550,7 @@ def psf_source_mask(self,sigma=5): m[i] = par_psf_source_mask(data[i],self.prf,sigma) return m * 1.0 - def _calc_qe(self,flux,bkg_smth,flux_e): + def _calc_qe(self,flux,bkg_smth):#,flux_e): ''' Calculate the effective quantum efficiency enhancement of the detector from scattered light. ''' @@ -646,7 +653,10 @@ def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False 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: - over_sub = (deepcopy(self.flux) - bkg_smth) < -self.eflux # -0.5 + 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) @@ -676,7 +686,7 @@ def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False # Calculate quantum efficiency if calc_qe: self.bkg = bkg_smth - qe = self._calc_qe(flux,bkg_smth,self.eflux) + qe = self._calc_qe(flux,bkg_smth)#,self.eflux) self.qe = qe bkg = bkg_smth * qe else: @@ -957,18 +967,20 @@ 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] + #eref = self.eflux[self.ref_ind] if self.parallel: ind = np.arange(len(f)) shifts = Parallel(n_jobs=self.num_cores)( - delayed(difference_shifts)(f[i],m,self.eflux[i],eref) for i in ind) + #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,self.eflux[i],eref) + shifts[i,:] = difference_shifts(f[i],m) sraw = deepcopy(shifts) shifts = Smooth_motion(shifts,self.tpf) @@ -1273,7 +1285,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 = np.nansum((self.eflux)*ap_tar,axis=(1,2))#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' @@ -1648,7 +1663,9 @@ def _psf_initialise(self,cutoutSize,loc,time_ind=None,ref=False): 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 + #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: return prf, cutout, ecutout def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): From 844c43e331ef7474fd58d5858db9f00c67d2b6a2 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 21 Aug 2025 09:51:23 +1200 Subject: [PATCH 25/31] fix --- tessreduce/tessreduce.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 0677af0..6937eac 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -1663,8 +1663,10 @@ def _psf_initialise(self,cutoutSize,loc,time_ind=None,ref=False): 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 - #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 + 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 From 0e4a6ac8312a0acd5607da2e26dc3443db7dce13 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 21 Aug 2025 11:07:56 +1200 Subject: [PATCH 26/31] scale of the differenced mask generation was not implemented correctly --- tessreduce/tessreduce.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 6937eac..afc4e9a 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -2144,12 +2144,12 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, 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*.5,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 From 26dcd83c68ff4e48c6a0ffb2a9c20cbb561ff150 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Sun, 14 Sep 2025 05:29:30 +1200 Subject: [PATCH 27/31] small shifts to the tessprf call --- tessreduce/tessreduce.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index afc4e9a..2fe2f76 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -1657,6 +1657,13 @@ def _psf_initialise(self,cutoutSize,loc,time_ind=None,ref=False): loc[0] = int(np.round(loc[0],0)) if isinstance(loc[1], (float, np.floating, np.float32, np.float64)): 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 prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row) # initialise psf kernel if ref: @@ -1754,6 +1761,13 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=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 From 235adbf1358fc04db557353f08e517b6c92e92a9 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Fri, 19 Sep 2025 15:33:08 +1200 Subject: [PATCH 28/31] last bad PRF call --- tessreduce/tessreduce.py | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 2fe2f76..53683bf 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -515,20 +515,35 @@ def psf_source_mask(self,sigma=5): """ + 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(np.round(loc[0],0)) + if isinstance(loc[1], (float, np.floating, np.float32, np.float64)): + 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 + # 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)) @@ -1783,7 +1798,7 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, 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) + localbkg_estimator=localbkg_estimator,xy_bounds=(0.05)) init = Table() init['x_init'] = [size//2] init['y_init'] = [size//2] From e32755f86a0d5bdd9c99080a6ae03e922310195e Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Tue, 23 Sep 2025 17:31:13 +1200 Subject: [PATCH 29/31] fixes with psf creation --- tessreduce/tessreduce.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 53683bf..ccefcd2 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -515,13 +515,8 @@ def psf_source_mask(self,sigma=5): """ - 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(np.round(loc[0],0)) - if isinstance(loc[1], (float, np.floating, np.float32, np.float64)): - loc[1] = int(np.round(loc[1])) + 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 @@ -2647,7 +2642,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. @@ -2694,7 +2689,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 From 302cf231b0324eeb2f5df30b8bcece85e57716de Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Sat, 8 Nov 2025 10:10:07 +1300 Subject: [PATCH 30/31] fixed bug with prf model indexing --- tessreduce/tessreduce.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index ccefcd2..d91b451 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -1674,8 +1674,13 @@ def _psf_initialise(self,cutoutSize,loc,time_ind=None,ref=False): col = 2090 if row > 2040: row = 2040 - - prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row) # initialise psf kernel + 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: @@ -2947,6 +2952,9 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None): 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 @@ -3008,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") From 69605400b2fbf39e7ccec05105cf7e5092a40a5b Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Tue, 11 Nov 2025 14:54:23 +1100 Subject: [PATCH 31/31] fixes to get ztf data and the sn_lookup --- tessreduce/ground_tools.py | 27 +++++++++++++------ tessreduce/helpers.py | 54 ++++++++++++++++++++------------------ 2 files changed, 47 insertions(+), 34 deletions(-) diff --git a/tessreduce/ground_tools.py b/tessreduce/ground_tools.py index 7ea9dea..53e05c8 100755 --- a/tessreduce/ground_tools.py +++ b/tessreduce/ground_tools.py @@ -78,7 +78,7 @@ def __init__(self,ra=None,dec=None,sn_name=None): self.ra = ra self.dec = dec self.sn_name = sn_name - + self.ztf_name = None # diags self.zp = -48.6 self.flux_type = 'mag' @@ -90,6 +90,17 @@ def __init__(self,ra=None,dec=None,sn_name=None): self.ps1 = None + def _get_ztf_name(self): + client = Alerce() + name = client.query_objects(ra=self.ra,dec=self.dec)['oid'].values + try: + name = name[0] + except: + pass + + self.ztf_name = name + + def __old_get_sn_name(self): """ If coordinates are known then get the transient name from OSC @@ -175,13 +186,13 @@ def get_ztf_data(self): """ Gets the ztf light curve data. First checks that the transient name is defined """ - if self.sn_name is None: - self.get_sn_name() - for name in self.cat_names: - if 'ZTF' in name: - ztf_name = name - if ztf_name is not None: - self.ztf = get_ztf(ztf_name) + if self.ztf_name is None: + self._get_ztf_name() + # for name in self.cat_names: + # if 'ZTF' in name: + # ztf_name = name + if self.ztf_name is not None: + self.ztf = get_ztf(self.ztf_name) return diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 774ca57..1918243 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -512,30 +512,32 @@ def sn_lookup(name,time='disc',buffer=0,print_table=True, df = False): tr_list : list list of ra, dec, and sector that can be put into tessreduce. """ - try: - url = 'https://api.astrocats.space/{}'.format(name) - response = requests.get(url) - json_acceptable_string = response.content.decode("utf-8").replace("'", "").split('\n')[0] - d = json.loads(json_acceptable_string) - if list(d.keys())[0] == 'message': - #print(d['message']) - - #return None - tns = True - else: - disc_t = d[name]['discoverdate'][0]['value'] - disc_t = Time(disc_t.replace('/','-')) - max_t = d[name]['maxdate'][0]['value'] - max_t = Time(max_t.replace('/','-')) - ra = d[name]['ra'][-1]['value'] - dec = d[name]['dec'][-1]['value'] - tns = False - except: - tns = True + tns = True + # try: + # url = 'https://api.astrocats.space/{}'.format(name) + # response = requests.get(url) + # json_acceptable_string = response.content.decode("utf-8").replace("'", "").split('\n')[0] + # d = json.loads(json_acceptable_string) + # if list(d.keys())[0] == 'message': + # #print(d['message']) + + # #return None + # tns = True + # else: + # disc_t = d[name]['discoverdate'][0]['value'] + # disc_t = Time(disc_t.replace('/','-')) + # max_t = d[name]['maxdate'][0]['value'] + # max_t = Time(max_t.replace('/','-')) + # ra = d[name]['ra'][-1]['value'] + # dec = d[name]['dec'][-1]['value'] + # tns = False + # except: + # tns = True if tns: #print('!! Open SNe Catalog down, using TNS !!') - name = name[name.index('2'):] - url = f'https://www.wis-tns.org/object/{name}' # hard coding in that the event is in the 2000s + if (name[:2].lower() == 'sn') | (name[:2].lower() == 'at'): + name = name[2:] + url = f'https://www.wis-tns.org/object/{name}' headers = {'User-Agent': 'Mozilla/5.0 (Macintosh; Intel Mac OS X 10_11_5) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/50.0.2661.102 Safari/537.36'} result = requests.get(url, headers=headers) if result.ok: @@ -544,11 +546,11 @@ def sn_lookup(name,time='disc',buffer=0,print_table=True, df = False): disc_t = Time(result.text.split('Discovery Date
')[-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