ra & dec in db (draft)

magicctapipe/scripts/lst1_magic/semi_automatic_scripts/database_production/set_ra_dec.py

@@ -0,0 +1,139 @@
+import pandas as pd
+import numpy as np
+import json
+from astropy.coordinates import SkyCoord
+
+
+def main():
+
+    """
+    Main function
+    """
+
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument(
+        "--begin-date",
+        "-b",
+        dest="begin",
+        type=int,
+        default=0,
+        help="First date to update database (YYYYMMDD)",
+    )
+    parser.add_argument(
+        "--end-date",
+        "-e",
+        dest="end",
+        type=int,
+        default=0,
+        help="End date to update database (YYYYMMDD)",
+    )
+    parser.add_argument(
+        "--dict-ra",
+        "-r",
+        dest="ra",
+        type=str,
+        default='./ra_dict.txt',
+        help="File with dictionary of RA coordinates for sources",
+    )
+    parser.add_argument(
+        "--dict-dec",
+        "-d",
+        dest="dec",
+        type=str,
+        default='./dec_dict.txt',
+        help="File with dictionary of Dec coordinates for sources",
+    )
+    parser.add_argument(
+        "--mc-dec",
+        "-m",
+        dest="dec_mc",
+        type=str,
+        default='./dec_mc.txt',
+        help="File with list of MC declinations",
+    )
+
+
+    args = parser.parse_args()
+    config_file = resource_file("database_config.yaml")
+
+    with open(
+        config_file, "rb"
+    ) as fc:  # "rb" mode opens the file in binary format for reading
+        config_dict = yaml.safe_load(fc)
+
+    LST_h5 = config_dict["database_paths"]["LST"]
+    LST_key = config_dict["database_keys"]["LST"]
+
+    df_LST = pd.read_hdf(
+        LST_h5,
+        key=LST_key,
+    ) 
+    if 'ra' not in df_LST:
+        df_LST['ra']=np.nan
+    if 'dec' not in df_LST:
+        df_LST['dec']=np.nan
+    if 'MC_dec' not in df_LST:
+        df_LST['MC_dec']=np.nan
+    df_LST_full=df_LST.copy(deep=True)
+    if args.begin != 0:
+        df_LST = df_LST[df_LST["DATE"].astype(int) >= args.begin]
+    if args.end != 0:
+        df_LST = df_LST[df_LST["DATE"].astype(int) <= args.end]
+
+    sources=np.unique(df_LST['source'])
+    with open(args.ra) as f: 
+        dict_ra = f.read()
+    with open(args.dec) as f: 
+        dict_dec = f.read()
+    with open(args.dec_mc) as f: 
+        mc_dec = f.read()
+    ra_dict = json.loads(dict_ra)
+    dec_dict = json.loads(dict_dec)
+    dec_mc=mc_dec.replace(',','').split('\n')
+    for src in sources:
+
+        try: 
+            coord=SkyCoord.from_name(src)
+            src_dec=coord.dec.degree
+            src_ra=coord.ra.degree
+        except astropy.coordinates.name_resolve.NameResolveError:
+            print (f"{i}: {src} not found in astropy. Looking to the dictionaries...")    
+            if src in ra_dict.keys():
+                src_ra=float(ra_dict.get(src))
+
+            else:
+                print (f"\t {i}: {src} RA not in the dictionaries. Please add it to the dictionaries") 
+                i+=1
+                continue  
+            if src in dec_dict.keys():
+                src_dec=float(dec_dict.get(src))
+            else:
+                print (f"\t {i}: {src} Dec not in the dictionaries. Please add it to the dictionaries") 
+                i+=1  
+                continue 
+            i+=1
+
+        df['dec'] = np.where(df['source'] == src, src_dec, df['dec'])
+        df['ra'] = np.where(df['source'] == src, src_ra, df['ra'])
+        df['MC_dec'] = np.where(df['source'] == src, float(dec_mc[np.argmin(np.abs(src_dec-dec_mc))]) , df['MC_dec'])
+    df_LST = pd.concat([df_LST, df_LST_full]).drop_duplicates(
+        subset="LST1_run", keep="first"
+    )
+    df_LST.to_hdf(
+        LST_h5,
+        key=LST_key,
+        mode="w",
+        min_itemsize={
+            "lstchain_versions": 20,
+            "last_lstchain_file": 90,
+            "processed_lstchain_file": 90,
+        },
+    )
+
+
+if __name__ == "__main__":
+    main()
+
+
+