Comparison of the original and new Chimera codes

Author: FionnStack

chimerapy/chimera.py

@@ -123,14 +123,25 @@ def filter_ch(mask: NDArray, map_obj: Map, min_area: Quantity["area"] = 1e4 * u.
     on_disk :
         Remove CHs that are not on the disk (above the limb)
     """
+    log.debug("Calculating area map")
     area_map, disk_mask = calculate_area_map(map_obj)
-
+    log.debug("Labeling regions")
     labeled_mask = measure.label(mask * disk_mask if on_disk else mask)
+    log.debug("Finding regions")
     regions = measure.regionprops(labeled_mask)
-    breakpoint()
+    log.debug("Filtering regions by pixel area size")
+
+    min_pixel_area = 100
+    filtered_regions_pixel_area = [region for region in regions if region.area >= min_pixel_area]
+
+    log.debug(f"Found {len(regions)} regions, {len(filtered_regions_pixel_area)} after pixel area filtering")
 
     filtered_regions = []
-    for region in regions:
+    log.debug(f"Found {len(filtered_regions_pixel_area)} regions")
+    #import ipdb; ipdb.set_trace()
+    #breakpoint()
+    for region in filtered_regions_pixel_area:
+        log.debug(f"Processing region {region.label} with area {region.area}")
         region_mask = labeled_mask == region.label
         contours = measure.find_contours(region_mask)
         if contours:
@@ -141,9 +152,6 @@ def filter_ch(mask: NDArray, map_obj: Map, min_area: Quantity["area"] = 1e4 * u.
             if region_surface_area >= min_area and not np.all(region_mask & disk_mask):
                 region.surface_area = region_surface_area
                 filtered_regions.append(region)
-            else:
-                log.debug(f"Removing CH region {region.label}")
-                labeled_mask[region_mask] = 0
 
     filtered_regions = sorted(filtered_regions, key=lambda region: region.surface_area, reverse=True)
 
@@ -156,7 +164,6 @@ def get_coronal_holes(filtered_regions, map_obj, labeled_mask):
     for region in filtered_regions:
         coords = region.coords
         world_coords = map_obj.pixel_to_world(coords[:, 1] * u.pix, coords[:, 0] * u.pix)
-        # heliographic_coords = world_coords.transform_to("heliographic_stonyhurst")
 
         wb = world_coords[np.nanargmax(world_coords.Tx)]
         eb = world_coords[np.nanargmin(world_coords.Tx)]
@@ -217,10 +224,13 @@ def map_threshold(im_map):
 
 
 def chimera(m171, m193, m211):
+    log.debug("Generating candidate mask")
     ch_mask = generate_candidate_mask(m171, m193, m211)
+    log.debug("Filtering regions")
     labeled_mask, filtered_regions = filter_ch(ch_mask, m171)
-
+    log.debug("Getting coronal holes")
     coronal_holes = get_coronal_holes(filtered_regions, m171, labeled_mask)
+    log.debug("Finished")
 
     for ch in coronal_holes:
         print(
@@ -237,7 +247,7 @@ def chimera(m171, m193, m211):
             f"N: {ch['nb'].transform_to('heliographic_stonyhurst').lat.value:.2f}, S:{ch['sb'].transform_to('heliographic_stonyhurst').lat.value:.2f}"
             f"N-S Extent = {ch['extent_lat']:.2f} °"
         )
-
+    return labeled_mask, ch_mask, coronal_holes
 
 def run_chimera(date):
     date = parse_time(date)

chimerapy/chimera_original.py

@@ -15,21 +15,26 @@
 
 import astropy.units as u
 from astropy import wcs
+from sunpy.map import Map
 from astropy.io import fits
 from astropy.modeling.models import Gaussian2D
 
+def imhmi(m171):
+    """Create a dummy HMI map for the original CHIMERA code."""
+    imhmi = Map(np.zeros_like(m171.data), m171.meta)
+    return imhmi
 
 def chimera_legacy(im171=None, im193=None, im211=None, imhmi=None):
     file_path = "./"
 
     if im171 is None:
-        im171 = glob.glob(file_path + "*171*.fts.gz")
-        im193 = glob.glob(file_path + "*193*.fts.gz")
-        im211 = glob.glob(file_path + "*211*.fts.gz")
-        imhmi = glob.glob(file_path + "*hmi*.fts.gz")
+        im171 = Map("C:\\Users\\aoife\\CHIMERApy\\chimerapy\\scripts\\downloaded_data\\aia.lev1.171A_2016_10_31T00_00_10.35Z.image_lev1.fits")
+        im193 = Map("C:\\Users\\aoife\\CHIMERApy\\chimerapy\\scripts\\downloaded_data\\aia.lev1.193A_2016_10_31T00_00_05.85Z.image_lev1.fits")
+        im211 = Map("C:\\Users\\aoife\\CHIMERApy\\chimerapy\\scripts\\downloaded_data\\aia.lev1.211A_2016_10_31T00_00_10.62Z.image_lev1.fits")
+        imhmi_map = imhmi(im171)
 
     circ, data, datb, datc, dattoarc, hedb, iarr, props, rs, slate, center, xgrid, ygrid = chimera(
-        im171, im193, im211, imhmi
+        im171, im193, im211, imhmi_map
     )
 
     # =====sets ident back to max value of iarr======
@@ -42,8 +47,8 @@ def chimera_legacy(im171=None, im193=None, im211=None, imhmi=None):
     return circ, data, datb, datc, dattoarc, hedb, iarr, props, rs, slate, center, xgrid, ygrid
 
 
-def chimera(im171, im193, im211, imhmi):
-    if im171 == [] or im193 == [] or im211 == [] or imhmi == []:
+def chimera(im171, im193, im211, imhmi_map):
+    if im171 == [] or im193 == [] or im211 == [] or imhmi_map == []:
         print("Not all required files present")
         sys.exit()
     # =====Reads in data and resizes images=====
@@ -61,8 +66,8 @@ def chimera(im171, im193, im211, imhmi):
     datc = fits.getdata(im211[0], ext=0) / (hedc["EXPTIME"])
     dn = RectBivariateSpline(x, x, datc, kx=1, ky=1)
     datc = dn(np.arange(0, 4096), np.arange(0, 4096))
-    hedm = fits.getheader(imhmi[0], hdu_number)
-    datm = fits.getdata(imhmi[0], ext=0)
+    hedm = fits.getheader(imhmi_map[0], hdu_number)
+    datm = fits.getdata(imhmi_map[0], ext=0)
     # dn = scipy.interpolate.interp2d(np.arange(4096), np.arange(4096), datm)
     # datm = dn(np.arange(0, 1024)*4, np.arange(0, 1024)*4)
     if hedm["crota1"] > 90:

chimerapy/scripts/comparison.py

@@ -0,0 +1,88 @@
+from sunpy.net import Fido, attrs as a
+import astropy.units as u
+import os
+from sunpy.map import Map
+import numpy as np
+import sys
+import chimerapy.chimera as chimera_new
+import chimerapy.chimera_original as chimera_old
+
+"""
+date = "2016/11/04"
+wavelength = 193 * u.Angstrom
+
+script_dir = os.path.dirname(os.path.abspath(__file__))
+save_path = os.path.join(script_dir, "downloaded_data")
+os.makedirs(save_path, exist_ok=True)
+
+results = Fido.search(
+    a.Time(f"{date} 00:00:00", f"{date} 00:00:30", near=True),
+    a.Instrument("AIA"),
+    a.Wavelength(wavelength)
+)
+
+if results:
+        downloaded_files = Fido.fetch(results[0], path=save_path)
+        if downloaded_files:
+            print(f"Downloaded file: {downloaded_files[0]}")
+        else:
+            print("Download failed.")
+else:
+    print("No files found for the specified date and wavelength.")
+"""
+
+def imhmi(m171):
+    """Create a dummy HMI map for the original CHIMERA code."""
+    imhmi = Map(np.zeros_like(m171.data), m171.meta)
+    return imhmi
+
+m171 = Map("C:\\Users\\aoife\\CHIMERApy\\chimerapy\\scripts\\downloaded_data\\aia.lev1.171A_2016_10_31T00_00_10.35Z.image_lev1.fits")
+m193 = Map("C:\\Users\\aoife\\CHIMERApy\\chimerapy\\scripts\\downloaded_data\\aia.lev1.193A_2016_10_31T00_00_05.85Z.image_lev1.fits")
+m211 = Map("C:\\Users\\aoife\\CHIMERApy\\chimerapy\\scripts\\downloaded_data\\aia.lev1.211A_2016_10_31T00_00_10.62Z.image_lev1.fits")
+imhmi_map = imhmi(m171)
+
+labeled_mask, ch_mask, coronal_holes = chimera_new.chimera(m171, m193, m211)
+
+if coronal_holes:
+    print("These are the new CHIMERA results")
+else:
+    print("No coronal holes found for new CHIMERA.")
+
+circ, data, datb, datc, dattoarc, hedb, iarr, props, rs, slate, center, xgrid, ygrid = chimera_old.chimera(m171, m193, m211, imhmi_map)
+
+print("Chimera_original processing complete.")
+
+num_coronal_holes = 0
+for i in range(1, props.shape[1]):
+    if props[0, i].isdigit():
+        num_coronal_holes += 1
+print(f"Chimera_original found {num_coronal_holes} coronal holes.")
+
+for i in range(1, num_coronal_holes + 1):
+    print(f"Coronal Hole {i}:")
+    print(f"  ID: {props[0, i]}")
+    print(f"  XCEN: {props[1, i]}")
+    print(f"  YCEN: {props[2, i]}")
+    print(f"  CENTROID: {props[3, i]}")
+    print(f"  X_EB: {props[4, i]}")
+    print(f"  Y_EB: {props[5, i]}")
+    print(f"  X_WB: {props[6, i]}")
+    print(f"  Y_WB: {props[7, i]}")
+    print(f"  X_NB: {props[8, i]}")
+    print(f"  Y_NB: {props[9, i]}")
+    print(f"  X_SB: {props[10, i]}")
+    print(f"  Y_SB: {props[11, i]}")
+    print(f"  WIDTH: {props[12, i]}")
+    print(f"  WIDTH: {props[13, i]}")
+    print(f"  AREA: {props[14, i]}")
+    print(f"  AREA%: {props[15, i]}")
+    print(f"  <B>: {props[16, i]}")
+    print(f"  <B+>: {props[17, i]}")
+    print(f"  <B->: {props[18, i]}")
+    print(f"  BMAX: {props[19, i]}")
+    print(f"  BMIN: {props[20, i]}")
+    print(f"  TOT_B+: {props[21, i]}")
+    print(f"  TOT_B-: {props[22, i]}")
+    print(f"  <PHI>: {props[23, i]}")
+    print(f"  <PHI+>: {props[24, i]}")
+    print(f"  <PHI->: {props[25, i]}")