mask filters related to observability, source brightness, galactic latitude implemented

Author: damien turpin

fink_filters/filter_svom_ft_tom/filter.py

@@ -0,0 +1,345 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Nov 21 14:37:49 2022
+
+@author: dt270490
+"""
+
+import numpy as np
+import healpy as hp
+from astroplan import Observer, FixedTarget
+from astroplan import is_observable, is_always_observable, months_observable
+from astroplan import AltitudeConstraint, AirmassConstraint, AtNightConstraint,\
+MoonSeparationConstraint
+from astropy.time import Time
+from astropy.coordinates import SkyCoord
+import astropy.units as u
+from astropy.table import Table
+import numpy as np
+import matplotlib.pyplot as plt
+from mhealpy import HealpixMap
+import mhealpy as hmap
+import pandas as pd
+import json
+import requests
+import io
+
+def get_test_alert():
+    
+    # Get latests 5 Early SN Ia candidates
+    r = requests.post(
+      'https://fink-portal.org/api/v1/latests',
+      json={
+        'class': 'Early SN Ia candidate',
+        'n': '2'
+      }
+    )
+    
+    # Format output in a DataFrame
+    pdf = pd.read_json(io.BytesIO(r.content))
+    return pdf
+
+def load_observatories(obs_path_file: str):
+    """
+    Load the observatories for which we want to crossmacth the observable 
+    skies
+
+    Parameters
+    ----------
+    obs_path_file : str
+        Path to the files where the observatory settings are stored.
+
+    Returns
+    -------
+    observatories : Dataframes
+        pandas Dataframe of the observatory configurations 
+        to feed the astroplan.Observer package.
+
+    """
+    observatories = pd.read_csv(obs_path_file, delimiter=";")
+    return observatories    
+    
+def simu_night_time_interval(ref_obs, ref_date: str, n_days: int, day_bin: int):
+    """
+    Compute the list of the date interval during which the observable skies
+    will be crossmatched
+
+    Parameters
+    ----------
+    ref_obs : astroplan.observer.Observer
+        astroplan.observer.Observer object for the Observatory chosen as the 
+        reference observatory.
+    ref_date : str
+        Date of reference to start the simulation.
+    n_days : int
+        Number of simulated days.
+    day_bin : int
+        Day interval between two simulated nights.
+
+    Returns
+    -------
+    None.
+
+    """  
+    # compute the start time of the nearest (compared to the detection time)
+    # night
+    date_start_night = ref_obs.twilight_evening_astronomical(ref_date,
+                                                                which='nearest')
+    # compute the morning time of the nearest (compared to the detection time)
+    # day
+    date_end_night = ref_obs.twilight_morning_astronomical(date_start_night,
+                                                                which='nearest')
+    # If the detection time is within the current night, use the detection time
+    # as a starting date and the end of the night as an ending date
+    if date_start_night < ref_date and date_end_night>ref_date:
+        time_ranges = Time([ref_date.iso, date_end_night.iso])
+    # If the detection time is during day time, use the night starting date
+    # as a starting date and the end of the night as an ending date
+    elif date_start_night > ref_date and date_end_night>ref_date:
+        time_ranges = Time([date_start_night.iso, date_end_night.iso])
+    # If the detection time is after the nearest night, use the next night 
+    # starting date as a starting date and the end of the night as an
+    # ending date
+    elif date_start_night < ref_date and date_end_night<ref_date:
+        date_start_night = ref_obs.twilight_evening_astronomical(ref_date,
+                                                                    which='next')
+        date_end_night = ref_obs.twilight_morning_astronomical(date_start_night,
+                                                                    which='next')
+        time_ranges = Time([date_start_night.iso, date_end_night.iso])
+    return time_ranges
+
+
+def build_targets(ras: np.array, decs: np.array):
+    """
+    Make the target objects to be ingested by the astroplan is_observable
+    functions
+
+    Parameters
+    ----------
+    ras : np.array
+        Numpy arrays of right ascencion in units of degrees.
+    decs : np.array
+        Numpy arrays of declinations in units of degrees.
+
+    Returns
+    -------
+    targets : 
+        List of targets for which we want to estimate the observability
+
+    """
+    ras_grid, decs_grid = np.meshgrid(ras, decs)
+
+    target_table = Table()
+    target_table["ra"] = np.reshape(ras_grid, ras_grid.shape[0] * ras_grid.shape[1])
+    target_table["dec"] = np.reshape(decs_grid, decs_grid.shape[0] * decs_grid.shape[1])
+
+    targets = FixedTarget(
+        coord=SkyCoord(ra=target_table["ra"] * u.deg, dec=target_table["dec"] * u.deg)
+    )
+
+    return targets, target_table
+
+
+def make_visibility_masks(constraints, observatory, targets, time_range: list):
+    """
+    Build the skymap mask regarding the obsverational constraints:  
+    True = the position is visible
+    False = the position is not visible
+
+    Parameters
+    ----------
+    constraints : TYPE
+        DESCRIPTION.
+    observatory : TYPE
+        DESCRIPTION.
+    targets : TYPE
+        DESCRIPTION.
+    time_range : list
+        time range during which the observability is estimated
+
+    Returns
+    -------
+    mask_visibility : np.array
+        Numpy arrays of booleans
+
+    """
+
+    mask_visibility = is_observable(
+        constraints, observatory, targets, time_range=time_range
+    )
+
+    return mask_visibility[0]
+
+def obs_visibility_pdf(ras:list,decs:list,det_time:list,
+                       observatories:pd.DataFrame,constraints:list):
+    """
+    This function builds a pandas DataFrame that collect all the visibility
+    flags of the observatories in the network related to the positions of 
+    optical transients and based on pre-defined observational constraints
+
+    Parameters
+    ----------
+    ras : list
+        list of right ascencion given in degrees.
+    decs : list
+        list of declination given in degrees.
+    det_time : list
+        list of first detection times expressed in Julian Date.
+    observatories : pd.DataFrame
+        Dataframe gathering the observatories informations.
+    constraints : list
+        list of astronomical constraints to be satisfied by each transient
+        candidates.
+
+    Returns
+    -------
+    df : pd.DataFrame
+        The Dataframe related to the observatories visibility for each candid.
+
+    """
+    mask_visibilities = []
+    for i in range(len(observatories)):
+        try:
+            observatory = Observer.at_site(observatories.obs_name[i])
+        except:
+            observatory = Observer(
+                longitude=observatories.longitude[i] * u.deg,
+                latitude=observatories.latitude[i] * u.deg,
+                elevation=observatories.elevation[i] * u.m,
+                name=observatories.obs_name[i],
+            )
+        mask_visibility = []
+        for k in range(len(pdf_test)):
+            targets = FixedTarget(
+                coord=SkyCoord(ra=pdf_test['i:ra'][k] * u.deg,
+                               dec=pdf_test['i:dec'][k]* u.deg)
+            )
+            time_ranges = simu_night_time_interval(
+                    observatory, Time(pdf_test['i:jd'][k],format='jd'), 1, 1)
+        
+            mask_visibility.append(make_visibility_masks(
+                        constraints, observatory, targets, time_ranges
+                    ))
+        mask_visibilities.append(mask_visibility)
+    # Build the DataFrame     
+    col_name = ['objectId','candid']+observatories.obs_name.values.tolist()
+    df = pd.DataFrame(columns=col_name)
+    df.objectId = pdf_test['i:objectId']
+    df.candid = pdf_test['i:candid']
+    
+    for i in range(len(observatories.obs_name)):
+       df[observatories.obs_name[i]] = mask_visibilities[i]
+    return df
+
+def brightness_mask(mag:list,mag_err:list,fid:list,mag_threshold:int):
+    """
+    Create a mask related to the brightness of the optical transient
+    candidates
+
+    Parameters
+    ----------
+    mag : list
+        list of the transient magnitudes.
+    mag_err : list
+        list of the transient magnitude errors.
+    fid : list
+        list of the filters used for the alert observation.
+    mag_threshold : int
+        threshold on the magnitude to keep a transient for follow-up.
+
+    Returns
+    -------
+    mask : pd.DataFrame
+        DataFrame masking the too faint candidates.
+
+    """
+    mask = pdf_test['i:magpsf']<mag_threshold
+    return mask.tolist()
+
+def extragal_mask(ras:list,decs:list,gal_lat_threshold:float):
+    """
+    Create a mask related to the galactic latitude of the optical transient
+    candidates
+
+    Parameters
+    ----------
+    ras : list
+        list of right ascencion given in degrees.
+    decs : list
+        list of declination given in degrees.
+    gal_lat_threshold : float
+        threshold on the galactic latitude to keep a transient for follow-up.
+
+    Returns
+    -------
+    mask : list
+        DataFrame masking the candidates too close to the galactic plane.
+
+    """
+    coord=SkyCoord(ra=ras * u.deg, dec=decs * u.deg)
+    mask = (coord.galactic.b.degree > gal_lat_threshold) |\
+        (coord.galactic.b.degree < -1*gal_lat_threshold)
+    return mask.tolist()
+
+def net_visibity_mask(pdf_cand_visibility:pd.DataFrame,
+                      mask_brightness:list,
+                      mask_extragal:list):
+    """
+    Build the final mask to filter out the transient candidates
+
+    Parameters
+    ----------
+    pdf_cand_visibility : pd.DataFrame
+        The Dataframe related to the observatories visibility for each candid.
+    mask_brightness : list
+        DataFrame masking the too faint candidates.
+    mask_extragal : list
+        DataFrame masking the candidates too close to the galactic plane.
+
+    Returns
+    -------
+    mask_final : list
+        the final mask to filter out the transient candidates
+
+    """
+    # The source has to visible at San Pedro Martir and Xinglong and NOT
+    mask_net_visibility = (pdf_cand_visibility['OAN-SPM'] == True) & \
+        (pdf_cand_visibility['Xinglong'] == True) &\
+            (pdf_cand_visibility['ORM'] == True)
+    
+    mask_final = mask_net_visibility & mask_brightness & mask_extragal
+    return mask_final.tolist()
+
+if __name__ == "__main__":
+    # Load the observatories info
+    obs_path_file = "../../conf/"
+    obs_filename = "svom_network.csv"
+    observatories = load_observatories(obs_path_file + obs_filename)
+    # Define the observationnal constraints
+    constraints = [
+        AltitudeConstraint(40 * u.deg, 80 * u.deg),
+        AirmassConstraint(2), MoonSeparationConstraint(min=20* u.deg)]#, AtNightConstraint.twilight_astronomical()]
+    #Load a sample of Fink alert tests
+    pdf_test = get_test_alert()
+    # Build Dataframe related to the observatories visibility for each candid
+    pdf_cand_visibility = obs_visibility_pdf(pdf_test['i:ra'],
+                                             pdf_test['i:dec'],
+                                             pdf_test['i:jd'],
+                                             observatories,
+                                             constraints)
+    # Add a column in the pdf_cand_visibility DataFrame to flag if 
+    # the source is bright enough for Colibri and NOT
+    mask_brightness = brightness_mask(pdf_test['i:magpsf'],
+                                      pdf_test['i:sigmapsf'],pdf_test['i:fid'],
+                                      18)
+    pdf_cand_visibility['mag_mask'] = mask_brightness
+    # Add a column in the pdf_cand_visibility DataFrame to flag if 
+    # the source is at high galactic latitutes
+    mask_extragal = extragal_mask(pdf_test['i:ra'],pdf_test['i:dec'],10)
+    pdf_cand_visibility['extragal_coord_mask'] = mask_extragal
+    # Make the final mask on every ZTF candidates
+    mask_final = net_visibity_mask(pdf_cand_visibility,mask_brightness,
+                                   mask_extragal)
+    #Add the final mask to the initial alert DataFrame
+    pdf_test['svom_mask'] = mask_final