Skip to content

Commit

Permalink
Merge pull request #47 from stscirij/fix_nogooddata
Browse files Browse the repository at this point in the history 
Add fix for HASP data sometimes creating a product with no good pixels
  • Loading branch information
@jotaylor
jotaylor authored Aug 25, 2023
2 parents bc3d3f5 + 85f2ac4 commit 97c2c79
Showing 1 changed file with 55 additions and 32 deletions.
87 changes: 55 additions & 32 deletions ullyses/coadd.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -191,8 +191,11 @@ def coadd(self):
weight = flux_weight[goodpixels]
net_counts = segment.data['net'][goodpixels] * segment.exptime
if self.instrument == 'COS':
variances = segment.data['variance_counts'][goodpixels] + segment.data['variance_bkg'][goodpixels] + segment.data['variance_flat'][goodpixels]

try:
variances = segment.data['variance_counts'][goodpixels] + segment.data['variance_bkg'][goodpixels] + segment.data['variance_flat'][goodpixels]
except KeyError:
print('WARNING: Cant get variance keywords for COS data, using gross counts instead')
variances = gcounts
flux = segment.data['flux'][goodpixels]
for i in range(npts):
self.sumnetcounts[indices[i]] = self.sumnetcounts[indices[i]] + net_counts[i]
Expand All @@ -203,30 +206,35 @@ def coadd(self):
if self.instrument == 'COS':
self.output_varsum[indices[i]] = self.output_varsum[indices[i]] + variances[i]
good_dq = np.where(self.output_exptime > 0.)
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
nonzeros = np.where(self.output_sumweight != 0)
if self.instrument == 'COS':
# Using the variances (which only COS has) gives spikes in the error when the flux goes negative.
self.output_varsum[nonzeros] = np.where(self.output_varsum[nonzeros] < 0.5, 0.5, self.output_varsum[nonzeros])
self.output_errors[nonzeros] = np.sqrt(self.output_varsum[nonzeros])
if len(good_dq[0]) == 0:
self.first_good_wavelength = None
self.last_good_wavelength = None
print("Warning: No good pixels in product")
else:
# For the moment calculate errors from the gross counts
self.output_errors[nonzeros] = np.sqrt(self.output_sumgcounts[nonzeros])
self.output_flux[nonzeros] = self.output_sumflux[nonzeros] / self.output_sumweight[nonzeros]
# Calculate the conversion from counts to flux units for errors
conversion = self.output_flux[nonzeros] / self.sumnetcounts[nonzeros]
# Clean out NaNs from where flux and net are zero
good = np.where(~np.isnan(conversion))
bad = np.where(np.isnan(conversion))
lenwl = len(wavelength)
interpolated_values = np.interp(self.output_wavelength[bad],
self.output_wavelength[good],
conversion[good])
conversion[bad] = interpolated_values
# Use conversion to calculate error in flux units (erg/cm^2/s/Ang)
self.output_errors[nonzeros] = self.output_errors[nonzeros] * conversion
self.signal_to_noise[nonzeros] = self.output_flux[nonzeros] / self.output_errors[nonzeros]
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
nonzeros = np.where(self.output_sumweight != 0)
if self.instrument == 'COS':
# Using the variances (which only COS has) gives spikes in the error when the flux goes negative.
self.output_varsum[nonzeros] = np.where(self.output_varsum[nonzeros] < 0.5, 0.5, self.output_varsum[nonzeros])
self.output_errors[nonzeros] = np.sqrt(self.output_varsum[nonzeros])
else:
# For the moment calculate errors from the gross counts
self.output_errors[nonzeros] = np.sqrt(self.output_sumgcounts[nonzeros])
self.output_flux[nonzeros] = self.output_sumflux[nonzeros] / self.output_sumweight[nonzeros]
# Calculate the conversion from counts to flux units for errors
conversion = self.output_flux[nonzeros] / self.sumnetcounts[nonzeros]
# Clean out NaNs from where flux and net are zero
good = np.where(~np.isnan(conversion))
bad = np.where(np.isnan(conversion))
lenwl = len(wavelength)
interpolated_values = np.interp(self.output_wavelength[bad],
self.output_wavelength[good],
conversion[good])
conversion[bad] = interpolated_values
# Use conversion to calculate error in flux units (erg/cm^2/s/Ang)
self.output_errors[nonzeros] = self.output_errors[nonzeros] * conversion
self.signal_to_noise[nonzeros] = self.output_flux[nonzeros] / self.output_errors[nonzeros]
return

def get_targname(self):
Expand Down Expand Up @@ -339,8 +347,13 @@ def coadd(self):
nonzeros = np.where(self.output_errors != 0.0)
self.signal_to_noise[nonzeros] = self.output_flux[nonzeros] / self.output_errors[nonzeros]
good_dq = np.where(self.output_exptime > 0.)
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
if len(good_dq[0]) > 0:
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
else:
self.first_good_wavelength = None
self.last_good_wavelength = None
print('No good data in product')
return

class CCDSegmentList(SegmentList):
Expand Down Expand Up @@ -395,7 +408,7 @@ def create_output_wavelength_grid(self):

self.delta_wavelength = max_delta_wavelength

wavegrid = np.arange(self.min_wavelength, self.max_wavelength, self.delta_wavelength)
wavegrid = np.arange(self.min_wavelength, self.max_wavelength+self.delta_wavelength, self.delta_wavelength)

self.output_wavelength = wavegrid

Expand Down Expand Up @@ -425,8 +438,13 @@ def coadd(self):
nonzeros = np.where(self.output_errors != 0.0)
self.signal_to_noise[nonzeros] = self.output_flux[nonzeros] / self.output_errors[nonzeros]
good_dq = np.where(self.output_exptime > 0.)
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
if len(good_dq[0]) > 0:
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
else:
self.first_good_wavelength = None
self.last_good_wavelength = None
print('No good data in product')
else:
nelements = len(self.output_wavelength)
self.output_sumsqerrors = np.zeros(nelements)
Expand All @@ -444,8 +462,13 @@ def coadd(self):
self.output_sumsqerrors[indices[i]] = self.output_sumsqerrors[indices[i]] + (err[i] * weight[i])**2
self.output_exptime[indices[i]] = self.output_exptime[indices[i]] + segment.exptime
good_dq = np.where(self.output_exptime > 0.)
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
if len(good_dq[0]) > 0:
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
else:
self.first_good_wavelength = None
self.last_good_wavelength = None
print('No good data in product')
nonzeros = np.where(self.output_sumweight != 0)
self.output_flux[nonzeros] = self.output_sumflux[nonzeros] / self.output_sumweight[nonzeros]
self.output_errors[nonzeros] = np.sqrt(self.output_sumsqerrors[nonzeros]) / self.output_sumweight[nonzeros]
Expand Down

0 comments on commit 97c2c79

Please sign in to comment.