diff --git a/VERSION b/VERSION index f3b5581..13a05ac 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -master-v2.1 \ No newline at end of file +master-v2.1.1 \ No newline at end of file diff --git a/dataval/dataval.py b/dataval/dataval.py index 6d5418a..2ba965c 100644 --- a/dataval/dataval.py +++ b/dataval/dataval.py @@ -472,13 +472,15 @@ def validate(self): logger.info('--------------------------------------------------------') #---------------------------------------------------------------------------------------------- - def basic(self, warn_errors_ratio=0.05): + def basic(self, errors_ratio_warn=0.05, errors_ratio_err=0.10): """ Perform basic checks of the TODO-file and the lightcurve files. Parameters: - warn_errors_ratio (float, optional): Fraction of photometry ERRORs to OK and WARNINGs + errors_ratio_warn (float, optional): Fraction of ERRORs to OK and WARNINGs to warn about. Default=5%. + errors_ratio_err (float, optional): Fraction of ERRORs to OK and WARNINGs + to throw error about. Default=10%. .. codeauthor:: Rasmus Handberg """ @@ -528,8 +530,9 @@ def basic(self, warn_errors_ratio=0.05): )) count_errors = self.cursor.fetchone()[0] ratio = count_errors/(count_good + count_errors) if count_good + count_errors > 0 else 0 - if ratio > warn_errors_ratio: - logger.warning(" CAMERA=%d, CCD=%d: High number of errors detected: %.2f%% (%d errors, %d good)", + if ratio > errors_ratio_warn or ratio > errors_ratio_err: + loglevel = logging.ERROR if ratio > errors_ratio_err else logging.WARNING + logger.log(loglevel, " CAMERA=%d, CCD=%d: High number of errors detected: %.2f%% (%d errors, %d good)", camera, ccd, 100*ratio, count_errors, count_good) else: logger.info(" CAMERA=%d, CCD=%d: %.2f%% (%d errors, %d good)", @@ -585,8 +588,9 @@ def basic(self, warn_errors_ratio=0.05): )) count_errors = self.cursor.fetchone()[0] ratio = count_errors/(count_good + count_errors) if count_good + count_errors > 0 else 0 - if ratio > warn_errors_ratio: - logger.warning(" CAMERA=%d, CCD=%d: High number of errors detected: %.2f%% (%d errors, %d good)", + if ratio > errors_ratio_warn or ratio > errors_ratio_err: + loglevel = logging.ERROR if ratio > errors_ratio_err else logging.WARNING + logger.log(loglevel, " CAMERA=%d, CCD=%d: High number of errors detected: %.2f%% (%d errors, %d good)", camera, ccd, 100*ratio, count_errors, count_good) else: logger.info(" CAMERA=%d, CCD=%d: %.2f%% (%d errors, %d good)", @@ -606,7 +610,8 @@ def basic(self, warn_errors_ratio=0.05): 'FileNotFoundError', 'sqlite3.%', 'TargetNotFoundError', # custom "error" set in photometry.TaskManager.save_result - 'TypeError' + 'TypeError', + 'Could not save lightcurve file' ] logger.info("Checking for specific errors...") diff --git a/dataval/mag2flux.py b/dataval/mag2flux.py index ce28e11..f7dddaf 100644 --- a/dataval/mag2flux.py +++ b/dataval/mag2flux.py @@ -27,14 +27,20 @@ def mag2flux(dval): logger = logging.getLogger('dataval') logger.info('Plotting Magnitude to Flux conversion...') - zp_grid = np.linspace(19.0, 22.0, 100) - mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 200) + # Get the maximum median flux to be used to set the plot limits: + max_mean_flux = dval.search_database( + select=['MAX(mean_flux) AS maxflux'], + search=["method_used='aperture'"]) + max_mean_flux = max_mean_flux[0]['maxflux'] + + zp_grid = np.linspace(19.0, 22.0, 300) + mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 300) # The interpolation is linear *in log* xmin = np.array([0, 1.5, 9, 12.6, 13, 14, 15, 16, 17, 18, 19]) ymin = np.array([8e7, 1.8e7, 12500, 250, 59, 5, 1, 1, 1, 1, 1]) min_bound_log = InterpolatedUnivariateSpline(xmin, np.log10(ymin), k=1, ext=3) - min_bound = lambda x: 10**(min_bound_log(x)) # noqa: E731 + min_bound = lambda x: np.clip(10**(min_bound_log(x)), 1, None) # noqa: E731 norm = colors.Normalize(vmin=0, vmax=1) fig2, ax2 = plt.subplots() @@ -125,13 +131,14 @@ def chi2(c): #ax31.plot(bin_centers2, 1.4826*3*bin_means2, color='g') # Add line with best fit, and the minimim bound: - ax1.plot(mags, 10**(-0.4*(mags - cc.x)), color='k', ls='--') + ax1.plot(mags, np.clip(10**(-0.4*(mags - cc.x)), 0, None), color='k', ls='--') ax1.plot(mags, min_bound(mags), 'r-') ax1.set_yscale('log') ax1.set_xlim(dval.tmag_limits[1], dval.tmag_limits[0]) + ax1.set_ylim(0.5, 2*max_mean_flux) ax1.set_xlabel('TESS magnitude') - ax1.set_ylabel('Median flux') + ax1.set_ylabel('Median flux (e$^-$/s)') ax1.xaxis.set_major_locator(MultipleLocator(2)) ax1.xaxis.set_minor_locator(MultipleLocator(1)) colorbar(im, ax=ax1, label='Contamination') diff --git a/dataval/noise_metrics.py b/dataval/noise_metrics.py index aea61aa..f6c5fa6 100644 --- a/dataval/noise_metrics.py +++ b/dataval/noise_metrics.py @@ -28,7 +28,15 @@ def noise_metrics(dval): logger = logging.getLogger('dataval') logger.info('Plotting Noise vs. Magnitude...') - mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 200) + table_noise = 'diagnostics_corr' if dval.corr else 'diagnostics' + factor = 1 if dval.corr else 1e6 + max_vals = dval.search_database(select=[ + 'MAX(' + table_noise + '.rms_hour) AS max_rms', + 'MAX(' + table_noise + '.ptp) AS max_ptp']) + max_rms = factor*max_vals[0]['max_rms'] + max_ptp = factor*max_vals[0]['max_ptp'] + + mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 300) # Colors for theoretical lines: # Define the colors directly here to avoid having to import @@ -45,30 +53,16 @@ def noise_metrics(dval): for cadence in dval.cadences: - if dval.corr: - factor = 1 - star_vals = dval.search_database( - select=[ - 'todolist.priority', - 'todolist.sector', - 'todolist.tmag', - 'diagnostics_corr.rms_hour', - 'diagnostics_corr.ptp', - 'diagnostics.contamination' - ], - search=f'cadence={cadence:d}') - else: - factor = 1e6 - star_vals = dval.search_database( - select=[ - 'todolist.priority', - 'todolist.sector', - 'todolist.tmag', - 'diagnostics.rms_hour', - 'diagnostics.ptp', - 'diagnostics.contamination' - ], - search=f'cadence={cadence:d}') + star_vals = dval.search_database( + select=[ + 'todolist.priority', + 'todolist.sector', + 'todolist.tmag', + table_noise + '.rms_hour', + table_noise + '.ptp', + 'diagnostics.contamination' + ], + search=f'cadence={cadence:d}') tmags = np.array([star['tmag'] for star in star_vals], dtype='float64') pri = np.array([star['priority'] for star in star_vals], dtype='int64') @@ -130,6 +124,9 @@ def noise_metrics(dval): axx.set_yscale('log') #axx.legend(loc='upper left') + ax1.set_ylim(0.3*np.min(tot_noise_rms), 2*max_rms) + ax2.set_ylim(0.3*np.min(tot_noise_ptp), 2*max_ptp) + colorbar(im1, ax=ax1, label='Contamination') colorbar(im2, ax=ax2, label='Contamination') diff --git a/dataval/noise_model.py b/dataval/noise_model.py index 397a0b5..c0e24ce 100644 --- a/dataval/noise_model.py +++ b/dataval/noise_model.py @@ -10,6 +10,7 @@ from astropy import units as u from astropy.coordinates import SkyCoord import scipy.interpolate as INT +from .utilities import mag2flux #-------------------------------------------------------------------------------------------------- def ZLnoise(gal_lat): @@ -55,28 +56,6 @@ def Pixinaperture(Tmag, cad=1800): return np.asarray(np.maximum(pixels, 3), dtype='int32') -#-------------------------------------------------------------------------------------------------- -def mean_flux_level(Tmag): - """ - Mean flux from TESS magnitude - - .. codeauthor:: Mikkel N. Lund - """ - # Magnitude system based on Sullivan et al. - #collecting_area = np.pi*(10.5/2)**2 # square cm - #Teff_list = np.array([2450, 3000, 3200, 3400, 3700, 4100, 4500, 5000, 5777, 6500, 7200, 9700]) # Based on Sullivan - #Flux_list = np.array([2.38, 1.43, 1.40, 1.38, 1.39, 1.41, 1.43, 1.45, 1.45, 1.48, 1.48, 1.56])*1e6 # photons per sec; Based on Sullivan - #Magn_list = np.array([306, -191, -202, -201, -174, -132, -101, -80, -69.5, -40, -34.1, 35])*1e-3 #Ic-Tmag (mmag) - - #Flux_int = INT.UnivariateSpline(Teff_list, Flux_list, k=1, s=0) - #Magn_int = INT.UnivariateSpline(Teff_list, Magn_list, k=1, s=0) - - #Imag = Magn_int(Teff)+Tmag - #Flux = 10**(-0.4*Imag) * Flux_int(Teff) * collecting_area - - Flux = 10**(-0.4*(Tmag - 20.54)) - return Flux - #-------------------------------------------------------------------------------------------------- def phot_noise(Tmag, timescale=3600, coord=None, sysnoise=60, Teff=5775, cadpix=1800): """ @@ -136,7 +115,7 @@ def phot_noise(Tmag, timescale=3600, coord=None, sysnoise=60, Teff=5775, cadpix= Flux_factor = np.sqrt(integrations * pixels) # Mean flux level in electrons per cadence - mean_level_ppm = mean_flux_level(Tmag) * timescale # electrons (based on measurement) #, Teff + mean_level_ppm = mag2flux(Tmag) * timescale # electrons (based on measurement) #, Teff # Shot noise shot_noise = 1e6/np.sqrt(mean_level_ppm) @@ -152,8 +131,10 @@ def phot_noise(Tmag, timescale=3600, coord=None, sysnoise=60, Teff=5775, cadpix= # Put individual components together in single table: noise_vals = np.column_stack((shot_noise, zodiacal_noise, read_noise, systematic_noise)) + noise_vals = np.clip(noise_vals, 0, None) # Calculate the total noise model by adding up the individual contributions: total_noise = np.sqrt(np.sum(noise_vals**2, axis=1)) + total_noise = np.clip(total_noise, 0, None) return total_noise, noise_vals # ppm per cadence diff --git a/dataval/pixinaperture.py b/dataval/pixinaperture.py index cb1c2bb..d4ab58f 100644 --- a/dataval/pixinaperture.py +++ b/dataval/pixinaperture.py @@ -26,7 +26,13 @@ def pixinaperture(dval): logger = logging.getLogger('dataval') logger.info('Plotting Pixels in aperture vs. Magnitude...') - mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 200) + mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 500) + + # Find the maximim mask size, used as a constant upper limit on the plots: + max_masksize = dval.search_database( + select=['MAX(diagnostics.mask_size) AS maxsize'], + search=["method_used='aperture'", 'diagnostics.mask_size IS NOT NULL']) + max_masksize = max_masksize[0]['maxsize'] # TODO: Does this need to be separated by cadence as well? # Are there significant differences beween 120s and 20s mask sizes? @@ -130,26 +136,33 @@ def pixinaperture(dval): # The interpolations are linear *in log space* which is the reason # for the stange lambda wrappers in the following code: + # Rounding to 13 decimal places is to avoid numerical noise causing single pixel jumps if datasource == 'ffi': # Minimum bound on FFI data: xmin = np.array([0, 2, 2.7, 3.55, 4.2, 4.8, 5.5, 6.8, 7.6, 8.4, 9.1, 10, 10.5, 11, 11.5, 11.6, 16]) ymin = np.array([30, 30, 30, 30, 29, 28, 27, 26, 25, 22, 20, 15, 11.15, 8, 5, 4, 4]) min_bound_log = InterpolatedUnivariateSpline(xmin, np.log10(ymin), k=1, ext=3) - min_bound = lambda x: np.floor(10**(min_bound_log(x))) # noqa: E731 + min_bound = lambda x: np.clip(np.floor(np.round(10**(min_bound_log(x)), decimals=13)), 4, None) # noqa: E731 # Maximum bound on FFI data: xmax = np.array([0, 2, 2.7, 3.55, 4.2, 4.8, 5.5, 6.8, 7.6, 8.4, 9.1, 10, 10.5, 11, 11.5, 12, 12.7, 13.3, 14, 14.5, 15, 16]) ymax = np.array([10000, 3200, 2400, 1400, 1200, 900, 800, 470, 260, 200, 170, 130, 120, 100, 94, 86, 76, 67, 59, 54, 50, 50]) max_bound_log = InterpolatedUnivariateSpline(xmax, np.log10(ymax), k=1, ext=3) - max_bound = lambda x: np.ceil(10**(max_bound_log(x))) # noqa: E731 + max_bound = lambda x: np.ceil(np.round(10**(max_bound_log(x)), decimals=13)) # noqa: E731 else: # Minimum bound on TPF data xmin = np.array([0, 2, 2.7, 3.55, 4.2, 4.8, 5.5, 6.8, 7.6, 8.4, 9.1, 10, 10.5, 11, 11.5, 11.6, 16, 19]) ymin = np.array([220, 200, 130, 70, 55, 43, 36, 30, 27, 22, 16, 10, 8, 6, 5, 4, 4, 4]) min_bound_log = InterpolatedUnivariateSpline(xmin, np.log10(ymin), k=1, ext=3) - min_bound = lambda x: np.floor(10**(min_bound_log(x))) # noqa: E731 + min_bound = lambda x: np.clip(np.floor(np.round(10**(min_bound_log(x)), decimals=13)), 4, None) # noqa: E731 max_bound = None + # Simple sanity check if the bounds are monotonical decresing: + if not np.all(np.diff(min_bound(mags)) <= 0): + raise RuntimeError("Minimum bound (" + datasource + ") is not monotonically decreasing") # pragma: no cover + if max_bound is not None and not np.all(np.diff(max_bound(mags)) <= 0): + raise RuntimeError("Maximum bound (" + datasource + ") is not monotonically decreasing") # pragma: no cover + # Plot limits: ax1.plot(mags, min_bound(mags), 'r-') if max_bound is not None: @@ -208,8 +221,7 @@ def pixinaperture(dval): #ax2.plot(mags, pix, color='k', ls='-') ax1.set_xlim(dval.tmag_limits) - ax1.set_ylim([0.99, np.nanmax(masksizes)+500]) - + ax1.set_ylim([0.99, max_masksize+500]) ax1.set_xlabel('TESS magnitude') ax1.set_ylabel('Pixels in aperture') xtick_major = np.median(np.diff(ax1.get_xticks())) diff --git a/dataval/utilities.py b/dataval/utilities.py index e98f49c..3a0ac43 100644 --- a/dataval/utilities.py +++ b/dataval/utilities.py @@ -221,18 +221,20 @@ def mad(x): return mad_to_sigma * nanmedian(np.abs(x - nanmedian(x))) #-------------------------------------------------------------------------------------------------- -def mag2flux(mag, zp=20.60654144): +def mag2flux(mag, zp=20.451): """ Convert from magnitude to flux using scaling relation from aperture photometry. This is an estimate. - The scaling is based on fast-track TESS data from sectors 1 and 2. + The default scaling is based on TASOC Data Release 5 from sectors 1-5. Parameters: - mag (float): Magnitude in TESS band. + mag (ndarray): Magnitude in TESS band. + zp (float): Zero-point to use in scaling. Default is estimated from + TASOC Data Release 5 from TESS sectors 1-5. Returns: - float: Corresponding flux value + ndarray: Corresponding flux value """ return np.clip(10**(-0.4*(mag - zp)), 0, None)