diff --git a/.gitignore b/.gitignore index 6d6f1ba..b76399a 100644 --- a/.gitignore +++ b/.gitignore @@ -261,3 +261,4 @@ package.json *.log .project +downloaded_data/ \ No newline at end of file diff --git a/chimerapy/chimera.py b/chimerapy/chimera.py index d9d4d39..b261651 100644 --- a/chimerapy/chimera.py +++ b/chimerapy/chimera.py @@ -136,7 +136,6 @@ def filter_ch(mask: NDArray, map_obj: Map, min_area: Quantity["area"] = 1e4 * u. 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) @@ -231,7 +230,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 ch_mask, labeled_mask + return ch_mask, labeled_mask, coronal_holes if __name__ == "__main__": diff --git a/chimerapy/chimera_original.py b/chimerapy/chimera_original.py index b20b4cb..216b720 100644 --- a/chimerapy/chimera_original.py +++ b/chimerapy/chimera_original.py @@ -49,7 +49,6 @@ def chimera(im171, im193, im211, imhmi): # =====Reads in data and resizes images===== ext_num = 0 if im171[0].endswith("fts.gz") else 1 x = np.arange(0, 1024) * 4 - hdu_number = 0 heda = fits.getheader(im171[0], ext_num) data = fits.getdata(im171[0], ext=ext_num) / (heda["EXPTIME"]) dn = RectBivariateSpline(x, x, data, kx=1, ky=1) @@ -62,7 +61,7 @@ def chimera(im171, im193, im211, imhmi): datc = fits.getdata(im211[0], ext=ext_num) / (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) + hedm = fits.getheader(imhmi[0], ext=0) datm = fits.getdata(imhmi[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) diff --git a/chimerapy/tests/test_compare.py b/chimerapy/tests/test_compare.py index 1281bcd..312369b 100644 --- a/chimerapy/tests/test_compare.py +++ b/chimerapy/tests/test_compare.py @@ -3,6 +3,7 @@ from astropy.utils.data import download_file +import matplotlib.pyplot as plt from chimerapy.chimera import chimera from chimerapy.chimera_original import chimera as chimera_original @@ -62,6 +63,6 @@ def test_compare(p171, p193, p211, pmag, m171, m193, m211): current = chimera(m171, m193, m211) # noqa F841 original = chimera_original(p171, p193, p211, pmag) # noqa F841 - # plt.imshow(original[6].reshape(1024, 4, 1024, 4).sum(axis=(1, 3))) - # plt.contour(current[1]) - # plt.imshow(m193.data, vmin=50, vmax=500) + plt.imshow(original[6].reshape(1024, 4, 1024, 4).sum(axis=(1, 3))) + plt.contour(current[1]) + plt.imshow(m193.data, vmin=50, vmax=500) diff --git a/scripts/comparison.py b/scripts/comparison.py new file mode 100644 index 0000000..6d2ac03 --- /dev/null +++ b/scripts/comparison.py @@ -0,0 +1,136 @@ +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 matplotlib.pyplot as plt +import chimerapy.chimera_original as chimera_old +import warnings +from astropy.io.fits.verify import VerifyWarning + +warnings.filterwarnings("ignore", category=VerifyWarning) + +""" +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.") +""" +path_171 = r"scripts\downloaded_data\AIA20161031_1156_0171.fits" +path_193 = r"scripts\downloaded_data\AIA20161031_1156_0193.fits" +path_211 = r"scripts\downloaded_data\AIA20161031_1156_0211.fits" + +m171 = Map(path_171) +m193 = Map(path_193) +m211 = Map(path_211) +imhmi = "https://solarmonitor.org/data/2016/09/22/fits/shmi/shmi_maglc_fd_20160922_094640.fts.gz" + +ch_mask, labeled_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([path_171], [path_193], [path_211], [imhmi]) + +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" : {props[16, i]}") + print(f" : {props[17, i]}") + print(f" : {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" : {props[23, i]}") + print(f" : {props[24, i]}") + print(f" : {props[25, i]}") + +old_mask = (iarr.reshape(1024, 4, 1024, 4).sum(axis=(1, 3))).astype(np.int8) +new_mask = (labeled_mask).astype(np.int8) +diff_map = old_mask - new_mask + +fig, axes = plt.subplots(1, 2, figsize=(10, 5)) + +axes[0].imshow(old_mask) +axes[0].set_title("Old Coronal Hole Mask") +axes[0].axis('off') + +axes[1].imshow(new_mask) +axes[1].set_title("New Coronal Hole Mask") +axes[1].axis('off') + +plt.tight_layout() + +fig, axes = plt.subplots(1, 2, figsize=(10, 5)) + +axes[0].imshow(m193.data, vmin=50, vmax=500) +axes[0].contour(old_mask, colors='w', linewidths=0.5) +axes[0].set_title("Old Coronal Hole Overlay") +axes[0].axis('off') + +axes[1].imshow(m193.data, vmin=50, vmax=500) +axes[1].contour(new_mask, colors='w', linewidths=0.5) +axes[1].set_title("New Coronal Hole Overlay") +axes[1].axis('off') + +plt.tight_layout() + +fig, axes = plt.subplots(1, 2, figsize=(10, 5)) + +axes[0].imshow(diff_map, cmap='bwr', vmin=-1, vmax=1) +axes[0].set_title("Difference Map (New(Blue) vs Old(Red))") +axes[0].axis('off') + +axes[1].imshow(m193.data, vmin=50, vmax=500) +axes[1].contour(old_mask, colors='r', linewidths=0.5) +axes[1].contour(new_mask, colors='w', linewidths=0.5) +axes[1].set_title("Both Coronal Hole Overlays New(White) vs Old(Red)") +axes[1].axis('off') + +plt.tight_layout() +plt.show() \ No newline at end of file