Continuum test

arccnet/data_generation/timeseries/drms_pipeline.py

@@ -29,9 +29,9 @@
 
     # Logging settings here.
     drms_log = logging.getLogger("drms")
-    drms_log.setLevel("WARNING")
+    drms_log.setLevel("ERROR")
     reproj_log = logging.getLogger("reproject.common")
-    reproj_log.setLevel("WARNING")
+    reproj_log.setLevel("ERROR")
     # May need to find a more robust solution with filters/exceptions for this.
     astropy_log.setLevel("ERROR")
     data_path = config["paths"]["data_folder"]
@@ -43,7 +43,8 @@
         duration=6,
         long_lim=65,
         types=["F1", "F2", "N1", "N2"],
-    )
+    )[0]
+
     cores = int(config["drms"]["cores"])
     with ProcessPoolExecutor(cores) as executor:
         for record in starts:
@@ -54,7 +55,9 @@
             patch_height = int(config["drms"]["patch_height"]) * u.pix
             patch_width = int(config["drms"]["patch_width"]) * u.pix
             try:
-                print(record["noaa_number"], record["goes_class"], record["start_time"], record["category"])
+                logging.info(
+                    f"{record['noaa_number']} {record['goes_class']} {record['start_time']} {record['category']}"
+                )
                 aia_maps, hmi_maps = drms_pipeline(
                     start_t=start,
                     end_t=end,
@@ -64,6 +67,9 @@
                     wavelengths=config["drms"]["wavelengths"],
                     sample=config["drms"]["sample"],
                 )
+                if len(aia_maps != 60):
+                    logging.warning("Bad Run detected - missing frames.")
+                    continue
 
                 hmi_proc = list(
                     tqdm(

arccnet/data_generation/timeseries/sdo_processing.py

@@ -1,4 +1,5 @@
 import os
+import re
 import sys
 import glob
 import shutil
@@ -29,9 +30,12 @@
 from arccnet.data_generation.utils.utils import save_compressed_map
 from arccnet.visualisation.data import mosaic_animate, mosaic_plot
 
+data_path = config["paths"]["data_folder"]
+qry_log = "bad_queries"
+
 warnings.simplefilter("ignore", RuntimeWarning)
 reproj_log = logging.getLogger("reproject.common")
-reproj_log.setLevel("WARNING")
+reproj_log.setLevel("ERROR")
 os.environ["JSOC_EMAIL"] = "[email protected]"
 data_path = config["paths"]["data_folder"]
 
@@ -51,6 +55,30 @@
 ]
 
 
+def bad_query(qry, data_path, name):
+    r"""
+    Logs bad drms queries to document for future reference.
+
+    Parameters
+    ----------
+        qry: `str`
+            Drms query which returns an empty dataframe.
+        data_path : `str`
+            Path to arccnet level 04 data.
+        name: `str`
+            Log filename.
+    """
+    logging.warning(f"Bad Query Detected - {qry}")
+    f_name = f"{data_path}/logs/{name}.txt"
+    if not os.path.exists(f_name):
+        file = open(f"{data_path}/logs/{name}.txt", "x")
+    file = open(f"{data_path}/logs/{name}.txt", "+")
+    entries = [row for row in file]
+    if qry not in entries:
+        file.write(qry)
+    file.close()
+
+
 def rand_select(table, size, types: list):
     r"""
     Randomly selects targets from provided table and list of data subsets.
@@ -127,7 +155,7 @@ def read_data(hek_path: str, srs_path: str, size: int, duration: int, long_lim:
     srs["srs_date"] = srs["target_time"].value
     srs["srs_date"] = [date.split("T")[0] for date in srs["srs_date"]]
     srs["srs_end_time"] = [(Time(time) + duration * u.hour) for time in srs["target_time"]]
-    print("Parsing Flares")
+    logging.info("Parsing Flares")
     flares["start_time"] = [time - (duration + 1) * u.hour for time in flares["start_time"]]
     flares["start_time"].format = "fits"
     flares["tb_date"] = flares["start_time"].value
@@ -146,7 +174,8 @@ def read_data(hek_path: str, srs_path: str, size: int, duration: int, long_lim:
         SkyCoord(lon * u.deg, lat * u.deg, obstime=t_time, observer="earth", frame=frames.HeliographicStonyhurst)
         for lat, lon, t_time in zip(srs["latitude"], srs["longitude"], srs["target_time"])
     ]
-    print("Parsing Active Regions")
+
+    logging.info("Parsing Active Regions")
     ar_cat, fl_cat = [], []
     for ar in srs:
         erl_time = Time(ar["target_time"]) - (duration + 1) * u.hour
@@ -295,14 +324,14 @@ def drms_pipeline(
     aia_query, aia_export = aia_query_export(hmi_query, aia_keys, wavelengths)
 
     hmi_dls, hmi_exs = l1_file_save(hmi_export, hmi_query, path)
-    cnt_dls, cnt_exs = l1_file_save(ic_query, ic_export, path)
+    cnt_dls, cnt_exs = l1_file_save(ic_export, ic_query, path)
     aia_dls, aia_exs = l1_file_save(aia_export, aia_query, path)
     img_exs = list(itertools.chain(aia_exs, cnt_exs))
 
     hmi_maps = sunpy.map.Map(hmi_exs)
     hmi_maps = add_fnames(hmi_maps, hmi_exs)
     img_maps = sunpy.map.Map(img_exs)
-    img_maps = add_fnames(img_maps, aia_exs)
+    img_maps = add_fnames(img_maps, img_exs)
     return img_maps, hmi_maps
 
 
@@ -369,10 +398,8 @@ def hmi_continuum_export(hmi_query, keys):
     """
     client = drms.Client()
     qstrs_ic = [f"hmi.Ic_noLimbDark_720s[{time}]{{continuum}}" for time in hmi_query["T_REC"]]
-
-    ic_value = [hmi_rec_find(qstr, keys, 3, 12) for qstr in qstrs_ic]
-    unpacked_ic = list(itertools.chain.from_iterable(ic_value))
-    joined_num = [str(num) for num in unpacked_ic]
+    ic_value = [hmi_rec_find(qstr, keys, 3, 12, cont=True) for qstr in qstrs_ic]
+    joined_num = [str(num) for num in ic_value]
     ic_num_str = str(joined_num).strip("[]")
     ic_comb_qstr = f"hmi.Ic_noLimbDark_720s[! recnum in ({ic_num_str}) !]{{continuum}}"
 
@@ -445,9 +472,21 @@ def hmi_rec_find(qstr, keys, retries, sample, cont=False):
     client = drms.Client()
     count = 0
     qry = client.query(ds=qstr, key=keys)
-    time = sunpy.time.parse_time(qry["T_REC"].values[0])
+    if qry.empty:
+        time = sunpy.time.parse_time(re.search(r"\[(.*?)\]", qstr).group(1))
+        qry["QUALITY"] = [10000]
+        logging.warning("Bad Query - HMI")
+        bad_query(qry, data_path, qry_log)
+
+    else:
+        time = sunpy.time.parse_time(qry["T_REC"].values[0])
     while qry["QUALITY"].values[0] != 0 and count <= retries:
         qry = client.query(ds=f"{series}[{time}]" + seg, key=keys)
+        if qry.empty:
+            time = sunpy.time.parse_time(re.search(r"\[(.*?)\]", qstr).group(1))
+            qry["QUALITY"] = [10000]
+            logging.warning("Bad Query - HMI")
+            bad_query(qry, data_path, qry_log)
         time = change_time(time, sample)
         count += 1
     return qry["*recnum*"].values[0]
@@ -483,6 +522,9 @@ def aia_rec_find(qstr, keys, retries, time_add):
             retry += 1
         if qry["QUALITY"].values[0] == 0:
             return qry["FSN"].values
+    else:
+        logging.warning("Bad Query - AIA")
+        bad_query(qry, data_path, qry_log)
 
 
 def l1_file_save(export, query, path):
@@ -804,8 +846,6 @@ def map_reproject(sdo_packed):
             The path location of the saved submap.
     """
     hmi_origin, sdo_path, ar_num, center = sdo_packed
-    lon = int(center.lon.value)
-    # sdo_map = sunpy.map.Map(sdo_packed.l2_map)
     sdo_map = sunpy.map.Map(sdo_path)
     with propagate_with_solar_surface():
         sdo_rpr = sdo_map.reproject_to(hmi_origin.wcs)
@@ -819,7 +859,7 @@ def map_reproject(sdo_packed):
     sdo_rpr.meta["wavelnth"] = sdo_map.meta["wavelnth"]
     sdo_rpr.meta["date-obs"] = sdo_map.meta["date-obs"]
     if sdo_rpr.dimensions[0].value < int(config["drms"]["patch_width"]):
-        sdo_rpr = pad_map(sdo_rpr, lon, config["drms"]["patch_width"])
+        sdo_rpr = pad_map(sdo_rpr, config["drms"]["patch_width"])
     save_compressed_map(sdo_rpr, fits_path, hdu_type=CompImageHDU, overwrite=True)
 
     return fits_path
@@ -843,11 +883,12 @@ def vid_match(table, name, path):
         output_file : `str`
             A string containing the path of the completed mosaic animation.
     """
-    # Check this carefully
     hmi_files = table["HMI files"].value
+    table["Wavelength"] = [int(wave) for wave in table["Wavelength"]]
     wvls = np.unique([table["Wavelength"].value])
+
     hmi_files = np.unique(hmi_files)
-    nrows, ncols = 4, 3  # define your subplot grid
+    nrows, ncols = 4, 3  # define subplot grid
     for file in range(len(hmi_files)):
         hmi = hmi_files[file]
         mosaic_plot(hmi, name, file, nrows, ncols, wvls, table, path)
@@ -912,6 +953,7 @@ def pad_map(map, targ_width):
             The new, padded submap.
     """
     x_dim = map.dimensions[0].value
+    targ_width = int(targ_width)
     diff = int(targ_width - x_dim)
     if map.center.Tx > 0:
         data = np.pad(map.data, ((0, 0), (diff, 0)), constant_values=np.nan)

arccnet/visualisation/data.py

@@ -522,16 +522,18 @@ def mosaic_plot(hmi, name, file, nrows, ncols, wvls, table, path):
             files = table[table["Wavelength"] == wv]
             files = files[files["HMI files"] == hmi]
             aia_files = files["AIA files"]
+            cmap = f"sdoaia{wv}"
+            if wv == 6173:
+                # keep yellow cmap
+                cmap = "sdoaia4500"
             try:
                 aia_map = Map(aia_files.value[0])
                 ax = fig.add_subplot(nrows, ncols, i + 1)
-                ax.imshow(np.sqrt(aia_map.data), cmap=f"sdoaia{wv}")
+                ax.imshow(np.sqrt(aia_map.data), cmap=cmap)
                 ax.text(0.05, 0.05, f"{wv} - {aia_map.date}", color="w", transform=ax.transAxes, fontsize=5)
             except IndexError:
-                aia_map = np.zeros(hmi_map.data.shape)
                 ax = fig.add_subplot(nrows, ncols, i + 1)
-                ax.imshow(np.sqrt(aia_map.data), cmap=f"sdoaia{wv}")
-                ax.text(0.05, 0.05, f"{wv} - MISSING", color="w", transform=ax.transAxes, fontsize=5)
+                ax.text(0.05, 0.05, f"{wv} - MISSING", color="black", transform=ax.transAxes, fontsize=5)
 
         else:
             ax = fig.add_subplot(nrows, ncols, i + 1)
@@ -581,6 +583,8 @@ def mosaic_animate(base_dir, name):
     for file in sorted_files:
         frame = cv2.imread(file)
         out.write(frame)
+        # Deletes frame after being added to animation to save data on unneeded frames.
+        os.remove(file)
     out.release()
     print(f"Video saved to {output_file}")
     return output_file