Skip to content

Commit c283a32

Browse files
CheerfulUserCheerfulUser
committed
winding back qe changes...
1 parent 44ab58d commit c283a32

File tree

2 files changed

+167
-118
lines changed

2 files changed

+167
-118
lines changed

‎tessreduce/helpers.py‎

Lines changed: 73 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -53,18 +53,18 @@
5353

5454
def strip_units(data):
5555
"""
56-
Removes the units off of data that was not in a NDarray, such as an astropy table. Returns an NDarray that has no units
57-
58-
Parameters:
59-
----------
60-
data: ArrayLike
61-
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.
62-
63-
Returns:
64-
-------
65-
data: ArrayLike
66-
Same shape as input data, but will not have any units
67-
"""
56+
Removes the units off of data that was not in a NDarray, such as an astropy table. Returns an NDarray that has no units
57+
58+
Parameters:
59+
----------
60+
data: ArrayLike
61+
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.
62+
63+
Returns:
64+
-------
65+
data: ArrayLike
66+
Same shape as input data, but will not have any units
67+
"""
6868

6969
if type(data) != np.ndarray:
7070
data = data.value
@@ -294,25 +294,25 @@ def image_sub(theta, image, ref, eimage, eref):
294294

295295
### TESTING CODE
296296
#def c_shift_image(img, shift):
297-
# return scipy.ndimage.shift(img, shift=shift, order=5, mode='nearest') # cubic interpolation
297+
# return scipy.ndimage.shift(img, shift=shift, order=5, mode='nearest') # cubic interpolation
298298
#
299299
#def cost_function(shift, ref_img, moving_img):
300-
# shifted = shift_image(moving_img, shift)
301-
# diff = ref_img - shifted
302-
# return np.sum(diff[2:-2,2:-2]**2) # Sum of Squared Differences (SSD)
300+
# shifted = shift_image(moving_img, shift)
301+
# diff = ref_img - shifted
302+
# return np.sum(diff[2:-2,2:-2]**2) # Sum of Squared Differences (SSD)
303303
#
304304
#def align_subpixel(ref_img, moving_img, initial_shift=(0,0)):
305-
# ref_img = ref_img.astype(np.float32)
306-
# moving_img = moving_img.astype(np.float32)
305+
# ref_img = ref_img.astype(np.float32)
306+
# moving_img = moving_img.astype(np.float32)
307307
#
308-
# # Minimize SSD by shifting moving_img
309-
# result = minimize(cost_function, initial_shift, args=(ref_img, moving_img), method='Powell')
308+
# # Minimize SSD by shifting moving_img
309+
# result = minimize(cost_function, initial_shift, args=(ref_img, moving_img), method='Powell')
310310
#
311-
# best_shift = result.x
312-
# aligned_img = c_shift_image(moving_img, best_shift)
311+
# best_shift = result.x
312+
# aligned_img = c_shift_image(moving_img, best_shift)
313313
#
314-
# print(f"Optimal shift (y, x): {best_shift}")
315-
# return aligned_img, best_shift
314+
# print(f"Optimal shift (y, x): {best_shift}")
315+
# return aligned_img, best_shift
316316
###
317317

318318

@@ -450,36 +450,36 @@ def smooth_zp(zp,time):
450450
return smoothed, err
451451

452452
def grads_rad(flux):
453-
"""
454-
Calculates the radius of the flux from the gradient of the flux, and the double gradient of the flux.
455-
456-
Parameters:
457-
----------
458-
flux: ArrayLike
459-
An array of flux values
460-
461-
Returns:
462-
-------
463-
rad: ArrayLike
464-
The radius of the fluxes
465-
"""
466-
rad = np.sqrt(np.gradient(flux)**2+np.gradient(np.gradient(flux))**2)
467-
return rad
453+
"""
454+
Calculates the radius of the flux from the gradient of the flux, and the double gradient of the flux.
455+
456+
Parameters:
457+
----------
458+
flux: ArrayLike
459+
An array of flux values
460+
461+
Returns:
462+
-------
463+
rad: ArrayLike
464+
The radius of the fluxes
465+
"""
466+
rad = np.sqrt(np.gradient(flux)**2+np.gradient(np.gradient(flux))**2)
467+
return rad
468468

469469
def grad_flux_rad(flux):
470470
"""
471-
Calculates the radius of the flux from the gradient of the flux.
472-
473-
Parameters:
474-
----------
475-
flux: ArrayLike
476-
An array of flux values
477-
478-
Returns:
479-
-------
480-
rad: ArrayLike
481-
The radius of the fluxes
482-
"""
471+
Calculates the radius of the flux from the gradient of the flux.
472+
473+
Parameters:
474+
----------
475+
flux: ArrayLike
476+
An array of flux values
477+
478+
Returns:
479+
-------
480+
rad: ArrayLike
481+
The radius of the fluxes
482+
"""
483483
rad = np.sqrt(flux**2+np.gradient(flux)**2)
484484
return rad
485485

@@ -1249,15 +1249,15 @@ def grad_clip(data,box_size=100):
12491249
12501250
"""
12511251
gradind = np.zeros_like(data)
1252-
1252+
12531253
for i in range(len(data)):
12541254
if i < box_size//2:
12551255
d = data[:i+box_size//2]
12561256
elif len(data) - i < box_size//2:
12571257
d = data[i-box_size//2:]
12581258
else:
12591259
d = data[i-box_size//2:i+box_size//2]
1260-
1260+
12611261
ind = np.isfinite(d)
12621262
d = d[ind]
12631263
if len(d) > 5:
@@ -1268,7 +1268,7 @@ def grad_clip(data,box_size=100):
12681268
gradind[i-box_size//2:][ind] = gind
12691269
else:
12701270
gradind[i-box_size//2:i+box_size//2][ind] = gind
1271-
1271+
12721272
gradind = gradind > 0
12731273
return gradind
12741274

@@ -1337,15 +1337,25 @@ def simulate_epsf(camera,ccd,sector,column,row,size=13,realizations=2000,oversam
13371337

13381338
return epsf
13391339

1340-
def parallel_photutils(cutout,e_cutout,psf_phot,init_params=None):
1341-
if np.nansum(abs(cutout)) > 0:
1342-
phot = psf_phot(cutout, error=e_cutout,init_params=init_params)
1343-
phot = phot.to_pandas()
1344-
return phot[['flux_fit']].values[0], phot[['flux_err']].values[0]
1345-
#return float(phot[['flux_fit']].values), float(phot[['flux_err']].values)
1346-
else:
1347-
return np.array([np.nan]), np.array([np.nan])
1348-
1340+
def parallel_photutils(cutout,e_cutout,psf_phot,init_params=None,return_pos=False):
1341+
if np.nansum(abs(cutout)) > 0:
1342+
phot = psf_phot(cutout, error=e_cutout,init_params=init_params)
1343+
phot = phot.to_pandas()
1344+
f = phot[['flux_fit']].values[0]; ef = phot[['flux_err']].values[0]
1345+
pos = phot[['x_fit','y_fit']].values
1346+
epos = phot[['x_err','y_err']].values
1347+
if return_pos:
1348+
return f, ef, pos, epos
1349+
else:
1350+
return f, ef
1351+
#return float(phot[['flux_fit']].values), float(phot[['flux_err']].values)
1352+
else:
1353+
if return_pos:
1354+
pos = np.array([np.nan,np.nan])
1355+
return np.array([np.nan]), np.array([np.nan]), pos, pos
1356+
else:
1357+
return np.array([np.nan]), np.array([np.nan])
1358+
13491359

13501360

13511361

‎tessreduce/tessreduce.py‎

Lines changed: 94 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -547,38 +547,53 @@ def _calc_qe(self,flux,bkg_smth,flux_e):
547547
Calculate the effective quantum efficiency enhancement of the detector from scattered light.
548548
'''
549549
norm = flux / bkg_smth
550-
straps = norm * ((self.mask & 4) > 0)
550+
straps = norm * ((self.mask & 4)>0)
551551

552552
straps[straps==0] = np.nan
553553
m,med,std = sigma_clipped_stats(straps,axis=(1,2),maxiters=10,mask_value=np.nan)
554554
masks = straps < (med + 2*std)[:,np.newaxis,np.newaxis]
555555
m = (np.nansum(masks,axis=0) > 0) * 1.
556556
m[m==0] = np.nan
557-
qe = np.ones_like(flux) * 1.
558-
if flux.shape[1] < 10000:
559-
ratio = (flux / bkg_smth) * m
560-
561-
m, med, std = sigma_clipped_stats(ratio,axis=1,sigma_upper=2,sigma_lower=3)
562-
#qe_1d = np.nanpercentile(ratio,10,axis=1)
563-
564-
#qe[:,:,:] = qe_1d[:,np.newaxis,:]
565-
qe[:,:,:] = med[:,np.newaxis,:]
566-
567-
else:
568-
index = np.arange(len(flux),dtype=int)
569-
if self.parallel:
570-
qe = Parallel(n_jobs=self.num_cores)(delayed(parallel_strap_fit)(flux[i],bkg_smth[i],flux_e[i],m) for i in index)
571-
else:
572-
qe = []
573-
for i in index:
574-
qe += [parallel_strap_fit(flux[i],bkg_smth[i],flux_e[i],m)]
575-
qe = np.array(qe)
576-
#filler = np.nanmedian(qe,axis=1)
577-
557+
ratio = flux / bkg_smth * m
558+
#m, med, std = sigma_clipped_stats(ratio,axis=1,sigma_upper=2)
559+
qe_1d = np.nanpercentile(ratio,10,axis=1)
560+
qe = np.ones_like(norm)
561+
qe[:,:,:] = qe_1d[:,np.newaxis,:]
562+
#qe[:,:,:] = med[:,np.newaxis,:]
578563
qe[np.isnan(qe)] = 1
579564
qe[qe<1] = 1
580565
return qe
581566

567+
# straps[straps==0] = np.nan
568+
# m,med,std = sigma_clipped_stats(straps,axis=(1,2),maxiters=10,mask_value=np.nan)
569+
# masks = straps < (med + 2*std)[:,np.newaxis,np.newaxis]
570+
# m = (np.nansum(masks,axis=0) > 0) * 1.
571+
# m[m==0] = np.nan
572+
# qe = np.ones_like(flux) * 1.
573+
# if flux.shape[1] < 10000:
574+
# ratio = (flux / bkg_smth) * m
575+
576+
# m, med, std = sigma_clipped_stats(ratio,axis=1,sigma_upper=2,sigma_lower=3)
577+
# #qe_1d = np.nanpercentile(ratio,10,axis=1)
578+
579+
# #qe[:,:,:] = qe_1d[:,np.newaxis,:]
580+
# qe[:,:,:] = med[:,np.newaxis,:]
581+
582+
# else:
583+
# index = np.arange(len(flux),dtype=int)
584+
# if self.parallel:
585+
# qe = Parallel(n_jobs=self.num_cores)(delayed(parallel_strap_fit)(flux[i],bkg_smth[i],flux_e[i],m) for i in index)
586+
# else:
587+
# qe = []
588+
# for i in index:
589+
# qe += [parallel_strap_fit(flux[i],bkg_smth[i],flux_e[i],m)]
590+
# qe = np.array(qe)
591+
# #filler = np.nanmedian(qe,axis=1)
592+
593+
# qe[np.isnan(qe)] = 1
594+
# qe[qe<1] = 1
595+
# return qe
596+
582597
def background(self,gauss_smooth=2,calc_qe=True,strap_iso=True,source_hunt=False,interpolate=True,rerun_negative=False):
583598
"""
584599
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):
16911706
pos[0,:] += xpos; pos[1,:] += ypos
16921707
return flux, pos
16931708

1694-
def psf_photutils(self,xPix,yPix,size=5,local_bkg=False,ref=False,return_pos=False):
1709+
def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None,
1710+
flux=None,eflux=None,ref=False,return_pos=False,
1711+
snap='brightest'):
16951712
from photutils.psf import PSFPhotometry
16961713
from astropy.table import Table
16971714
rad = size // 2
1698-
if self.epsf is None:
1699-
col = self.tpf.column - int(self.size//2) + yPix # find column and row, when specifying location on a *say* 90x90 px cutout
1700-
row = self.tpf.row - int(self.size//2) + xPix
1701-
self.epsf = simulate_epsf(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row)
1715+
if epsf is None:
1716+
if self.epsf is None:
1717+
col = self.tpf.column - int(self.size//2) + yPix # find column and row, when specifying location on a *say* 90x90 px cutout
1718+
row = self.tpf.row - int(self.size//2) + xPix
1719+
self.epsf = simulate_epsf(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row)
1720+
epsf = self.epsf
1721+
17021722
if local_bkg:
17031723
localbkg_estimator = LocalBackground(1.5, 7, bkgstat)
17041724
else:
17051725
localbkg_estimator = None
1706-
if ref:
1707-
flux = self.flux + self.ref
1708-
else:
1726+
if flux is None:
17091727
flux = self.flux
1710-
cutouts = flux[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1]
1711-
ecutouts = self.eflux[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1]
1728+
if (flux.shape[1] < size) | (flux.shape[2] < size):
1729+
e = 'Image dimensions must be larger than the cutout size'
1730+
raise ValueError(e)
1731+
if (xPix is None) | (yPix is None):
1732+
xPix = flux.shape[2]//2
1733+
yPix = flux.shape[1]//2
1734+
if eflux is None:
1735+
eflux = self.eflux
1736+
if ref:
1737+
flux = flux + self.ref
17121738

1713-
weight = np.abs(np.nansum((cutouts/ecutouts)[:,1:4,1:4],axis=(1,2)))
1714-
weight[np.isnan(weight)] = 0
1715-
ind = np.argmax(weight)
1739+
cutouts = flux[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1]
1740+
ecutouts = eflux[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1]
17161741
fit_shape = (size, size)
1717-
# initial position fit on brightest frame
1718-
psfphot = PSFPhotometry(self.epsf, fit_shape, finder=None,aperture_radius=1.5,
1742+
psfphot = PSFPhotometry(epsf, fit_shape, finder=None,aperture_radius=1.5,
17191743
localbkg_estimator=localbkg_estimator)
17201744
init = Table()
17211745
init['x_init'] = [size//2]
17221746
init['y_init'] = [size//2]
1723-
phot = psfphot(cutouts[ind], error=ecutouts[ind],init_params=init)
1724-
1725-
# fit with the best position
1726-
init2 = Table()
1727-
init['x_init'] = phot['x_fit']
1728-
init['y_init'] = phot['y_fit']
1729-
flux = np.zeros(len(self.flux)) * np.nan
1730-
eflux = np.zeros(len(self.flux)) * np.nan
1731-
psfphot2 = PSFPhotometry(self.epsf, fit_shape, finder=None,aperture_radius=1.5,
1732-
xy_bounds=(0.05),localbkg_estimator=localbkg_estimator)
1733-
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))))
1734-
f = np.array(f).flatten()
1735-
ef = np.array(ef).flatten()
1736-
1737-
phot = phot.to_pandas()
1738-
pos = phot[['x_fit','y_fit']].values + np.array([xPix,yPix]) - size//2
1739-
epos = phot[['x_err','y_err']].values
1747+
if snap.lower() != 'all':
1748+
if snap.lower() == 'brightest':
1749+
weight = np.abs(np.nansum((cutouts/ecutouts)[:,1:4,1:4],axis=(1,2)))
1750+
weight[np.isnan(weight)] = 0
1751+
ind = np.argmax(weight)
1752+
rcut = cutouts[ind]
1753+
ecut = ecutouts[ind]
1754+
elif snap.lower() == 'ref':
1755+
rcut = self.ref[:,yPix-rad:yPix+rad+1,xPix-rad:xPix+rad+1]
1756+
ecut = ecutouts[self.ref_ind]
1757+
1758+
1759+
phot = psfphot(rcut, error=ecut,init_params=init)
1760+
# fit with the best position
1761+
init2 = Table()
1762+
init['x_init'] = phot['x_fit']
1763+
init['y_init'] = phot['y_fit']
1764+
flux = np.zeros(len(self.flux)) * np.nan
1765+
eflux = np.zeros(len(self.flux)) * np.nan
1766+
psfphot2 = PSFPhotometry(epsf, fit_shape, finder=None,aperture_radius=1.5,
1767+
xy_bounds=(0.05),localbkg_estimator=localbkg_estimator)
1768+
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))))
1769+
f = np.array(f).flatten()
1770+
ef = np.array(ef).flatten()
1771+
phot = phot.to_pandas()
1772+
pos = phot[['x_fit','y_fit']].values + np.array([xPix,yPix]) - size//2
1773+
epos = phot[['x_err','y_err']].values
1774+
else:
1775+
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))))
1776+
pos['x_fit'] += xPix - size//2
1777+
pos['y_fit'] += xPix - size//2
1778+
17401779
if return_pos:
17411780
return f, ef, pos, epos
17421781
else:
@@ -1821,7 +1860,7 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa
18211860
bkg = self.bkg[:,int(yPix),int(xPix)]
18221861
#lowbkg = bkg < np.nanpercentile(bkg,16)
18231862
#weight = np.abs(np.nansum(cutouts[:,int(yPix)-1:int(yPix)+2,int(xPix)-1:int(xPix)+2],axis=(1,2))) / bkg
1824-
weight = np.abs(np.nansum((cutouts / ecuts)[:,int(yPix)-1:int(yPix)+2,int(xPix)-1:int(xPix)+2],axis=(1,2)))
1863+
weight = np.abs(np.nansum((cutouts / ecutouts)[:,int(yPix)-1:int(yPix)+2,int(xPix)-1:int(xPix)+2],axis=(1,2)))
18251864
weight[np.isnan(weight)] = 0
18261865
ind = np.argmax(weight)
18271866
#ind = np.where(cutouts==np.nanmax(cutouts))[0][0]

0 commit comments

Comments
 (0)