Code refactor, add new FILTER kw, fix small bugs (#163)

iop4admin/modeladmins/astrosource.py

@@ -239,7 +239,8 @@ def remove_field_stars_from_panstarrs(self, request, queryset):
         main_src = queryset.first()
 
         field_stars = main_src.calibrators.filter(name__startswith="PanSTARRS")
+        n_field_stars = len(field_stars) # count before deletion
         field_stars.delete()
-        logger.info(f"Removed {len(field_stars)} PanSTARRS field stars for {main_src.name}")
+        logger.info(f"Removed {n_field_stars} PanSTARRS field stars for {main_src.name}")
 
-        messages.success(request, f"Removed {len(field_stars)} PanSTARRS field stars for {main_src.name}")
+        messages.success(request, f"Removed {n_field_stars} PanSTARRS field stars for {main_src.name}")

iop4admin/modeladmins/reducedfit.py

@@ -64,7 +64,7 @@ def exptime(self, obj):
     @admin.display(description='SRCS IN FIELD')
     def get_targets_in_field(self, obj):
 
-        cat_targets = list(obj.sources_in_field.filter(Q(srctype=SRCTYPES.BLAZAR) | Q(srctype=SRCTYPES.STAR)).values_list('name', flat=True))
+        cat_targets = list(obj.sources_in_field.filter(is_calibrator=False).values_list('name', flat=True))
 
         if len(cat_targets) > 0:
             return cat_targets
@@ -73,6 +73,8 @@ def get_targets_in_field(self, obj):
             kw_obj_val = obj.rawfit.header['OBJECT']
         except FileNotFoundError:
             return format_html(f"<i>rawfit not found</i>")
+        except KeyError:
+            return format_html(f"<i>no OBJECT keyword</i>")
 
         guessed_target = AstroSource.objects.filter(Q(name__icontains=kw_obj_val) | Q(other_names__icontains=kw_obj_val)).values_list('name', flat=True)
 

iop4lib/instruments/cafos.py

@@ -15,6 +15,7 @@
 from iop4lib.enums import *
 from .instrument import Instrument
 from iop4lib.telescopes import CAHAT220
+from iop4lib.utils import filter_zero_points, calibrate_photopolresult
 
 # logging
 import logging
@@ -382,8 +383,7 @@ def compute_relative_polarimetry(cls, polarimetry_group):
                     mag_zp_err = np.nan
                 else:
                     mag_zp = mag_known - mag_inst
-                    # mag_zp_err = math.sqrt(mag_known_err ** 2 + mag_inst_err ** 2)
-                    mag_zp_err = math.fabs(mag_inst_err) # do not add error on literature magnitude
+                    mag_zp_err = math.sqrt(mag_known_err ** 2 + mag_inst_err ** 2)
             else:
                 mag_zp = None
                 mag_zp_err = None
@@ -408,37 +408,16 @@ def compute_relative_polarimetry(cls, polarimetry_group):
             logger.error("No results could be computed for this group.")
             return
 
-        # 3. Get average zero point from zp of all calibrators in the group
-
-        calib_mag_zp_array = np.array([result.mag_zp or np.nan for result in photopolresult_L if result.astrosource.is_calibrator]) # else it fills with None also and the dtype becomes object
-        calib_mag_zp_array = calib_mag_zp_array[~np.isnan(calib_mag_zp_array)]
-
-        calib_mag_zp_array_err = np.array([result.mag_zp_err or np.nan for result in photopolresult_L if result.astrosource.is_calibrator])
-        calib_mag_zp_array_err = calib_mag_zp_array_err[~np.isnan(calib_mag_zp_array_err)]
-
-        if len(calib_mag_zp_array) == 0:
-            logger.error(f"Can not compute magnitude during relative photo-polarimetry without any calibrators for this reduced fit.")
-
-        zp_avg = np.nanmean(calib_mag_zp_array)
-        zp_std = np.nanstd(calib_mag_zp_array)
-
-        zp_err = np.sqrt(np.nansum(calib_mag_zp_array_err ** 2)) / len(calib_mag_zp_array_err)
-        zp_err = math.sqrt(zp_err ** 2 + zp_std ** 2)
-
-        # 4. Compute the calibrated magnitudes for non-calibrators in the group using the averaged zero point
+        # 3. Compute the calibrated magnitudes for non-calibrators in the group using the averaged zero point
 
         for result in photopolresult_L:
 
             if result.astrosource.is_calibrator:
                 continue
 
-            result.mag_zp = zp_avg
-            result.mag_zp_err = zp_err
-
-            result.mag = result.mag_inst + zp_avg
-            result.mag_err = math.sqrt(result.mag_inst_err ** 2 + zp_err ** 2)
+            logger.debug(f"calibrating {result}")
 
-            result.save()
+            calibrate_photopolresult(result, photopolresult_L)
 
         # 5. Save results
         for result in photopolresult_L:

iop4lib/instruments/dipol.py

@@ -254,7 +254,7 @@ def apply_masters(cls, reducedfit):
 
 
     @classmethod
-    def get_header_hintobject(self, rawfit):
+    def get_header_hintobject(self, rawfit: 'RawFit'):
         r""" Overriden for DIPOL, which are using the convention for the other_names field. 
         
         The regex used has been obtained from the notebook checking all keywords.
@@ -431,7 +431,7 @@ def build_wcs(cls, reducedfit: 'ReducedFit', summary_kwargs : dict = None, metho
             redf_phot = ReducedFit.objects.filter(instrument=reducedfit.instrument,
                                                   sources_in_field__in=[reducedfit.header_hintobject], 
                                                   obsmode=OBSMODES.PHOTOMETRY, 
-                                                  flags__has=ReducedFit.FLAGS.BUILT_REDUCED).first()
+                                                  flags__has=ReducedFit.FLAGS.BUILT_REDUCED).order_by('-juliandate').first()
             try:
                 n_expected_simbad_sources = len(get_simbad_sources(reducedfit.header_hintobject.coord, radius=(reducedfit.width*cls.arcsec_per_pix*u.arcsec)))
             except Exception as e:
@@ -587,7 +587,7 @@ def _build_wcs_for_polarimetry_images_photo_quads(cls, redf: 'ReducedFit', summa
         redf_phot = ReducedFit.objects.filter(instrument=redf_pol.instrument, 
                                               sources_in_field__in=[target_src], 
                                               obsmode=OBSMODES.PHOTOMETRY, 
-                                              flags__has=ReducedFit.FLAGS.BUILT_REDUCED).first()
+                                              flags__has=ReducedFit.FLAGS.BUILT_REDUCED).order_by('-juliandate').first()
 
         if redf_phot is None:
             logger.error(f"No astro-calibrated photometry field found for {redf_pol}.")

iop4lib/instruments/instrument.py

@@ -21,6 +21,7 @@
 
 # iop4lib imports
 from iop4lib.enums import *
+from iop4lib.utils import filter_zero_points, calibrate_photopolresult
 
 # logging
 import logging
@@ -684,40 +685,9 @@ def compute_relative_photometry(cls, redf: 'ReducedFit') -> None:
             if result.astrosource.is_calibrator:
                 continue
 
-            # 3.a Average the zero points
+            logger.debug(f"{redf}: calibrating {result}")
 
-            # get all the computed photopolresults that calibrate this source
-            calibrator_results_L = [r for r in photopolresult_L if r.astrosource.calibrates.filter(pk=result.astrosource.pk).exists()]
-
-            # Create an array with nan instead of None (this avoids the dtype becoming object)
-            calib_mag_zp_array = np.array([result.mag_zp or np.nan for result in calibrator_results_L if result.astrosource.is_calibrator]) 
-            calib_mag_zp_array = calib_mag_zp_array[~np.isnan(calib_mag_zp_array)]
-
-            calib_mag_zp_array_err = np.array([result.mag_zp_err or np.nan for result in calibrator_results_L if result.astrosource.is_calibrator])
-            calib_mag_zp_array_err = calib_mag_zp_array_err[~np.isnan(calib_mag_zp_array_err)]
-
-            if len(calib_mag_zp_array) == 0:
-                logger.error(f"{redf}: can not perform relative photometry on source {result.astrosource.name}, no calibrator zero-points found.")
-                # [result.delete() for result in redf.photopolresults.all()]
-                continue
-
-            zp_avg = np.nanmean(calib_mag_zp_array)
-            zp_std = np.nanstd(calib_mag_zp_array)
-
-            zp_err = math.sqrt(np.sum(calib_mag_zp_array_err**2)) / len(calib_mag_zp_array_err)
-            zp_err = math.sqrt(zp_std**2 + zp_err**2)
-
-            # 3.b Compute the calibrated magnitude
-
-            # save the zp (to be) used
-            result.mag_zp = zp_avg
-            result.mag_zp_err = zp_err
-
-            # compute the calibrated magnitude
-            result.mag = zp_avg + result.mag_inst
-            result.mag_err = math.sqrt(result.mag_inst_err**2 + zp_err**2)
-
-            result.save()
+            calibrate_photopolresult(result, photopolresult_L)
 
         # 5. Save the results
 
@@ -731,7 +701,7 @@ def compute_relative_photometry(cls, redf: 'ReducedFit') -> None:
 
 
     @classmethod
-    def estimate_common_apertures(cls, reducedfits, reductionmethod=None, fit_boxsize=None, search_boxsize=(90,90), fwhm_min=2, fwhm_max=50):
+    def estimate_common_apertures(cls, reducedfits, reductionmethod=None, fit_boxsize=None, search_boxsize=(90,90), fwhm_min=2, fwhm_max=50, fwhm_default=3.5):
         r"""estimate an appropriate common aperture for a list of reduced fits.
         
         It fits the target source profile in the fields and returns some multiples of the fwhm which are used as the aperture and as the inner and outer radius of the annulus for local bkg estimation).
@@ -771,11 +741,10 @@ def estimate_common_apertures(cls, reducedfits, reductionmethod=None, fit_boxsiz
             mean_fwhm = np.mean(fwhm_L)
         else:
             logger.error(f"Could not find an appropriate aperture for Reduced Fits {[redf.id for redf in reducedfits]}, using standard fwhm of 3.5px")
-            mean_fwhm = 3.5
+            mean_fwhm = fwhm_default
 
         sigma = mean_fwhm / (2*np.sqrt(2*math.log(2)))
-        r = sigma
 
-        return 6.0*r, 7.0*r, 15.0*r, {'mean_fwhm':mean_fwhm, 'sigma':sigma}
+        return 3.0*sigma, 5.0*sigma, 9.0*sigma, {'mean_fwhm':mean_fwhm, 'sigma':sigma}
 
 

iop4lib/instruments/osn_cameras.py

@@ -16,6 +16,7 @@
 from iop4lib.enums import *
 from .instrument import Instrument
 from iop4lib.telescopes import OSNT090, OSNT150
+from iop4lib.utils import filter_zero_points, calibrate_photopolresult
 
 # logging
 import logging
@@ -114,7 +115,7 @@ def classify_obsmode_rawfit(cls, rawfit):
             raise ValueError("Cannot classify obsmode if band is ERROR.")
 
         if rawfit.band == BANDS.R:
-            if rawfit.header['FILTER'] == "R":
+            if rawfit.header['FILTER'] == "R" or rawfit.header['FILTER'] == "R Johnson":
                 rawfit.obsmode = OBSMODES.PHOTOMETRY
             else:
                 logger.debug("Band is R, but FILTER is not exactly R, for OSN this must mean it is polarimetry. Trying to extract angle from FILTER keyword.")
@@ -407,37 +408,16 @@ def compute_relative_polarimetry(cls, polarimetry_group):
             photopolresult_L.append(result)
 
 
-        # 3. Get average zero point from zp of all calibrators in the group
-
-        calib_mag_zp_array = np.array([result.mag_zp or np.nan for result in photopolresult_L if result.astrosource.is_calibrator]) # else it fills with None also and the dtype becomes object
-        calib_mag_zp_array = calib_mag_zp_array[~np.isnan(calib_mag_zp_array)]
-
-        calib_mag_zp_array_err = np.array([result.mag_zp_err or np.nan for result in photopolresult_L if result.astrosource.is_calibrator])
-        calib_mag_zp_array_err = calib_mag_zp_array_err[~np.isnan(calib_mag_zp_array_err)]
-
-        if len(calib_mag_zp_array) == 0:
-            logger.error(f"Can not compute magnitude during relative photo-polarimetry without any calibrators for this reduced fit.")
-
-        zp_avg = np.nanmean(calib_mag_zp_array)
-        zp_std = np.nanstd(calib_mag_zp_array)
-
-        zp_err = np.sqrt(np.nansum(calib_mag_zp_array_err ** 2)) / len(calib_mag_zp_array_err)
-        zp_err = math.sqrt(zp_err ** 2 + zp_std ** 2)
-
-        # 4. Compute the calibrated magnitudes for non-calibrators in the group using the averaged zero point
+        # 3. Compute the calibrated magnitudes for non-calibrators in the group using the averaged zero point
 
         for result in photopolresult_L:
 
             if result.astrosource.is_calibrator:
                 continue
 
-            result.mag_zp = zp_avg
-            result.mag_zp_err = zp_err
-
-            result.mag = result.mag_inst + zp_avg
-            result.mag_err = math.sqrt(result.mag_inst_err ** 2 + zp_err ** 2)
+            logger.debug(f"calibrating {result}")
 
-            result.save()
+            calibrate_photopolresult(result, photopolresult_L)
 
         # 5. Save results
         for result in photopolresult_L:

iop4lib/iop4.py

@@ -104,16 +104,20 @@ def process_astrosource(args):
     from iop4lib.db import ReducedFit, AstroSource, PhotoPolResult, Epoch
 
     astrosource = AstroSource.objects.get(name=args.astrosource)
+
     qs_all = ReducedFit.objects.filter(epoch__night__gte=args.date_start, epoch__night__lte=args.date_end)
 
     logger.info(f"Found {qs_all.count()} reduced fits between {args.date_start} and {args.date_end}.")
 
-    # filter files that have identified this source (in sources_in_field) or have this source as header_hintobject
-
-    redfL = list()
-    for redf in qs_all:
-        if redf.sources_in_field.filter(name=args.astrosource).exists() or redf.header_hintobject == astrosource:
-            redfL.append(redf)
+    if args.only_sources_in_field:
+        qs_all = qs_all.filter(sources_in_field__name=args.astrosource)
+        redfL = list(qs_all)
+    else:
+        # filter files that have identified this source (in sources_in_field) or have this source as header_hintobject
+        redfL = list()
+        for redf in qs_all:
+            if redf.sources_in_field.filter(name=args.astrosource).exists() or redf.header_hintobject == astrosource:
+                redfL.append(redf)
 
     qs_redf = ReducedFit.objects.filter(pk__in=[redf.pk for redf in redfL])
 
@@ -128,6 +132,12 @@ def process_astrosource(args):
         redfL_failed = [redf for redf in redfL if redf.has_flag(ReducedFit.FLAGS.ERROR_ASTROMETRY)]
         Epoch.reduce_reducedfits(redfL_failed)
 
+    if args.check_sources_in_field_again:
+            for reducedfit in redfL:
+                sources_in_field = AstroSource.get_sources_in_field(fit=reducedfit)    
+                logger.debug(f"{reducedfit}: found {len(sources_in_field)} sources in field.")        
+                reducedfit.sources_in_field.set(sources_in_field, clear=True)
+
     if args.recompute:
         qs_res = PhotoPolResult.objects.filter(reducedfits__in=redfL)
         n_res_before = qs_res.count()
@@ -359,6 +369,8 @@ def main():
     parser.add_argument('--astrosource', type=str, default=None, help='<Optional> Select files only of this source')
     parser.add_argument('--recompute', action='store_true', help='<Optional> Recompute photometry and polarimetry results')
     parser.add_argument('--retry-failed', action='store_true', help='<Optional> Retry failed reduced fits')
+    parser.add_argument('--only-sources-in-field', action='store_true', help='<Optional> Only process files that have this source in sources_in_field')
+    parser.add_argument('--check-sources-in-field-again', action='store_true', help='<Optional> Check sources_in_field again')
 
     # other options
     parser.add_argument('--skip-remote-file-list', action='store_true', help='<Optional> Skip remote file list check')

iop4lib/iop4_night_summary.html

@@ -93,6 +93,15 @@ <h2>Overview</h2>
 
             <h2>Summary of results (band R only)</h2>
 
+            {% if sources_without_calibrators %}
+                <p style="color:red;">
+                    Some sources do not have calibrators:
+                    {% for source in sources_without_calibrators %}
+                        {{ source.name }}{% if not forloop.last %}, {% endif %}
+                    {% endfor %}
+                </p>
+            {% endif %}
+
             <p>Some of these plots might include the previous existing night for the source for comparison.</p>
 
             {% for srcname, imgb64 in results_summary_images.items %}

iop4lib/iop4_night_summary.py

@@ -72,6 +72,9 @@ def gather_context(args):
 
     sources = AstroSource.objects.exclude(is_calibrator=True).filter(photopolresults__epoch__night=args.date).distinct()
 
+    # check if some of these sources lack calibrators
+    sources_without_calibrators = [source for source in sources if not source.calibrators.exists()]
+
     if sources:
         instruments = [instrument.name for instrument in Instrument.get_known()]
         colors = [mplt.colormaps['tab10'](i) for i in range(len(instruments))]
@@ -160,6 +163,7 @@ def gather_context(args):
     context['epochs'] = epochs    
     context['sources'] = sources
     context['results_summary_images'] = results_summary_images
+    context['sources_without_calibrators'] = sources_without_calibrators
     context['args'] = args
 
     return context