Merge pull request #53 from duxfrederic/minor_dev

- added option to skip background subtraction
- common star footprint star query fixed
- switching from panstarrs dr2 to dr1 due to problems in querying dr2. 
- oversized keywords handling for lsst headers
- gaia star position plot: plotting roi as well

Author: duxfrederic

lightcurver/pipeline/example_config_file/config.yaml

@@ -70,11 +70,12 @@ telescope:
 ############################################# Background estimation ####################################################
 
 # We need to very precisely subtract the background. The quality of the photometry depends immensely on this.
-# `sep` does a great job at this, but we offer the possibility to first mask the sources before estimating the
-# background (in case of crowded fields)
+do_background_subtraction: true  # set to false if images are already background subtracted.
+# `sep` does a great job at background subraction, but we offer the possibility to first mask the sources before
+# background estimation (in case of crowded fields)
 mask_sources_before_background: false
 # also, if you have strong gradients, you might want to allow your image to be subdivided in more boxes.
-background_estimation_n_boxes: 10
+background_estimation_n_boxes: 3
 # (this is how many boxes we divide the side of the image in, so in total there'll be the square of this number boxes)
 
 

lightcurver/plotting/sources_plotting.py

@@ -96,6 +96,8 @@ def plot_footprints_with_stars(footprint_arrays, stars, save_path=None):
     """
     fig, ax = plot_footprints(footprint_arrays, common_footprint=None, largest_footprint=None, save_path=None)
     for _, star in stars.iterrows():
+        if star['name'] == 'roi':
+            ax.plot(star['ra'], star['dec'], 'o', color='red', markersize=10, mfc='None')
         ax.plot(star['ra'], star['dec'], 'o', color='red', markersize=5, mfc='None')
         ax.text(star['ra'], star['dec'], star['name'], fontsize=8, ha='right')
 

lightcurver/processes/frame_importation.py

@@ -77,29 +77,47 @@ def process_new_frame(fits_file, user_config):
     cutout_data *= mjd_gain_exptime_dict['gain'] / mjd_gain_exptime_dict['exptime']
     # unit: electron per second
 
-    # ok, now subtract sky!
+    # ok, now estimate and subtract sky!
     cutout_data_sub, bkg = subtract_background(cutout_data,
                                                mask_sources_first=user_config['mask_sources_before_background'],
                                                n_boxes=user_config['background_estimation_n_boxes'])
+    # we keep doing this because it is a nice way of getting the background RMS.
+    # however, if the user specifies that the background is already subtracted, just use the 
+    # original data.
+    if not user_config['do_background_subtraction']:
+        cutout_data_sub = cutout_data
+    # read the background stats:
     sky_level_electron_per_second = float(bkg.globalback)
     background_rms_electron_per_second = float(bkg.globalrms)
     logger.info(f'  file {fits_file}: background estimated.')
 
-    # before we write, let's keep as much as we can from our previous header
-    new_header = wcs.to_header()
-    # then copy non-WCS entries from the original header
-    for key in header:
-        if key not in new_header and not key.startswith('WCSAXES') and not key.startswith(
-                'CRPIX') and not key.startswith('CRVAL') and not key.startswith('CDELT') and not key.startswith(
-                'CTYPE') and not key.startswith('CUNIT') and not key.startswith('CD') and key not in ['COMMENT',
-                                                                                                      'HISTORY']:
-            if key.strip():  # some headers are weird
-                new_header[key] = header[key]
+    new_header = header.copy()
+    # here we branch: if already plate solved, we keep the header intact.
+    # if not, we delete wcs keywords and populate an initial one from astropy
+    if not user_config['already_plate_solved']:
+        wcs_header = wcs.to_header()
+        # then copy non-WCS entries from the original header
+        for key in header:
+            if (
+                     key in wcs_header
+                  or key.startswith('WCSAXES')
+                  or key.startswith('CRPIX')
+                  or key.startswith('CRVAL')
+                  or key.startswith('CDELT')
+                  or key.startswith('CTYPE')
+                  or key.startswith('CUNIT')
+                  or key.startswith('CD')
+                  or key.startswith('A_')
+                  or key.startswith('B_')
+            ):
+                del new_header[key]
+        #
+        new_header.update(wcs_header)
     # now we can write the file
     write_file = user_config['workdir'] / copied_image_relpath
     logger.info(f'    Writing file: {write_file}')
     fits.writeto(write_file, cutout_data_sub.astype(np.float32),
-                 header=new_header, overwrite=True)
+                 header=new_header, overwrite=True, output_verify="silentfix")
     # and find sources
     # (do plots if toggle set)
     do_plot = user_config.get('source_extraction_do_plots', False)

lightcurver/processes/psf_modelling.py

@@ -173,12 +173,17 @@ def model_all_psfs():
         psf_plots_dir.mkdir(exist_ok=True, parents=True)
         frame_name = Path(frame['image_relpath']).stem
         seeing = frame['seeing_pixels'] * frame['pixel_scale']
+        # what's the fitted fwhm here?
+        kwargs_moffat = result['kwargs_psf']['kwargs_moffat']
+        fwhm_moffat_arcseconds = 0.5 * (kwargs_moffat['fwhm_x'] + kwargs_moffat['fwhm_y']) * frame['pixel_scale']
+        fwhm_moffat_arcseconds = float(fwhm_moffat_arcseconds.item()) # from jax array to float
         loss_history = result['adabelief_extra_fields']['loss_history']
+        diagnostic_text = f"{frame_name}\nseeing estimation: {seeing:.02f}\nseeing moffat: {fwhm_moffat_arcseconds:.02f}"
         plot_psf_diagnostic(datas=datas, noisemaps=noisemaps, residuals=result['residuals'],
                             full_psf=result['full_psf'],
                             loss_curve=loss_history,
                             masks=masks, names=names,
-                            diagnostic_text=f"{frame_name}\nseeing: {seeing:.02f}",
+                            diagnostic_text=diagnostic_text,
                             save_path=psf_plots_dir / f"{frame['id']}_{frame_name}.jpg")
 
         # now we can do the bookkeeping stuff
@@ -202,10 +207,11 @@ def model_all_psfs():
         end_change = np.nanmax(loss_history[loss_index:]) - np.nanmin(loss_history[loss_index:])
         relative_loss_differential = float(end_change / initial_change)
         insert_params = (frame['id'], float(result['chi2']), relative_loss_differential, psf_ref,
-                         combined_footprint_hash, user_config['subsampling_factor'])
+                         combined_footprint_hash, user_config['subsampling_factor'], fwhm_moffat_arcseconds)
         insert_query = f"""
         REPLACE INTO PSFs 
-        (frame_id, chi2, relative_loss_differential, psf_ref, combined_footprint_hash, subsampling_factor)
+           (frame_id, chi2, relative_loss_differential, psf_ref, 
+            combined_footprint_hash, subsampling_factor, fwhm_moffat_arcseconds)
         VALUES ({','.join(['?'] * len(insert_params))})
         """
         execute_sqlite_query(insert_query, insert_params, is_select=False)

lightcurver/processes/roi_file_preparation.py

@@ -108,8 +108,12 @@ def fetch_and_adjust_zeropoints(combined_footprint_hash):
     zp_scatter_normalized = zeropoints_data['adjusted_zeropoint'].std()
     global_zp = zeropoints_data['adjusted_zeropoint'].median()
 
-    message = "The scatter in zeropoints before normalizing is lower than after normalizing? Not normal, investigate."
-    #assert zp_scatter_normalized < zp_scatter_not_normalized, message
+    if zp_scatter_normalized > zp_scatter_not_normalized:
+        # This shouldn't happen in well controlled case. (normalising should result in a very stable zeropoint)
+        # If might happen harmlessly for very low frame numbers.
+        message = "The scatter in zeropoints before normalizing is lower than after normalizing? Not normal, investigate."
+        logger = logging.getLogger('lightcurver.roi_file_preparation')
+        logger.warning(message)
 
     return global_zp, zp_scatter_normalized
 

lightcurver/processes/star_photometry.py

@@ -231,9 +231,9 @@ def update_star_fluxes(flux_data):
 
 def do_star_photometry():
     """
-    Here we do a joint deconvolution of all the frames at once, for each star.
+    Here we do a joint modelling of all the frames at once, for each star.
     This is equivalent to PSF photometry, but we do it in the same way
-    as the final deconvolution to eliminate potential systematics.
+    as the final modelling to eliminate potential systematics.
     Returns:
 
     """
@@ -325,11 +325,11 @@ def do_star_photometry():
             starlet_global_background=user_config['star_photometry_starlet_global_background']
         )
         # ok, plot the diagnostic
-        plot_deconv_dir = user_config['plots_dir'] / 'deconvolutions' / str(combined_footprint_hash)
-        plot_deconv_dir.mkdir(exist_ok=True, parents=True)
+        plot_star_modelling_dir = user_config['plots_dir'] / 'star_modelling' / str(combined_footprint_hash)
+        plot_star_modelling_dir.mkdir(exist_ok=True, parents=True)
         loss_history = result['loss_curve']
 
-        plot_file = plot_deconv_dir / f"{time_now}_joint_deconv_star_{star['name']}.jpg"
+        plot_file = plot_star_modelling_dir / f"{time_now}_joint_modelling_star_{star['name']}.jpg"
 
         kwargs_plot = {
             'datas': data,

lightcurver/processes/star_querying.py

@@ -2,6 +2,7 @@
 from astropy.coordinates import SkyCoord
 import json
 import numpy as np
+import pandas as pd
 
 from ..utilities.footprint import load_combined_footprint_from_db, get_frames_hash
 from ..structure.user_config import get_user_config
@@ -120,5 +121,9 @@ def query_gaia_stars():
     results = execute_sqlite_query(query)
     polygon_list = [np.array(json.loads(result[1])) for result in results]
     save_path = user_config['plots_dir'] / 'footprints_with_gaia_stars.jpg'
-    plot_footprints_with_stars(footprint_arrays=polygon_list, stars=stars_table.to_pandas(), save_path=save_path)
+    stars_table = stars_table.to_pandas()
+    roi_coord = user_config['ROI_SkyCoord']
+    roi_row = {'name': 'roi', 'ra': roi_coord.ra.value, 'dec': roi_coord.dec.value}
+    stars_table = pd.concat([stars_table, pd.DataFrame([roi_row])], ignore_index=True)
+    plot_footprints_with_stars(footprint_arrays=polygon_list, stars=stars_table, save_path=save_path)
     logger.info(f'Plot of the queried reference Gaia stars saved at {save_path}.')

lightcurver/structure/database.py

@@ -385,6 +385,7 @@ def initialize_database(db_path=None):
                       psf_ref TEXT, -- convention: sorted concatenation of all star names used in the model.
                       subsampling_factor INTEGER,  -- we do starred (pixelated) PSFs.
                       relative_loss_differential REAL, -- absolute change in last 10% of the iterations vs beginning
+                      fwhm_moffat_arcseconds REAL DEFAULT NULL, -- keeping track of image quality
                       FOREIGN KEY (frame_id) REFERENCES frames(id),
                       FOREIGN KEY (combined_footprint_hash) REFERENCES combined_footprint(hash),
                       PRIMARY KEY (combined_footprint_hash, frame_id, psf_ref)

lightcurver/utilities/absolute_magnitudes_from_panstarrs.py

@@ -73,7 +73,7 @@ def search_panstarrs_around_coordinates(gaia_id):
     coord = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, frame='icrs')
     radius = 1.5 * u.arcsecond  # this is generous given the magnitude of the proper motion of most stars.
     try:
-        result = Catalogs.query_region(coord, radius=radius, catalog="PanSTARRS", data_release="dr2")
+        result = Catalogs.query_region(coord, radius=radius, catalog="PanSTARRS", data_release="dr1")
     except ValueError as e:
         # The MAST API occasionally returns 'None' strings in float columns, which astroquery
         # fails to convert. Treat this as no result found.

lightcurver/utilities/footprint.py

@@ -22,7 +22,7 @@ def get_combined_footprint_hash(user_config, frames_id_list):
     """
     if user_config['star_selection_strategy'] != 'ROI_disk':
         # then it depends on the frames we're considering.
-        frames_hash = get_frames_hash(frames_id_list)
+        return get_frames_hash(frames_id_list)
     else:
         # if ROI_disk, it does not depend on the frames: unique region defined by its radius.
         return hash(user_config['ROI_disk_radius_arcseconds'])