diff --git a/chimerapy/chimera.py b/chimerapy/chimera.py index 4462553..098bb9e 100644 --- a/chimerapy/chimera.py +++ b/chimerapy/chimera.py @@ -1,8 +1,4 @@ -""" - -""" - - +"""Package for Coronal Hole Identification Algorithm""" import glob import sys @@ -11,63 +7,115 @@ import mahotas import matplotlib.pyplot as plt import numpy as np +import scipy import scipy.interpolate +import sunpy import sunpy.map from astropy import wcs from astropy.io import fits from astropy.modeling.models import Gaussian2D +from astropy.visualization import astropy_mpl_style from skimage.util import img_as_ubyte +plt.style.use(astropy_mpl_style) -def chimera_legacy(): - file_path = "./" +"""loading in the images as fits files""" - 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") +file_path = "./" - circ, data, datb, datc, dattoarc, hedb, iarr, props, rs, slate, center, xgrid, ygrid = chimera( - im171, im193, im211, imhmi - ) - # =====sets ident back to max value of iarr====== - - # ident = ident - 1 - np.savetxt("ch_summary.txt", props, fmt="%s") - - plot_tricolor(data, datb, datc, xgrid, ygrid, slate) - plot_mask(slate, iarr, circ, rs, dattoarc, center, xgrid, ygrid, hedb) - - -def chimera(im171, im193, im211, imhmi): - if im171 == [] or im193 == [] or im211 == [] or imhmi == []: - print("Not all required files present") - sys.exit() - # =====Reads in data and resizes images===== - x = np.arange(0, 1024) * 4 - hdu_number = 0 - heda = fits.getheader(im171[0], hdu_number) - data = fits.getdata(im171[0], ext=0) / (heda["EXPTIME"]) - dn = scipy.interpolate.interp2d(x, x, data) - data = dn(np.arange(0, 4096), np.arange(0, 4096)) - hedb = fits.getheader(im193[0], hdu_number) - datb = fits.getdata(im193[0], ext=0) / (hedb["EXPTIME"]) - dn = scipy.interpolate.interp2d(x, x, datb) - datb = dn(np.arange(0, 4096), np.arange(0, 4096)) - hedc = fits.getheader(im211[0], hdu_number) - datc = fits.getdata(im211[0], ext=0) / (hedc["EXPTIME"]) - dn = scipy.interpolate.interp2d(x, x, datc) - datc = dn(np.arange(0, 4096), np.arange(0, 4096)) - hedm = fits.getheader(imhmi[0], hdu_number) - 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) - if hedm["crota1"] > 90: - datm = np.rot90(np.rot90(datm)) - # =====Specifies solar radius and calculates conversion value of pixel to arcsec===== - s = np.shape(data) - rs = heda["rsun"] +im171 = glob.glob(file_path + "*171*.fts") +im193 = glob.glob(file_path + "*193*.fts") +im211 = glob.glob(file_path + "*211*.fts") +imhmi = glob.glob(file_path + "*hmi*.fts") + +"""ensure that all images are present""" + +if im171 == [] or im193 == [] or im211 == [] or imhmi == []: + print("Not all required files present") + sys.exit() + + +def rescale_aia(image: np.array, orig_res: int, desired_res: int): + """ + Rescale the input aia image dimensions. + + Parameters + ---------- + image: 'np.array' + orig_res: 'int' + desired_res: 'int + + Returns + ------- + 'np.array' + """ + + if desired_res > orig_res: + scaled_array = np.linspace(start=0, stop=desired_res, num=orig_res) + dn = scipy.interpolate.RectBivariateSpline( + scaled_array, scaled_array, fits.getdata(image[0], 0) / (fits.getheader(image[0], 0)["EXPTIME"]) + ) + return dn(np.arange(0, desired_res), np.arange(0, desired_res)) + elif desired_res < orig_res: + scaled_array = np.linspace(start=0, stop=orig_res, num=desired_res) + dn = scipy.interpolate.RectBivariateSpline(scaled_array, scaled_array, fits.getdata(image[0], 0)) + return dn(np.arange(0, desired_res), np.arange(0, desired_res)) + + +def rescale_hmi(image: np.array, orig_res: int, desired_res: int): + """ + Rescale the input hmi image dimensions. + + Parameters + ---------- + image: 'np.array' + orig_res: 'int' + desired_res: 'int + + Returns + ------- + 'np.array' + """ + if desired_res > orig_res: + scaled_array = np.linspace(start=0, stop=desired_res, num=orig_res) + dn = scipy.interpolate.RectBivariateSpline(scaled_array, scaled_array, fits.getdata(image[0], ext=0)) + return dn(np.arange(0, desired_res), np.arange(0, desired_res)) + elif desired_res < orig_res: + scaled_array = np.linspace(start=0, stop=orig_res, num=desired_res) + dn = scipy.interpolate.RectBivariateSpline(scaled_array, scaled_array, fits.getdata(image[0], ext=0)) + return dn(np.arange(0, desired_res), np.arange(0, desired_res)) + + +"""defining data arrays which are used in later steps""" + +data = rescale_aia(im171, 1024, 4096) +datb = rescale_aia(im193, 1024, 4096) +datc = rescale_aia(im211, 1024, 4096) +datm = rescale_hmi(imhmi, 1024, 4096) + + +def filter(aiaa: np.array, aiab: np.array, aiac: np.array, aiam: np.array): + """ + Defines headers and filters aia arrays to meet header requirements + + Parameters + ---------- + aiaa: 'np.array' + aiab: 'np.array' + aiac: 'np.array' + aiam: 'np.array' + + Returns + ------- + 'np.array' + + """ + global heda, hedb, hedc, hedm, datm + heda = fits.getheader(aiaa[0], 0) + hedb = fits.getheader(aiab[0], 0) + hedc = fits.getheader(aiac[0], 0) + hedm = fits.getheader(aiam[0], 0) if hedb["ctype1"] != "solar_x ": hedb["ctype1"] = "solar_x " hedb["ctype2"] = "solar_y " @@ -90,136 +138,355 @@ def chimera(im171, im193, im211, imhmi): hedc["crpix1"] * 4.0, hedc["crpix2"] * 4.0, ) + if hedm["crota1"] > 90: + datm = np.rot90(np.rot90(datm)) + + +filter(im171, im193, im211, imhmi) + + +def remove_neg(aiaa: np.array, aiab: np.array, aiac: np.array): + """ + Removes negative values from arrays + + Parameters + ---------- + aiaa: 'np.array' + aiab: 'np.array' + aiac: 'np.array' + + Returns + ------- + 'np.array' + + """ + global data, datb, datc + data[np.where(data <= 0)] = 0 + datb[np.where(datb <= 0)] = 0 + datc[np.where(datc <= 0)] = 0 + + +remove_neg(im171, im193, im211) + +"""defines the shape of the arrays as "s" and "rs" as the solar radius""" +s = np.shape(data) +rs = heda["rsun"] + + +def pix_arc(aia: np.array): + global dattoarc dattoarc = heda["cdelt1"] + global conver + conver = ((s[0]) / 2) * dattoarc / hedm["cdelt1"] - (s[1] / 2) + global convermul convermul = dattoarc / hedm["cdelt1"] - # =====Alternative coordinate systems===== - hdul = fits.open(im171[0]) - hdul[0].header["CUNIT1"] = "arcsec" - hdul[0].header["CUNIT2"] = "arcsec" - aia = sunpy.map.Map(hdul[0].data, hdul[0].header) - adj = 4096.0 / aia.dimensions[0].value + + +pix_arc(im171) + + +def to_helio(image: np.array): + """ + Converts arrays to the Heliographic Stonyhurst coordinate system + + Parameters + ---------- + image: 'np.array' + + Returns + ------- + 'np.array' + + """ + + aia = sunpy.map.Map(image) + adj = 4096 / aia.dimensions[0].value x, y = (np.meshgrid(*[np.arange(adj * v.value) for v in aia.dimensions]) * u.pixel) / adj + print(x, y) + global hpc hpc = aia.pixel_to_world(x, y) + global hg hg = hpc.transform_to(sunpy.coordinates.frames.HeliographicStonyhurst) + global csys csys = wcs.WCS(hedb) - # =======setting up arrays to be used============ - ident = 1 - iarr = np.zeros((s[0], s[1]), dtype=np.byte) - offarr, slate = np.array(iarr), np.array(iarr) - bmcool = np.zeros((s[0], s[1]), dtype=np.float32) - cand, bmmix, bmhot = np.array(bmcool), np.array(bmcool), np.array(bmcool) - circ = np.zeros((s[0], s[1]), dtype=int) - # =======creation of a 2d gaussian for magnetic cut offs=========== - r = (s[1] / 2.0) - 450 - xgrid, ygrid = np.meshgrid(np.arange(s[0]), np.arange(s[1])) - center = [int(s[1] / 2.0), int(s[1] / 2.0)] - w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 > r**2) - y, x = np.mgrid[0:4096, 0:4096] - garr = Gaussian2D(1, s[0] / 2, s[1] / 2, 2000 / 2.3548, 2000 / 2.3548)(x, y) - garr[w] = 1.0 - # ======creation of array for CH properties========== - props = np.zeros((26, 30), dtype="", - "", - "", - "BMAX", - "BMIN", - "TOT_B+", - "TOT_B-", - "", - "", - "", - ) - props[:, 1] = ( - "num", - '"', - '"', - "H deg", - '"', - '"', - '"', - '"', - '"', - '"', - '"', - '"', - "H deg", - "deg", - "Mm^2", - "%", - "G", - "G", - "G", - "G", - "G", - "G", - "G", - "Mx", - "Mx", - "Mx", - ) - # =====removes negative data values===== - data[np.where(data <= 0)] = 0 - datb[np.where(datb <= 0)] = 0 - datc[np.where(datc <= 0)] = 0 - # ============make a multi-wavelength image for contours================== - with np.errstate(divide="ignore"): - t0 = np.log10(datc) - t1 = np.log10(datb) - t2 = np.log10(data) - t0[np.where(t0 < 0.8)] = 0.8 - t0[np.where(t0 > 2.7)] = 2.7 - t1[np.where(t1 < 1.4)] = 1.4 - t1[np.where(t1 > 3.0)] = 3.0 - t2[np.where(t2 < 1.2)] = 1.2 - t2[np.where(t2 > 3.9)] = 3.9 - t0 = np.array(((t0 - 0.8) / (2.7 - 0.8)) * 255, dtype=np.float32) - t1 = np.array(((t1 - 1.4) / (3.0 - 1.4)) * 255, dtype=np.float32) - t2 = np.array(((t2 - 1.2) / (3.9 - 1.2)) * 255, dtype=np.float32) - # ====create 3 segmented bitmasks===== + + +to_helio(im171) + +"""Setting up arrays to be used in later processing""" +ident = 1 +iarr = np.zeros((s[0], s[1]), dtype=np.byte) +bmcool = np.zeros((s[0], s[1]), dtype=np.float32) +offarr, slate = np.array(iarr), np.array(iarr) +cand, bmmix, bmhot = np.array(bmcool), np.array(bmcool), np.array(bmcool) +circ = np.zeros((s[0], s[1]), dtype=int) + +"""creation of a 2d gaussian for magnetic cut offs""" +r = (s[1] / 2.0) - 450 +xgrid, ygrid = np.meshgrid(np.arange(s[0]), np.arange(s[1])) +center = [int(s[1] / 2.0), int(s[1] / 2.0)] +w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 > r**2) +y, x = np.mgrid[0:4096, 0:4096] +garr = Gaussian2D(1, s[0] / 2, s[1] / 2, 2000 / 2.3548, 2000 / 2.3548)(x, y) +garr[w] = 1.0 + +"""creates sub-arrays of props to isolate column of index 0 and column of index 1""" +props = np.zeros((26, 30), dtype="", + "", + "", + "BMAX", + "BMIN", + "TOT_B+", + "TOT_B-", + "", + "", + "", +) +props[:, 1] = ( + "num", + '"', + '"', + "H°", + '"', + '"', + '"', + '"', + '"', + '"', + '"', + '"', + "H°", + "°", + "Mm^2", + "%", + "G", + "G", + "G", + "G", + "G", + "G", + "G", + "Mx", + "Mx", + "Mx", +) + +"""define threshold values in log space""" + +with np.errstate(divide="ignore"): + t0 = np.log10(datc) + t1 = np.log10(datb) + t2 = np.log10(data) + + +class Bounds: + """Mixin to change and define array boundaries and slopes""" + + def __init__(self, upper, lower, slope): + self.upper = upper + self.lower = lower + self.slope = slope + + def new_u(self, new_upper): + self.upper = new_upper + + def new_l(self, new_lower): + self.lower = new_lower + + def new_s(self, new_slope): + self.slope = new_slope + + +t0b = Bounds(0.8, 2.7, 255) +t1b = Bounds(1.4, 3.0, 255) +t2b = Bounds(1.2, 3.9, 255) + + +def threshold(tval: np.array): + """ + Threshold arrays based on desired boundaries + + Parameters + ---------- + tval: 'np.array' + + Returns + ------- + 'np.array' + + """ + global t0, t1, t2 + if tval.all() == t0.all(): + t0[np.where(t0 < t0b.upper)] = t0b.upper + t0[np.where(t0 > t0b.lower)] = t0b.lower + if tval.all() == t1.all(): + t1[np.where(t1 < t1b.upper)] = t1b.upper + t1[np.where(t1 > t1b.lower)] = t2b.lower + if tval.all() == t2.all(): + t2[np.where(t2 < t2b.upper)] = t2b.upper + t2[np.where(t2 > t2b.lower)] = t2b.lower + else: + print("Must input valid logarithmic arrays") + + +threshold(t0) +threshold(t1) +threshold(t2) + + +def set_contour(tval: np.array): + """Sets contour values for bounded arrays + + Parameters + ---------- + tval: 'np.array' + + Returns + ------- + 'np.array' + + """ + global t0, t1, t2 + if tval.all() == t0.all(): + t0 = np.array(((t0 - t0b.upper) / (t0b.lower - t0b.upper)) * t0b.slope, dtype=np.float32) + elif tval.all() == t1.all(): + t1 = np.array(((t1 - t1b.upper) / (t1b.lower - t1b.upper)) * t1b.slope, dtype=np.float32) + elif tval.all() == t2.all(): + t2 = np.array(((t2 - t2b.upper) / (t2b.lower - t2b.upper)) * t2b.slope, dtype=np.float32) + + +set_contour(t0) +set_contour(t1) +set_contour(t2) + + +def create_mask(): + """ + Creates 3 segmented bitmasks + + Returns + ------- + 'np.array' + + """ + + global t0, t1, t2, bmmix, bmhot, bmcool with np.errstate(divide="ignore", invalid="ignore"): bmmix[np.where(t2 / t0 >= ((np.mean(data) * 0.6357) / (np.mean(datc))))] = 1 bmhot[np.where(t0 + t1 < (0.7 * (np.mean(datb) + np.mean(datc))))] = 1 bmcool[np.where(t2 / t1 >= ((np.mean(data) * 1.5102) / (np.mean(datb))))] = 1 - # ====logical conjunction of 3 segmentations======= + + +create_mask() + + +def conjunction(): + """ + Creates a conjunction of 3 segmentations + + Returns + ------- + 'np.array' + + """ + global bmhot, bmcool, bmmix, cand cand = bmcool * bmmix * bmhot - # ====plot tricolour image with lon/lat conotours======= - # ======removes off detector mis-identifications========== + + +conjunction() + + +def misid(): + """ + Removes off-detector mis-identification + + Returns + ------- + 'np.array' + + """ + global s, r, w, circ, cand r = (s[1] / 2.0) - 100 w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 <= r**2) circ[w] = 1.0 cand = cand * circ - # =======Separates on-disk and off-limb CHs=============== + + +misid() + + +def on_off(): + """ + Seperates on-disk and off-limb coronal holes + + Returns + ------- + 'np.array' + + """ + global circ, cand circ[:] = 0 r = (rs / dattoarc) - 10 - w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 <= r**2) - circ[w] = 1.0 + inside = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 <= r**2) + circ[inside] = 1.0 r = (rs / dattoarc) + 40 - w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 >= r**2) - circ[w] = 1.0 + outside = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 >= r**2) + circ[outside] = 1.0 cand = cand * circ - # ====open file for property storage===== - # =====contours the identified datapoints======= + + +on_off() + + +def contours(): + """ + Contours the identified datapoints + + Returns + ------- + cand: 'np.array' + cont: 'tuple' + heir: 'np.array' + + """ + global cand, cont, heir cand = np.array(cand, dtype=np.uint8) cont, heir = cv2.findContours(cand, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) - # ======sorts contours by size============ + + +contours() + + +def sort(): + """ + Sorts the contours by size + + Returns + ------- + reord: 'list' + tmp: 'list' + cont: 'list' + + """ + global sizes, reord, tmp, cont sizes = [] for i in range(len(cont)): sizes = np.append(sizes, len(cont[i])) @@ -228,257 +495,449 @@ def chimera(im171, im193, im211, imhmi): for i in range(len(cont)): tmp[i] = cont[reord[i]] cont = list(tmp) - # =====cycles through contours========= - for i in range(len(cont)): - x = np.append(x, len(cont[i])) - # =====only takes values of minimum surface length and calculates area====== - if len(cont[i]) <= 100: - continue - area = 0.5 * np.abs( - np.dot(cont[i][:, 0, 0], np.roll(cont[i][:, 0, 1], 1)) - - np.dot(cont[i][:, 0, 1], np.roll(cont[i][:, 0, 0], 1)) - ) - arcar = area * (dattoarc**2) - if arcar > 1000: - # =====finds centroid======= - - # chpts = len(cont[i]) - cent = [np.mean(cont[i][:, 0, 0]), np.mean(cont[i][:, 0, 1])] - - # ===remove quiet sun regions encompassed by coronal holes====== - - if ( - cand[ - np.max(cont[i][:, 0, 0]) + 1, - cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], - ] - > 0 - ) and ( - iarr[ - np.max(cont[i][:, 0, 0]) + 1, - cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], - ] - > 0 +sort() + + +# =====cycles through contours========= + + +def extent(i, ypos, xpos, hg, cont): + """ + Finds coronal hole extent in latitude and longitude + + Parameters + ---------- + i: 'int' + ypos: 'astropy.units.quantity.Quantity' + xpos: 'astropy.units.quantity.Quantity' + hg: 'astropy.coordinates.sky_coordinate.SkyCoord' + cont: 'list' + + Returns + ------- + maxxlon: 'astropy.coordinates.angles.core.Longitude' + minxlon: 'astropy.coordinates.angles.core.Longitude' + centlat: 'astropy.coordinates.angles.core.Latitude' + centlon: 'astropy.coordinates.angles.core.Longitude' + + """ + global maxxlat, maxxlon, maxylat, maxylon, minxlon, minylat, minylon, minxlat + maxxlat = hg.lat[ + cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], + np.max(cont[i][:, 0, 0]), + ] + maxxlon = hg.lon[ + cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], + np.max(cont[i][:, 0, 0]), + ] + maxylat = hg.lat[ + np.max(cont[i][:, 0, 1]), + cont[i][np.where(cont[i][:, 0, 1] == np.max(cont[i][:, 0, 1]))[0][0], 0, 0], + ] + maxylon = hg.lon[ + np.max(cont[i][:, 0, 1]), + cont[i][np.where(cont[i][:, 0, 1] == np.max(cont[i][:, 0, 1]))[0][0], 0, 0], + ] + minxlat = hg.lat[ + cont[i][np.where(cont[i][:, 0, 0] == np.min(cont[i][:, 0, 0]))[0][0], 0, 1], + np.min(cont[i][:, 0, 0]), + ] + minxlon = hg.lon[ + cont[i][np.where(cont[i][:, 0, 0] == np.min(cont[i][:, 0, 0]))[0][0], 0, 1], + np.min(cont[i][:, 0, 0]), + ] + minylat = hg.lat[ + np.min(cont[i][:, 0, 1]), + cont[i][np.where(cont[i][:, 0, 1] == np.min(cont[i][:, 0, 1]))[0][0], 0, 0], + ] + minylon = hg.lon[ + np.min(cont[i][:, 0, 1]), + cont[i][np.where(cont[i][:, 0, 1] == np.min(cont[i][:, 0, 1]))[0][0], 0, 0], + ] + + # =====CH centroid in lat/lon======= + + centlat = hg.lat[int(ypos), int(xpos)] + centlon = hg.lon[int(ypos), int(xpos)] + return maxxlon, minxlon, centlat, centlon + + +def coords(i, csys, cont): + """ + Finds coordinates of CH boundaries + + Parameters + ---------- + i: 'int' + csys: 'astropy.wcs.wcs.WCS' + cont: 'list' + + Returns + ------- + Ywb: 'np.array' + Xwb: 'np.array' + Yeb: 'np.array' + Xeb: 'np.array' + Ynb: 'np.array' + Xnb: 'np.array' + Ysb: 'np.array' + Xsb: 'np.array' + """ + global Ywb, Xwb, Yeb, Xeb, Ynb, Xnb, Ysb, Xsb + Ywb, Xwb = csys.all_pix2world( + cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], + np.max(cont[i][:, 0, 0]), + 0, + ) + Yeb, Xeb = csys.all_pix2world( + cont[i][np.where(cont[i][:, 0, 0] == np.min(cont[i][:, 0, 0]))[0][0], 0, 1], + np.min(cont[i][:, 0, 0]), + 0, + ) + Ynb, Xnb = csys.all_pix2world( + np.max(cont[i][:, 0, 1]), + cont[i][np.where(cont[i][:, 0, 1] == np.max(cont[i][:, 0, 1]))[0][0], 0, 0], + 0, + ) + Ysb, Xsb = csys.all_pix2world( + np.min(cont[i][:, 0, 1]), + cont[i][np.where(cont[i][:, 0, 1] == np.min(cont[i][:, 0, 1]))[0][0], 0, 0], + 0, + ) + + return Ywb, Xwb, Yeb, Xeb, Ynb, Xnb, Ysb, Xsb + + +def ins_prop( + datm, + rs, + ident, + props, + arcar, + arccent, + pos, + npix, + trummar, + centlat, + centlon, + mB, + mBpos, + mBneg, + Ywb, + Xwb, + Yeb, + Xeb, + Ynb, + Xnb, + Ysb, + Xsb, + width, + eastl, + westl, +): + """ + Insertion of CH properties into property array + + Parameters + ---------- + datm: 'np.array' + rs: 'float' + ident: 'int' + props: 'np.array' + arcar: 'np.float64' + arccent: 'list' + pos: 'np.array' + npix: 'list' + trummar: 'np.float64' + centlat: 'str' + centlon: 'str' + mB: 'np.float64' + mBpos: 'np.float64' + mBneg: 'np.float64' + + Returns + ------- + props[0, ident + 1]: 'str' + props[1, ident + 1]: 'str' + props[2, ident + 1]: 'str' + props[3, ident + 1]: 'str' + props[4, ident + 1]: 'str' + props[5, ident + 1]: 'str' + props[6, ident + 1]: 'str' + props[7, ident + 1]: 'str' + props[8, ident + 1]: 'str' + props[9, ident + 1]: 'str' + props[10, ident + 1]: 'str' + props[11, ident + 1]: 'str' + props[12, ident + 1]: 'str' + props[13, ident + 1]: 'str' + props[14, ident + 1]: 'str' + props[15, ident + 1]: 'str' + props[16, ident + 1]: 'str' + props[17, ident + 1]: 'str' + props[18, ident + 1]: 'str' + props[19, ident + 1]: 'str' + props[20, ident + 1]: 'str' + tbpos: 'np.float64' + props[21, ident + 1]: 'str' + tbneg: 'np.float64' + props[22, ident + 1]: 'str' + props[23, ident + 1]: 'str' + props[24, ident + 1]: 'str' + props[25, ident + 1]: 'str' + + """ + props[0, ident + 1] = str(ident) + props[1, ident + 1] = str(np.round(arccent[0])) + props[2, ident + 1] = str(np.round(arccent[1])) + props[3, ident + 1] = str(centlon + centlat) + props[4, ident + 1] = str(np.round(Xeb)) + props[5, ident + 1] = str(np.round(Yeb)) + props[6, ident + 1] = str(np.round(Xwb)) + props[7, ident + 1] = str(np.round(Ywb)) + props[8, ident + 1] = str(np.round(Xnb)) + props[9, ident + 1] = str(np.round(Ynb)) + props[10, ident + 1] = str(np.round(Xsb)) + props[11, ident + 1] = str(np.round(Ysb)) + props[12, ident + 1] = str(eastl + "-" + westl) + props[13, ident + 1] = str(width) + props[14, ident + 1] = f"{trummar/1e+12:.1e}" + props[15, ident + 1] = str(np.round((arcar * 100 / (np.pi * (rs**2))), 1)) + props[16, ident + 1] = str(np.round(mB, 1)) + props[17, ident + 1] = str(np.round(mBpos, 1)) + props[18, ident + 1] = str(np.round(mBneg, 1)) + props[19, ident + 1] = str(np.round(np.max(npix[1]), 1)) + props[20, ident + 1] = str(np.round(np.min(npix[1]), 1)) + tbpos = np.sum(datm[pos[:, 0], pos[:, 1]][np.where(datm[pos[:, 0], pos[:, 1]] > 0)]) + props[21, ident + 1] = f"{tbpos:.1e}" + tbneg = np.sum(datm[pos[:, 0], pos[:, 1]][np.where(datm[pos[:, 0], pos[:, 1]] < 0)]) + props[22, ident + 1] = f"{tbneg:.1e}" + props[23, ident + 1] = f"{mB*trummar*1e+16:.1e}" + props[24, ident + 1] = f"{mBpos*trummar*1e+16:.1e}" + props[25, ident + 1] = f"{mBneg*trummar*1e+16:.1e}" + + +"""Cycles through contours""" + +for i in range(len(cont)): + x = np.append(x, len(cont[i])) + + """only takes values of minimum surface length and calculates area""" + + if len(cont[i]) <= 100: + continue + area = 0.5 * np.abs( + np.dot(cont[i][:, 0, 0], np.roll(cont[i][:, 0, 1], 1)) + - np.dot(cont[i][:, 0, 1], np.roll(cont[i][:, 0, 0], 1)) + ) + arcar = area * (dattoarc**2) + if arcar > 1000: + """finds centroid""" + + chpts = len(cont[i]) + cent = [np.mean(cont[i][:, 0, 0]), np.mean(cont[i][:, 0, 1])] + + """remove quiet sun regions encompassed by coronal holes""" + if ( + cand[ + np.max(cont[i][:, 0, 0]) + 1, + cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], + ] + > 0 + ) and ( + iarr[ + np.max(cont[i][:, 0, 0]) + 1, + cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], + ] + > 0 + ): + mahotas.polygon.fill_polygon(np.array(list(zip(cont[i][:, 0, 1], cont[i][:, 0, 0]))), slate) + print(slate) + iarr[np.where(slate == 1)] = 0 + slate[:] = 0 + + else: + """Create a simple centre point if coronal hole regions is not quiet""" + + arccent = csys.all_pix2world(cent[0], cent[1], 0) + + """classifies off limb CH regions""" + + if (((arccent[0] ** 2) + (arccent[1] ** 2)) > (rs**2)) or ( + np.sum(np.array(csys.all_pix2world(cont[i][0, 0, 0], cont[i][0, 0, 1], 0)) ** 2) > (rs**2) ): + mahotas.polygon.fill_polygon(np.array(list(zip(cont[i][:, 0, 1], cont[i][:, 0, 0]))), offarr) + else: + """classifies on disk coronal holes""" + mahotas.polygon.fill_polygon(np.array(list(zip(cont[i][:, 0, 1], cont[i][:, 0, 0]))), slate) - iarr[np.where(slate == 1)] = 0 + poslin = np.where(slate == 1) slate[:] = 0 + print(poslin) + + """create an array for magnetic polarity""" + + pos = np.zeros((len(poslin[0]), 2), dtype=np.uint) + pos[:, 0] = np.array((poslin[0] - (s[0] / 2)) * convermul + (s[1] / 2), dtype=np.uint) + pos[:, 1] = np.array((poslin[1] - (s[0] / 2)) * convermul + (s[1] / 2), dtype=np.uint) + npix = list( + np.histogram( + datm[pos[:, 0], pos[:, 1]], + bins=np.arange( + np.round(np.min(datm[pos[:, 0], pos[:, 1]])) - 0.5, + np.round(np.max(datm[pos[:, 0], pos[:, 1]])) + 0.6, + 1, + ), + ) + ) + npix[0][np.where(npix[0] == 0)] = 1 + npix[1] = npix[1][:-1] + 0.5 - else: - # ====create a simple centre point====== - - arccent = csys.all_pix2world(cent[0], cent[1], 0) + wh1 = np.where(npix[1] > 0) + wh2 = np.where(npix[1] < 0) - # ====classifies off limb CH regions======== + """Filters magnetic cutoff values by area""" - if (((arccent[0] ** 2) + (arccent[1] ** 2)) > (rs**2)) or ( - np.sum(np.array(csys.all_pix2world(cont[i][0, 0, 0], cont[i][0, 0, 1], 0)) ** 2) > (rs**2) + if ( + np.absolute((np.sum(npix[0][wh1]) - np.sum(npix[0][wh2])) / np.sqrt(np.sum(npix[0]))) + <= 10 + and arcar < 9000 ): - mahotas.polygon.fill_polygon( - np.array(list(zip(cont[i][:, 0, 1], cont[i][:, 0, 0]))), offarr - ) - else: - # =====classifies on disk coronal holes======= + continue + if ( + np.absolute(np.mean(datm[pos[:, 0], pos[:, 1]])) < garr[int(cent[0]), int(cent[1])] + and arcar < 40000 + ): + continue + iarr[poslin] = ident - mahotas.polygon.fill_polygon( - np.array(list(zip(cont[i][:, 0, 1], cont[i][:, 0, 0]))), slate - ) - poslin = np.where(slate == 1) - slate[:] = 0 - - # ====create an array for magnetic polarity======== - - pos = np.zeros((len(poslin[0]), 2), dtype=np.uint) - pos[:, 0] = np.array((poslin[0] - (s[0] / 2)) * convermul + (s[1] / 2), dtype=np.uint) - pos[:, 1] = np.array((poslin[1] - (s[0] / 2)) * convermul + (s[1] / 2), dtype=np.uint) - npix = list( - np.histogram( - datm[pos[:, 0], pos[:, 1]], - bins=np.arange( - np.round(np.min(datm[pos[:, 0], pos[:, 1]])) - 0.5, - np.round(np.max(datm[pos[:, 0], pos[:, 1]])) + 0.6, - 1, - ), - ) - ) - npix[0][np.where(npix[0] == 0)] = 1 - npix[1] = npix[1][:-1] + 0.5 - - wh1 = np.where(npix[1] > 0) - wh2 = np.where(npix[1] < 0) - - # =====magnetic cut offs dependent on area========= - - if ( - np.absolute((np.sum(npix[0][wh1]) - np.sum(npix[0][wh2])) / np.sqrt(np.sum(npix[0]))) - <= 10 - and arcar < 9000 - ): - continue - if ( - np.absolute(np.mean(datm[pos[:, 0], pos[:, 1]])) < garr[int(cent[0]), int(cent[1])] - and arcar < 40000 - ): - continue - iarr[poslin] = ident - - # ====create an accurate center point======= - - ypos = np.sum((poslin[0]) * np.absolute(hg.lat[poslin])) / np.sum( - np.absolute(hg.lat[poslin]) - ) - xpos = np.sum((poslin[1]) * np.absolute(hg.lon[poslin])) / np.sum( - np.absolute(hg.lon[poslin]) - ) + """create an accurate center point""" - arccent = csys.all_pix2world(xpos, ypos, 0) - - # ======calculate average angle coronal hole is subjected to====== - - dist = np.sqrt((arccent[0] ** 2) + (arccent[1] ** 2)) - ang = np.arcsin(dist / rs) - - # =====calculate area of CH with minimal projection effects====== - - trupixar = abs(area / np.cos(ang)) - truarcar = trupixar * (dattoarc**2) - trummar = truarcar * ((6.96e08 / rs) ** 2) - - # ====find CH extent in latitude and longitude======== - - # maxxlat = hg.lat[ - # cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], - # np.max(cont[i][:, 0, 0]), - # ] - maxxlon = hg.lon[ - cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], - np.max(cont[i][:, 0, 0]), - ] - # maxylat = hg.lat[ - # np.max(cont[i][:, 0, 1]), - # cont[i][np.where(cont[i][:, 0, 1] == np.max(cont[i][:, 0, 1]))[0][0], 0, 0], - # ] - # maxylon = hg.lon[ - # np.max(cont[i][:, 0, 1]), - # cont[i][np.where(cont[i][:, 0, 1] == np.max(cont[i][:, 0, 1]))[0][0], 0, 0], - # ] - # minxlat = hg.lat[ - # cont[i][np.where(cont[i][:, 0, 0] == np.min(cont[i][:, 0, 0]))[0][0], 0, 1], - # np.min(cont[i][:, 0, 0]), - # ] - minxlon = hg.lon[ - cont[i][np.where(cont[i][:, 0, 0] == np.min(cont[i][:, 0, 0]))[0][0], 0, 1], - np.min(cont[i][:, 0, 0]), - ] - # minylat = hg.lat[ - # np.min(cont[i][:, 0, 1]), - # cont[i][np.where(cont[i][:, 0, 1] == np.min(cont[i][:, 0, 1]))[0][0], 0, 0], - # ] - # minylon = hg.lon[ - # np.min(cont[i][:, 0, 1]), - # cont[i][np.where(cont[i][:, 0, 1] == np.min(cont[i][:, 0, 1]))[0][0], 0, 0], - # ] - - # =====CH centroid in lat/lon======= - - centlat = hg.lat[int(ypos), int(xpos)] - centlon = hg.lon[int(ypos), int(xpos)] - - # ====calculate the mean magnetic field===== - - mB = np.mean(datm[pos[:, 0], pos[:, 1]]) - mBpos = np.sum(npix[0][wh1] * npix[1][wh1]) / np.sum(npix[0][wh1]) - mBneg = np.sum(npix[0][wh2] * npix[1][wh2]) / np.sum(npix[0][wh2]) - - # =====finds coordinates of CH boundaries======= - - Ywb, Xwb = csys.all_pix2world( - cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], - np.max(cont[i][:, 0, 0]), - 0, - ) - Yeb, Xeb = csys.all_pix2world( - cont[i][np.where(cont[i][:, 0, 0] == np.min(cont[i][:, 0, 0]))[0][0], 0, 1], - np.min(cont[i][:, 0, 0]), - 0, - ) - Ynb, Xnb = csys.all_pix2world( - np.max(cont[i][:, 0, 1]), - cont[i][np.where(cont[i][:, 0, 1] == np.max(cont[i][:, 0, 1]))[0][0], 0, 0], - 0, - ) - Ysb, Xsb = csys.all_pix2world( - np.min(cont[i][:, 0, 1]), - cont[i][np.where(cont[i][:, 0, 1] == np.min(cont[i][:, 0, 1]))[0][0], 0, 0], - 0, - ) + ypos = np.sum((poslin[0]) * np.absolute(hg.lat[poslin])) / np.sum(np.absolute(hg.lat[poslin])) + xpos = np.sum((poslin[1]) * np.absolute(hg.lon[poslin])) / np.sum(np.absolute(hg.lon[poslin])) + + arccent = csys.all_pix2world(xpos, ypos, 0) + + """calculate average angle coronal hole is subjected to""" + + dist = np.sqrt((arccent[0] ** 2) + (arccent[1] ** 2)) + ang = np.arcsin(dist / rs) + + """calculate area of CH with minimal projection effects""" + + trupixar = abs(area / np.cos(ang)) + truarcar = trupixar * (dattoarc**2) + trummar = truarcar * ((6.96e08 / rs) ** 2) + + """find CH extent in lattitude and longitude""" + + maxxlon, minxlon, centlat, centlon = extent(i, ypos, xpos, hg, cont) + + """caluclate the mean magnetic field""" + + mB = np.mean(datm[pos[:, 0], pos[:, 1]]) + mBpos = np.sum(npix[0][wh1] * npix[1][wh1]) / np.sum(npix[0][wh1]) + mBneg = np.sum(npix[0][wh2] * npix[1][wh2]) / np.sum(npix[0][wh2]) + + """finds coordinates of CH boundaries""" - width = round(maxxlon.value) - round(minxlon.value) - - if minxlon.value >= 0.0: - eastl = "W" + str(int(np.round(minxlon.value))) - else: - eastl = "E" + str(np.absolute(int(np.round(minxlon.value)))) - if maxxlon.value >= 0.0: - westl = "W" + str(int(np.round(maxxlon.value))) - else: - westl = "E" + str(np.absolute(int(np.round(maxxlon.value)))) - - if centlat >= 0.0: - centlat = "N" + str(int(np.round(centlat.value))) - else: - centlat = "S" + str(np.absolute(int(np.round(centlat.value)))) - if centlon >= 0.0: - centlon = "W" + str(int(np.round(centlon.value))) - else: - centlon = "E" + str(np.absolute(int(np.round(centlon.value)))) - - # ====insertions of CH properties into property array===== - - props[0, ident + 1] = str(ident) - props[1, ident + 1] = str(np.round(arccent[0])) - props[2, ident + 1] = str(np.round(arccent[1])) - props[3, ident + 1] = str(centlon + centlat) - props[4, ident + 1] = str(np.round(Xeb)) - props[5, ident + 1] = str(np.round(Yeb)) - props[6, ident + 1] = str(np.round(Xwb)) - props[7, ident + 1] = str(np.round(Ywb)) - props[8, ident + 1] = str(np.round(Xnb)) - props[9, ident + 1] = str(np.round(Ynb)) - props[10, ident + 1] = str(np.round(Xsb)) - props[11, ident + 1] = str(np.round(Ysb)) - props[12, ident + 1] = str(eastl + "-" + westl) - props[13, ident + 1] = str(width) - props[14, ident + 1] = f"{trummar / 1e+12:.1e}" - props[15, ident + 1] = str(np.round((arcar * 100 / (np.pi * (rs**2))), 1)) - props[16, ident + 1] = str(np.round(mB, 1)) - props[17, ident + 1] = str(np.round(mBpos, 1)) - props[18, ident + 1] = str(np.round(mBneg, 1)) - props[19, ident + 1] = str(np.round(np.max(npix[1]), 1)) - props[20, ident + 1] = str(np.round(np.min(npix[1]), 1)) - tbpos = np.sum(datm[pos[:, 0], pos[:, 1]][np.where(datm[pos[:, 0], pos[:, 1]] > 0)]) - props[21, ident + 1] = f"{tbpos:.1e}" - tbneg = np.sum(datm[pos[:, 0], pos[:, 1]][np.where(datm[pos[:, 0], pos[:, 1]] < 0)]) - props[22, ident + 1] = f"{tbneg:.1e}" - props[23, ident + 1] = f"{mB * trummar * 1e+16:.1e}" - props[24, ident + 1] = f"{mBpos * trummar * 1e+16:.1e}" - props[25, ident + 1] = f"{mBneg * trummar * 1e+16:.1e}" - - # =====sets up code for next possible coronal hole===== - - ident = ident + 1 - return circ, data, datb, datc, dattoarc, hedb, iarr, props, rs, slate, center, xgrid, ygrid + Ywb, Xwb, Yeb, Xeb, Ynb, Xnb, Ysb, Xsb = coords(i, csys, cont) + + width = round(maxxlon.value) - round(minxlon.value) + + if minxlon.value >= 0.0: + eastl = "W" + str(int(np.round(minxlon.value))) + else: + eastl = "E" + str(np.absolute(int(np.round(minxlon.value)))) + if maxxlon.value >= 0.0: + westl = "W" + str(int(np.round(maxxlon.value))) + else: + westl = "E" + str(np.absolute(int(np.round(maxxlon.value)))) + + if centlat >= 0.0: + centlat = "N" + str(int(np.round(centlat.value))) + else: + centlat = "S" + str(np.absolute(int(np.round(centlat.value)))) + if centlon >= 0.0: + centlon = "W" + str(int(np.round(centlon.value))) + else: + centlon = "E" + str(np.absolute(int(np.round(centlon.value)))) + + """insertions of CH properties into property array""" + + ins_prop( + datm, + rs, + ident, + props, + arcar, + arccent, + pos, + npix, + trummar, + centlat, + centlon, + mB, + mBpos, + mBneg, + Ywb, + Xwb, + Yeb, + Xeb, + Ynb, + Xnb, + Ysb, + Xsb, + width, + eastl, + westl, + ) + """sets up code for next possible coronal hole""" + + ident = ident + 1 + +"""sets ident back to max value of iarr""" + +ident = ident - 1 + +"""stores all CH properties in a text file""" +np.savetxt("ch_summary.txt", props, fmt="%s") def rescale01(arr, cmin=None, cmax=None, a=0, b=1): + """ + Rescales array + + Parameters + ---------- + arr: 'np.arr' + cmin: 'np.float' + cmax: 'np.float' + a: 'int' + b: 'int' + + Returns + ------- + np.array + + """ if cmin or cmax: arr = np.clip(arr, cmin, cmax) return (b - a) * ((arr - np.min(arr)) / (np.max(arr) - np.min(arr))) + a -def plot_tricolor(data, datb, datc, xgrid, ygrid, slate): +def plot_tricolor(): + """ + Plots a tricolor mask of image data + + Returns + ------- + plot: 'matplotlib.image.AxesImage' + + """ + tricolorarray = np.zeros((4096, 4096, 3)) data_a = img_as_ubyte(rescale01(np.log10(data), cmin=1.2, cmax=3.9)) @@ -491,17 +950,30 @@ def plot_tricolor(data, datb, datc, xgrid, ygrid, slate): fig, ax = plt.subplots(figsize=(10, 10)) - plt.imshow(tricolorarray, origin="lower") # , extent = ) + plt.imshow(tricolorarray, origin="lower") plt.contour(xgrid, ygrid, slate, colors="white", linewidths=0.5) plt.savefig("tricolor.png") plt.close() -def plot_mask(slate, iarr, circ, rs, dattoarc, center, xgrid, ygrid, hedb): +def plot_mask(slate=slate): + """ + Plots the contour mask + + Parameters + ---------- + slate: 'np.array' + + Returns + ------- + plot: 'matplotlib.image.AxesImage' + + """ + chs = np.where(iarr > 0) slate[chs] = 1 slate = np.array(slate, dtype=np.uint8) - # cont, heir = cv2.findContours(slate, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) + cont, heir = cv2.findContours(slate, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) circ[:] = 0 r = rs / dattoarc @@ -518,4 +990,7 @@ def plot_mask(slate, iarr, circ, rs, dattoarc, center, xgrid, ygrid, hedb): plt.contour(xgrid, ygrid, circ, colors="black", linewidths=1.0) plt.savefig("CH_mask_" + hedb["DATE"] + ".png", transparent=True) - plt.close() + + +plot_tricolor() +plot_mask() diff --git a/chimerapy/chimera_copy.py b/chimerapy/chimera_copy.py new file mode 100644 index 0000000..1e4bc71 --- /dev/null +++ b/chimerapy/chimera_copy.py @@ -0,0 +1,160 @@ +import copy + + +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.colors import LogNorm +import sunpy +import sunpy.map +import astropy.units as u +from astropy.coordinates import SkyCoord +from astropy.units import UnitsError +from astropy.visualization import make_lupton_rgb +from sunpy.map import Map, all_coordinates_from_map, pixelate_coord_path, sample_at_coords + +m171 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0171.fits") +m193 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0193.fits") +m211 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0211.fits") + +def segmenting_plots(scale_cold: float, scale_warm: float, scale_hot: float): + m171 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0171.fits") + m193 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0193.fits") + m211 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0211.fits") + + m171_orig = copy.deepcopy(m171) + m193_orig = copy.deepcopy(m193) + m211_orig = copy.deepcopy(m211) + + # Since the data are taken at similar times neglect any coordinate changes so just use 171 maps coordinates + coords = all_coordinates_from_map(m171) + disk_mask = (coords.Tx**2 + coords.Ty**2) ** 0.5 < m171.rsun_obs + + m171 = m171 / m171.exposure_time + m193 = m193 / m193.exposure_time + m211 = m211 / m211.exposure_time + + xx = np.linspace(0, 300, 500) + + fig, axes = plt.subplot_mosaic( + [["cool_hist"], ["warm_hist"], ["hot_hist"]], + layout="constrained", + figsize=(3, 5), + ) + + # 171 v 193 + cool_counts, *cool_bins = axes["cool_hist"].hist2d( + m171.data[disk_mask].flatten(), + m193.data[disk_mask].flatten(), + bins=150, + range=[[0, 300], [0, 300]], + norm=LogNorm(), + density=True, + ) + #Finding the indices of nonzero counts + non_zero_cool = np.where(cool_counts > 0) + + #Finding the corresponding minimum y-index + min_y_index = np.min(non_zero_cool[1]) + + #Map the index to the actual y-value + min_y_cool= cool_bins[1][min_y_index] + axes["cool_hist"].set_facecolor("k") + axes["cool_hist"].plot(xx, (xx**scale_cold) + min_y_cool, "w") + + # 171 v 211 + warm_counts, *warm_bins = axes["warm_hist"].hist2d( + m171.data[disk_mask].flatten(), + m211.data[disk_mask].flatten(), + bins=250, + range=[[0, 300], [0, 300]], + norm=LogNorm(), + density=True, + ) + #Finding the indices of nonzero counts + non_zero_warm = np.where(warm_counts > 0) + + #Finding the corresponding minimum y-index + min_y_index = np.min(non_zero_warm[1]) + + #Map the index to the actual y-value + min_y_warm= warm_bins[1][min_y_index] + axes["warm_hist"].set_ylim(0, 100) + axes["warm_hist"].set_facecolor("k") + axes["warm_hist"].plot(xx, (xx**scale_warm) + min_y_warm, "w") + + # 193 v 311 + hot_counts, *hot_bins = axes["hot_hist"].hist2d( + m193.data[disk_mask].flatten(), + m211.data[disk_mask].flatten(), + bins=250, + range=[[0, 300], [0, 300]], + norm=LogNorm(), + density=True, + ) + #Finding the indices of nonzero counts + non_zero_hot = np.where(hot_counts > 0) + + #Finding the corresponding minimum y-index + min_y_index = np.min(non_zero_hot[1]) + + #Map the index to the actual y-value + min_y_hot= hot_bins[1][min_y_index] + + axes["hot_hist"].set_ylim(0, 100) + axes["hot_hist"].set_facecolor("k") + axes["hot_hist"].plot(xx, (xx**scale_hot) + min_y_hot, "w") + + plt.show() + +segmenting_plots(.7, .6, .7) + +def clip_scale(map1: sunpy.map, clip1: float, clip2: float, clip3: float, scale: float): + map_clipped = np.clip(np.log10(map1.data), clip1, clip3) + map_clipped_scaled = ((map_clipped - clip1) / clip2) * scale + return map_clipped_scaled + +cs_171 = clip_scale(m171, 1.2, 2.7, 3.9, 255) +cs_193 = clip_scale(m193, 1.4, 1.6, 3.0, 255) +cs_211 = clip_scale(m211, .8, 1.9, 2.7, 255) + + +def create_mask(scale1: float, scale2: float, scale3: float): + mask1 = (cs_171 / cs_211) >= (np.mean(m171.data) * scale1) / np.mean(m211.data) + mask2 = (cs_211 + cs_193) < (scale2 * (np.mean(m193.data) + np.mean(m211.data))) + mask3 = (cs_171 / cs_193) >= ((np.mean(m171.data) * scale3) / np.mean(m193.data)) + return mask1, mask2, mask3 + +mask211_171, mask211_193, mask171_193 = create_mask(.6357, .7, 1.5102) + +tri_color_img = make_lupton_rgb(cs_171, cs_193, cs_211, Q=10, stretch=50) +comb_mask = mask211_171 * mask211_193 * mask171_193 + + +fig, axes = plt.subplot_mosaic([["tri", "comb_mask"]], layout="constrained", figsize=(6, 3)) + +axes["tri"].imshow(tri_color_img, origin="lower") +axes["comb_mask"].imshow(comb_mask, origin="lower") + +plt.show() + +def create_contours(mask1: np.ndarray, mask2: np.ndarray, mask3: np.ndarray): + mask_map = Map(((mask1 * mask2 * mask3).astype(int), m171.meta)) + try: + contours = mask_map.contour(0.5 / u.s) + except UnitsError: + contours = mask_map.contour(50 * u.percent) + + contours = sorted(contours, key=lambda x: x.size, reverse=True) + + fig, axes = plt.subplot_mosaic( + [["seg"]], layout="constrained", figsize=(6, 3), subplot_kw={"projection": m171} + ) + m171.plot(axes=axes["seg"]) + axes["seg"].imshow(tri_color_img) + + # For the moment just plot to top 5 contours based on "size" for contour + for contour in contours[:6]: + axes["seg"].plot_coord(contour, color="w", linewidth=0.5) + plt.show() + +create_contours(mask211_171, mask211_193, mask171_193) \ No newline at end of file diff --git a/chimerapy/chimera_legacy.py b/chimerapy/chimera_legacy.py new file mode 100644 index 0000000..4462553 --- /dev/null +++ b/chimerapy/chimera_legacy.py @@ -0,0 +1,521 @@ +""" + +""" + + +import glob +import sys + +import astropy.units as u +import cv2 +import mahotas +import matplotlib.pyplot as plt +import numpy as np +import scipy.interpolate +import sunpy.map +from astropy import wcs +from astropy.io import fits +from astropy.modeling.models import Gaussian2D +from skimage.util import img_as_ubyte + + +def chimera_legacy(): + file_path = "./" + + 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") + + circ, data, datb, datc, dattoarc, hedb, iarr, props, rs, slate, center, xgrid, ygrid = chimera( + im171, im193, im211, imhmi + ) + + # =====sets ident back to max value of iarr====== + + # ident = ident - 1 + np.savetxt("ch_summary.txt", props, fmt="%s") + + plot_tricolor(data, datb, datc, xgrid, ygrid, slate) + plot_mask(slate, iarr, circ, rs, dattoarc, center, xgrid, ygrid, hedb) + + +def chimera(im171, im193, im211, imhmi): + if im171 == [] or im193 == [] or im211 == [] or imhmi == []: + print("Not all required files present") + sys.exit() + # =====Reads in data and resizes images===== + x = np.arange(0, 1024) * 4 + hdu_number = 0 + heda = fits.getheader(im171[0], hdu_number) + data = fits.getdata(im171[0], ext=0) / (heda["EXPTIME"]) + dn = scipy.interpolate.interp2d(x, x, data) + data = dn(np.arange(0, 4096), np.arange(0, 4096)) + hedb = fits.getheader(im193[0], hdu_number) + datb = fits.getdata(im193[0], ext=0) / (hedb["EXPTIME"]) + dn = scipy.interpolate.interp2d(x, x, datb) + datb = dn(np.arange(0, 4096), np.arange(0, 4096)) + hedc = fits.getheader(im211[0], hdu_number) + datc = fits.getdata(im211[0], ext=0) / (hedc["EXPTIME"]) + dn = scipy.interpolate.interp2d(x, x, datc) + datc = dn(np.arange(0, 4096), np.arange(0, 4096)) + hedm = fits.getheader(imhmi[0], hdu_number) + 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) + if hedm["crota1"] > 90: + datm = np.rot90(np.rot90(datm)) + # =====Specifies solar radius and calculates conversion value of pixel to arcsec===== + s = np.shape(data) + rs = heda["rsun"] + if hedb["ctype1"] != "solar_x ": + hedb["ctype1"] = "solar_x " + hedb["ctype2"] = "solar_y " + if heda["cdelt1"] > 1: + heda["cdelt1"], heda["cdelt2"], heda["crpix1"], heda["crpix2"] = ( + heda["cdelt1"] / 4.0, + heda["cdelt2"] / 4.0, + heda["crpix1"] * 4.0, + heda["crpix2"] * 4.0, + ) + hedb["cdelt1"], hedb["cdelt2"], hedb["crpix1"], hedb["crpix2"] = ( + hedb["cdelt1"] / 4.0, + hedb["cdelt2"] / 4.0, + hedb["crpix1"] * 4.0, + hedb["crpix2"] * 4.0, + ) + hedc["cdelt1"], hedc["cdelt2"], hedc["crpix1"], hedc["crpix2"] = ( + hedc["cdelt1"] / 4.0, + hedc["cdelt2"] / 4.0, + hedc["crpix1"] * 4.0, + hedc["crpix2"] * 4.0, + ) + dattoarc = heda["cdelt1"] + convermul = dattoarc / hedm["cdelt1"] + # =====Alternative coordinate systems===== + hdul = fits.open(im171[0]) + hdul[0].header["CUNIT1"] = "arcsec" + hdul[0].header["CUNIT2"] = "arcsec" + aia = sunpy.map.Map(hdul[0].data, hdul[0].header) + adj = 4096.0 / aia.dimensions[0].value + x, y = (np.meshgrid(*[np.arange(adj * v.value) for v in aia.dimensions]) * u.pixel) / adj + hpc = aia.pixel_to_world(x, y) + hg = hpc.transform_to(sunpy.coordinates.frames.HeliographicStonyhurst) + csys = wcs.WCS(hedb) + # =======setting up arrays to be used============ + ident = 1 + iarr = np.zeros((s[0], s[1]), dtype=np.byte) + offarr, slate = np.array(iarr), np.array(iarr) + bmcool = np.zeros((s[0], s[1]), dtype=np.float32) + cand, bmmix, bmhot = np.array(bmcool), np.array(bmcool), np.array(bmcool) + circ = np.zeros((s[0], s[1]), dtype=int) + # =======creation of a 2d gaussian for magnetic cut offs=========== + r = (s[1] / 2.0) - 450 + xgrid, ygrid = np.meshgrid(np.arange(s[0]), np.arange(s[1])) + center = [int(s[1] / 2.0), int(s[1] / 2.0)] + w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 > r**2) + y, x = np.mgrid[0:4096, 0:4096] + garr = Gaussian2D(1, s[0] / 2, s[1] / 2, 2000 / 2.3548, 2000 / 2.3548)(x, y) + garr[w] = 1.0 + # ======creation of array for CH properties========== + props = np.zeros((26, 30), dtype="", + "", + "", + "BMAX", + "BMIN", + "TOT_B+", + "TOT_B-", + "", + "", + "", + ) + props[:, 1] = ( + "num", + '"', + '"', + "H deg", + '"', + '"', + '"', + '"', + '"', + '"', + '"', + '"', + "H deg", + "deg", + "Mm^2", + "%", + "G", + "G", + "G", + "G", + "G", + "G", + "G", + "Mx", + "Mx", + "Mx", + ) + # =====removes negative data values===== + data[np.where(data <= 0)] = 0 + datb[np.where(datb <= 0)] = 0 + datc[np.where(datc <= 0)] = 0 + # ============make a multi-wavelength image for contours================== + with np.errstate(divide="ignore"): + t0 = np.log10(datc) + t1 = np.log10(datb) + t2 = np.log10(data) + t0[np.where(t0 < 0.8)] = 0.8 + t0[np.where(t0 > 2.7)] = 2.7 + t1[np.where(t1 < 1.4)] = 1.4 + t1[np.where(t1 > 3.0)] = 3.0 + t2[np.where(t2 < 1.2)] = 1.2 + t2[np.where(t2 > 3.9)] = 3.9 + t0 = np.array(((t0 - 0.8) / (2.7 - 0.8)) * 255, dtype=np.float32) + t1 = np.array(((t1 - 1.4) / (3.0 - 1.4)) * 255, dtype=np.float32) + t2 = np.array(((t2 - 1.2) / (3.9 - 1.2)) * 255, dtype=np.float32) + # ====create 3 segmented bitmasks===== + with np.errstate(divide="ignore", invalid="ignore"): + bmmix[np.where(t2 / t0 >= ((np.mean(data) * 0.6357) / (np.mean(datc))))] = 1 + bmhot[np.where(t0 + t1 < (0.7 * (np.mean(datb) + np.mean(datc))))] = 1 + bmcool[np.where(t2 / t1 >= ((np.mean(data) * 1.5102) / (np.mean(datb))))] = 1 + # ====logical conjunction of 3 segmentations======= + cand = bmcool * bmmix * bmhot + # ====plot tricolour image with lon/lat conotours======= + # ======removes off detector mis-identifications========== + r = (s[1] / 2.0) - 100 + w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 <= r**2) + circ[w] = 1.0 + cand = cand * circ + # =======Separates on-disk and off-limb CHs=============== + circ[:] = 0 + r = (rs / dattoarc) - 10 + w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 <= r**2) + circ[w] = 1.0 + r = (rs / dattoarc) + 40 + w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 >= r**2) + circ[w] = 1.0 + cand = cand * circ + # ====open file for property storage===== + # =====contours the identified datapoints======= + cand = np.array(cand, dtype=np.uint8) + cont, heir = cv2.findContours(cand, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) + # ======sorts contours by size============ + sizes = [] + for i in range(len(cont)): + sizes = np.append(sizes, len(cont[i])) + reord = sizes.ravel().argsort()[::-1] + tmp = list(cont) + for i in range(len(cont)): + tmp[i] = cont[reord[i]] + cont = list(tmp) + # =====cycles through contours========= + for i in range(len(cont)): + x = np.append(x, len(cont[i])) + + # =====only takes values of minimum surface length and calculates area====== + + if len(cont[i]) <= 100: + continue + area = 0.5 * np.abs( + np.dot(cont[i][:, 0, 0], np.roll(cont[i][:, 0, 1], 1)) + - np.dot(cont[i][:, 0, 1], np.roll(cont[i][:, 0, 0], 1)) + ) + arcar = area * (dattoarc**2) + if arcar > 1000: + # =====finds centroid======= + + # chpts = len(cont[i]) + cent = [np.mean(cont[i][:, 0, 0]), np.mean(cont[i][:, 0, 1])] + + # ===remove quiet sun regions encompassed by coronal holes====== + + if ( + cand[ + np.max(cont[i][:, 0, 0]) + 1, + cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], + ] + > 0 + ) and ( + iarr[ + np.max(cont[i][:, 0, 0]) + 1, + cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], + ] + > 0 + ): + mahotas.polygon.fill_polygon(np.array(list(zip(cont[i][:, 0, 1], cont[i][:, 0, 0]))), slate) + iarr[np.where(slate == 1)] = 0 + slate[:] = 0 + + else: + # ====create a simple centre point====== + + arccent = csys.all_pix2world(cent[0], cent[1], 0) + + # ====classifies off limb CH regions======== + + if (((arccent[0] ** 2) + (arccent[1] ** 2)) > (rs**2)) or ( + np.sum(np.array(csys.all_pix2world(cont[i][0, 0, 0], cont[i][0, 0, 1], 0)) ** 2) > (rs**2) + ): + mahotas.polygon.fill_polygon( + np.array(list(zip(cont[i][:, 0, 1], cont[i][:, 0, 0]))), offarr + ) + else: + # =====classifies on disk coronal holes======= + + mahotas.polygon.fill_polygon( + np.array(list(zip(cont[i][:, 0, 1], cont[i][:, 0, 0]))), slate + ) + poslin = np.where(slate == 1) + slate[:] = 0 + + # ====create an array for magnetic polarity======== + + pos = np.zeros((len(poslin[0]), 2), dtype=np.uint) + pos[:, 0] = np.array((poslin[0] - (s[0] / 2)) * convermul + (s[1] / 2), dtype=np.uint) + pos[:, 1] = np.array((poslin[1] - (s[0] / 2)) * convermul + (s[1] / 2), dtype=np.uint) + npix = list( + np.histogram( + datm[pos[:, 0], pos[:, 1]], + bins=np.arange( + np.round(np.min(datm[pos[:, 0], pos[:, 1]])) - 0.5, + np.round(np.max(datm[pos[:, 0], pos[:, 1]])) + 0.6, + 1, + ), + ) + ) + npix[0][np.where(npix[0] == 0)] = 1 + npix[1] = npix[1][:-1] + 0.5 + + wh1 = np.where(npix[1] > 0) + wh2 = np.where(npix[1] < 0) + + # =====magnetic cut offs dependent on area========= + + if ( + np.absolute((np.sum(npix[0][wh1]) - np.sum(npix[0][wh2])) / np.sqrt(np.sum(npix[0]))) + <= 10 + and arcar < 9000 + ): + continue + if ( + np.absolute(np.mean(datm[pos[:, 0], pos[:, 1]])) < garr[int(cent[0]), int(cent[1])] + and arcar < 40000 + ): + continue + iarr[poslin] = ident + + # ====create an accurate center point======= + + ypos = np.sum((poslin[0]) * np.absolute(hg.lat[poslin])) / np.sum( + np.absolute(hg.lat[poslin]) + ) + xpos = np.sum((poslin[1]) * np.absolute(hg.lon[poslin])) / np.sum( + np.absolute(hg.lon[poslin]) + ) + + arccent = csys.all_pix2world(xpos, ypos, 0) + + # ======calculate average angle coronal hole is subjected to====== + + dist = np.sqrt((arccent[0] ** 2) + (arccent[1] ** 2)) + ang = np.arcsin(dist / rs) + + # =====calculate area of CH with minimal projection effects====== + + trupixar = abs(area / np.cos(ang)) + truarcar = trupixar * (dattoarc**2) + trummar = truarcar * ((6.96e08 / rs) ** 2) + + # ====find CH extent in latitude and longitude======== + + # maxxlat = hg.lat[ + # cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], + # np.max(cont[i][:, 0, 0]), + # ] + maxxlon = hg.lon[ + cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], + np.max(cont[i][:, 0, 0]), + ] + # maxylat = hg.lat[ + # np.max(cont[i][:, 0, 1]), + # cont[i][np.where(cont[i][:, 0, 1] == np.max(cont[i][:, 0, 1]))[0][0], 0, 0], + # ] + # maxylon = hg.lon[ + # np.max(cont[i][:, 0, 1]), + # cont[i][np.where(cont[i][:, 0, 1] == np.max(cont[i][:, 0, 1]))[0][0], 0, 0], + # ] + # minxlat = hg.lat[ + # cont[i][np.where(cont[i][:, 0, 0] == np.min(cont[i][:, 0, 0]))[0][0], 0, 1], + # np.min(cont[i][:, 0, 0]), + # ] + minxlon = hg.lon[ + cont[i][np.where(cont[i][:, 0, 0] == np.min(cont[i][:, 0, 0]))[0][0], 0, 1], + np.min(cont[i][:, 0, 0]), + ] + # minylat = hg.lat[ + # np.min(cont[i][:, 0, 1]), + # cont[i][np.where(cont[i][:, 0, 1] == np.min(cont[i][:, 0, 1]))[0][0], 0, 0], + # ] + # minylon = hg.lon[ + # np.min(cont[i][:, 0, 1]), + # cont[i][np.where(cont[i][:, 0, 1] == np.min(cont[i][:, 0, 1]))[0][0], 0, 0], + # ] + + # =====CH centroid in lat/lon======= + + centlat = hg.lat[int(ypos), int(xpos)] + centlon = hg.lon[int(ypos), int(xpos)] + + # ====calculate the mean magnetic field===== + + mB = np.mean(datm[pos[:, 0], pos[:, 1]]) + mBpos = np.sum(npix[0][wh1] * npix[1][wh1]) / np.sum(npix[0][wh1]) + mBneg = np.sum(npix[0][wh2] * npix[1][wh2]) / np.sum(npix[0][wh2]) + + # =====finds coordinates of CH boundaries======= + + Ywb, Xwb = csys.all_pix2world( + cont[i][np.where(cont[i][:, 0, 0] == np.max(cont[i][:, 0, 0]))[0][0], 0, 1], + np.max(cont[i][:, 0, 0]), + 0, + ) + Yeb, Xeb = csys.all_pix2world( + cont[i][np.where(cont[i][:, 0, 0] == np.min(cont[i][:, 0, 0]))[0][0], 0, 1], + np.min(cont[i][:, 0, 0]), + 0, + ) + Ynb, Xnb = csys.all_pix2world( + np.max(cont[i][:, 0, 1]), + cont[i][np.where(cont[i][:, 0, 1] == np.max(cont[i][:, 0, 1]))[0][0], 0, 0], + 0, + ) + Ysb, Xsb = csys.all_pix2world( + np.min(cont[i][:, 0, 1]), + cont[i][np.where(cont[i][:, 0, 1] == np.min(cont[i][:, 0, 1]))[0][0], 0, 0], + 0, + ) + + width = round(maxxlon.value) - round(minxlon.value) + + if minxlon.value >= 0.0: + eastl = "W" + str(int(np.round(minxlon.value))) + else: + eastl = "E" + str(np.absolute(int(np.round(minxlon.value)))) + if maxxlon.value >= 0.0: + westl = "W" + str(int(np.round(maxxlon.value))) + else: + westl = "E" + str(np.absolute(int(np.round(maxxlon.value)))) + + if centlat >= 0.0: + centlat = "N" + str(int(np.round(centlat.value))) + else: + centlat = "S" + str(np.absolute(int(np.round(centlat.value)))) + if centlon >= 0.0: + centlon = "W" + str(int(np.round(centlon.value))) + else: + centlon = "E" + str(np.absolute(int(np.round(centlon.value)))) + + # ====insertions of CH properties into property array===== + + props[0, ident + 1] = str(ident) + props[1, ident + 1] = str(np.round(arccent[0])) + props[2, ident + 1] = str(np.round(arccent[1])) + props[3, ident + 1] = str(centlon + centlat) + props[4, ident + 1] = str(np.round(Xeb)) + props[5, ident + 1] = str(np.round(Yeb)) + props[6, ident + 1] = str(np.round(Xwb)) + props[7, ident + 1] = str(np.round(Ywb)) + props[8, ident + 1] = str(np.round(Xnb)) + props[9, ident + 1] = str(np.round(Ynb)) + props[10, ident + 1] = str(np.round(Xsb)) + props[11, ident + 1] = str(np.round(Ysb)) + props[12, ident + 1] = str(eastl + "-" + westl) + props[13, ident + 1] = str(width) + props[14, ident + 1] = f"{trummar / 1e+12:.1e}" + props[15, ident + 1] = str(np.round((arcar * 100 / (np.pi * (rs**2))), 1)) + props[16, ident + 1] = str(np.round(mB, 1)) + props[17, ident + 1] = str(np.round(mBpos, 1)) + props[18, ident + 1] = str(np.round(mBneg, 1)) + props[19, ident + 1] = str(np.round(np.max(npix[1]), 1)) + props[20, ident + 1] = str(np.round(np.min(npix[1]), 1)) + tbpos = np.sum(datm[pos[:, 0], pos[:, 1]][np.where(datm[pos[:, 0], pos[:, 1]] > 0)]) + props[21, ident + 1] = f"{tbpos:.1e}" + tbneg = np.sum(datm[pos[:, 0], pos[:, 1]][np.where(datm[pos[:, 0], pos[:, 1]] < 0)]) + props[22, ident + 1] = f"{tbneg:.1e}" + props[23, ident + 1] = f"{mB * trummar * 1e+16:.1e}" + props[24, ident + 1] = f"{mBpos * trummar * 1e+16:.1e}" + props[25, ident + 1] = f"{mBneg * trummar * 1e+16:.1e}" + + # =====sets up code for next possible coronal hole===== + + ident = ident + 1 + return circ, data, datb, datc, dattoarc, hedb, iarr, props, rs, slate, center, xgrid, ygrid + + +def rescale01(arr, cmin=None, cmax=None, a=0, b=1): + if cmin or cmax: + arr = np.clip(arr, cmin, cmax) + return (b - a) * ((arr - np.min(arr)) / (np.max(arr) - np.min(arr))) + a + + +def plot_tricolor(data, datb, datc, xgrid, ygrid, slate): + tricolorarray = np.zeros((4096, 4096, 3)) + + data_a = img_as_ubyte(rescale01(np.log10(data), cmin=1.2, cmax=3.9)) + data_b = img_as_ubyte(rescale01(np.log10(datb), cmin=1.4, cmax=3.0)) + data_c = img_as_ubyte(rescale01(np.log10(datc), cmin=0.8, cmax=2.7)) + + tricolorarray[..., 0] = data_c / np.max(data_c) + tricolorarray[..., 1] = data_b / np.max(data_b) + tricolorarray[..., 2] = data_a / np.max(data_a) + + fig, ax = plt.subplots(figsize=(10, 10)) + + plt.imshow(tricolorarray, origin="lower") # , extent = ) + plt.contour(xgrid, ygrid, slate, colors="white", linewidths=0.5) + plt.savefig("tricolor.png") + plt.close() + + +def plot_mask(slate, iarr, circ, rs, dattoarc, center, xgrid, ygrid, hedb): + chs = np.where(iarr > 0) + slate[chs] = 1 + slate = np.array(slate, dtype=np.uint8) + # cont, heir = cv2.findContours(slate, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) + + circ[:] = 0 + r = rs / dattoarc + w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 <= r**2) + circ[w] = 1.0 + + plt.figure(figsize=(10, 10)) + plt.xlim(143, 4014) + plt.ylim(143, 4014) + plt.scatter(chs[1], chs[0], marker="s", s=0.0205, c="black", cmap="viridis", edgecolor="none", alpha=0.2) + plt.gca().set_aspect("equal", adjustable="box") + plt.axis("off") + plt.contour(xgrid, ygrid, slate, colors="black", linewidths=0.5) + plt.contour(xgrid, ygrid, circ, colors="black", linewidths=1.0) + + plt.savefig("CH_mask_" + hedb["DATE"] + ".png", transparent=True) + plt.close() diff --git a/chimerapy/tests/test_chimera.py b/chimerapy/tests/test_chimera.py index e08da98..714c112 100644 --- a/chimerapy/tests/test_chimera.py +++ b/chimerapy/tests/test_chimera.py @@ -1,9 +1,20 @@ import os -from pathlib import Path +import warnings -from parfive import Downloader +import mahotas +import numpy as np +from astropy.io import fits +from astropy.modeling.models import Gaussian2D -from chimerapy.chimera import chimera_legacy +from chimerapy.chimera import * +from chimerapy.chimera import (Bounds, Xeb, Xnb, Xsb, Xwb, Yeb, Ynb, Ysb, Ywb, + ang, arcar, arccent, area, cent, centlat, + centlon, chpts, cont, coords, csys, data, datb, + datc, datm, dist, eastl, extent, filter, hg, + ins_prop, mB, mBneg, mBpos, npix, pos, + remove_neg, rescale_aia, rescale_hmi, + set_contour, sort, threshold, truarcar, trummar, + trupixar, westl, width, xpos, ypos) INPUT_FILES = { "aia171": "https://solarmonitor.org/data/2016/09/22/fits/saia/saia_00171_fd_20160922_103010.fts.gz", @@ -13,17 +24,398 @@ } -def test_chimera(tmp_path): - Downloader.simple_download(INPUT_FILES.values(), path=tmp_path) - os.chdir(tmp_path) - chimera_legacy() +def img_present(): + assert im171 != [] or im193 != [] or im211 != [] or imhmi != [], "Not all required files present" - test_summary_file = Path(__file__).parent / "test_ch_summary.txt" - with test_summary_file.open("r") as f: - test_summary_text = f.read() - ch_summary_file = tmp_path / "ch_summary.txt" - with ch_summary_file.open("r") as f: - ch_summary_text = f.read() +def rest_rescale(): + global data, datb, datc, datm + data = rescale_aia(im171, 1024, 4096) + datb = rescale_aia(im193, 1024, 4096) + datc = rescale_aia(im211, 1024, 4096) + datm = rescale_hmi(imhmi, 1024, 4096) + assert ( + len(data) == len(datb) == len(datc) == len(datm) == 4096 + ), "Array size does not match desired array size" - assert ch_summary_text == test_summary_text + +heda = fits.getheader(im171[0], 0) +hedb = fits.getheader(im193[0], 0) +hedc = fits.getheader(im211[0], 0) +hedm = fits.getheader(imhmi[0], 0) + + +def test_filter(): + initial_values = { + 'hedb["ctyple1"]': hedb["ctype1"], + 'hedb["ctype2"]': hedb["ctype2"], + 'heda["cdelt1"]': heda["cdelt1"], + 'heda["cdelt2"]': heda["cdelt2"], + 'heda["crpix1"]': heda["crpix1"], + 'heda["crpix2"]': heda["crpix2"], + 'hedb["cdelt1"]': hedb["cdelt1"], + 'hedb["cdelt2"]': hedb["cdelt2"], + 'hedb["crpix1"]': hedb["crpix1"], + 'hedb["crpix2"]': hedb["crpix2"], + 'hedc["cdelt1"]': hedc["cdelt1"], + 'hedc["cdelt2"]': hedc["cdelt2"], + 'hedc["crpix1"]': hedc["crpix1"], + 'hedc["crpix2"]': hedc["crpix2"], + "datm": datm, + } + filter(im171, im193, im211, imhmi) + final_values = { + 'hedb["ctyple1"]': hedb["ctype1"], + 'hedb["ctype2"]': hedb["ctype2"], + 'heda["cdelt1"]': heda["cdelt1"], + 'heda["cdelt2"]': heda["cdelt2"], + 'heda["crpix1"]': heda["crpix1"], + 'heda["crpix2"]': heda["crpix2"], + 'hedb["cdelt1"]': hedb["cdelt1"], + 'hedb["cdelt2"]': hedb["cdelt2"], + 'hedb["crpix1"]': hedb["crpix1"], + 'hedb["crpix2"]': hedb["crpix2"], + 'hedc["cdelt1"]': hedc["cdelt1"], + 'hedc["cdelt2"]': hedc["cdelt2"], + 'hedc["crpix1"]': hedc["crpix1"], + 'hedc["crpix2"]': hedc["crpix2"], + "datm": datm, + } + if initial_values == final_values: + raise Warning("No filtering occured - ensure filtering conditions were not met") + + +def test_neg(): + remove_neg(im171, im193, im211) + for num in im171: + assert num >= 0, "Array still contains negative number" + for num in im193: + assert num >= 0, "Array still contains negative number" + for num in im211: + assert num >= 0, "Array still contains negative number" + + +def tremove_neg(): + assert len(data[data < 0]) == 0, "Data contains negative values" + assert len(datb[datb < 0]) == 0, "Data contains negative values" + assert len(datc[datc < 0]) == 0, "Data contains negative values" + + +def s_rs(): + global s, rs + s = np.shape(data) + rs = heda["rsun"] + assert s == (4096, 4096), "Incorrect data shape" + if rs < 970 or rs > 975: + warnings.warn("Solar radius may be inconsistant with accepted value (~973)") + + +ident = 1 +iarr = np.zeros((s[0], s[1]), dtype=np.byte) +bmcool = np.zeros((s[0], s[1]), dtype=np.float32) +offarr, slate = np.array(iarr), np.array(iarr) +cand, bmmix, bmhot = np.array(bmcool), np.array(bmcool), np.array(bmcool) +circ = np.zeros((s[0], s[1]), dtype=int) + +r = (s[1] / 2.0) - 450 +xgrid, ygrid = np.meshgrid(np.arange(s[0]), np.arange(s[1])) +center = [int(s[1] / 2.0), int(s[1] / 2.0)] +w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 > r**2) +y, x = np.mgrid[0:4096, 0:4096] +garr = Gaussian2D(1, s[0] / 2, s[1] / 2, 2000 / 2.3548, 2000 / 2.3548)(x, y) +garr[w] = 1.0 + +props = np.zeros((26, 30), dtype="", + "", + "", + "BMAX", + "BMIN", + "TOT_B+", + "TOT_B-", + "", + "", + "", +) +props[:, 1] = ( + "num", + '"', + '"', + "H°", + '"', + '"', + '"', + '"', + '"', + '"', + '"', + '"', + "H°", + "°", + "Mm^2", + "%", + "G", + "G", + "G", + "G", + "G", + "G", + "G", + "Mx", + "Mx", + "Mx", +) + +assert len(props) == 26, "Incorrect property array size" + + +def test_bounds(): + t0b = Bounds(0.8, 2.7, 255) + t1b = Bounds(1.4, 3.0, 255) + t2b = Bounds(1.2, 3.9, 255) + assert t0b.upper > t0b.lower, "Upper bound must be greater than lower bound" + assert t1b.upper > t0b.lower, "Upper bound must be greater than lower bound" + assert t2b.upper > t0b.lower, "Upper bound must be greater than lower bound" + + +with np.errstate(divide="ignore"): + t0 = np.log10(datc) + t1 = np.log10(datb) + t2 = np.log10(data) + + +def test_threshold(): + init = {"t0": t0, "t1": t1, "t2": t2} + assert init["t0"] != threshold(t0), "Data was not bounded by threshold values" + assert init["t1"] != threshold(t0), "Data was not bounded by threshold values" + assert init["t2"] != threshold(t0), "Data was not bounded by threshold values" + + +t0b = Bounds(0.8, 2.7, 255) +t1b = Bounds(1.4, 3.0, 255) +t2b = Bounds(1.2, 3.9, 255) + + +def test_contour(): + init = {"t0": t0, "t1": t1, "t2": t2} + assert init["t0"] != set_contour(t0), "Data was not bounded by threshold values" + assert init["t1"] != set_contour(t0), "Data was not bounded by threshold values" + assert init["t2"] != set_contour(t0), "Data was not bounded by threshold values" + + +def has_dupl(arr): + seen = set() + for num in arr: + if num in seen: + return True + seen.add(num) + return False + + +def test_dupl(): + global sizes, reord, tmp, cont + sort() + assert not has_dupl(reord), "Sorted list should contain no duplicates" + assert not has_dupl(tmp), "Sorted list should contain no duplicates" + assert not has_dupl(cont), "Sorted list should contain no duplicates" + + +for i in range(len(cont)): + x = np.append(x, len(cont[i])) + + +def test_extent(): + i_maxxlat = None + i_maxxlon = None + i_maxylat = None + i_maxylon = None + i_minxlon = None + i_minylat = None + i_minylon = None + i_minxlat = None + extent(i, ypos, xpos, hg, cont) + global maxxlat, maxxlon, maxylat, maxylon, minxlon, minylat, minylon, minxlat + assert i_maxxlat != maxxlat, "maxxlat not created successfully" + assert i_maxxlon != maxxlon, "maxxlon not created successfully" + assert i_maxylat != maxylat, "maxylat not created successfully" + assert i_maxylon != maxylon, "maxylon not created successfully" + assert i_minxlon != minxlon, "minxlon not created successfully" + assert i_minylat != minylat, "minylat not created successfully" + assert i_minylon != minylon, "minylon not created successfully" + assert i_minxlat != minxlat, "minxlat not created successfully" + + +def test_coords(): + i_Ywb = None + i_Xwb = None + i_Yeb = None + i_Xeb = None + i_Ynb = None + i_Xnb = None + i_Ysb = None + i_Xsb = None + coords(i, csys, cont) + assert i_Ywb != Ywb, "Ywb not created successfully" + assert i_Xwb != Xwb, "Xwb not created successfully" + assert i_Yeb != Yeb, "Yeb not created successfully" + assert i_Xeb != Xeb, "Xeb not created successfully" + assert i_Ynb != Ynb, "Ynb not created successfully" + assert i_Xnb != Xnb, "Xnb not created successfully" + assert i_Ysb != Ysb, "Ysb not created successfully" + assert i_Xsb != Xsb, "Xsb not created successfully" + + +def test_props(): + ins_prop( + datm, + rs, + ident, + props, + arcar, + arccent, + pos, + npix, + trummar, + centlat, + centlon, + mB, + mBpos, + mBneg, + Ywb, + Xwb, + Yeb, + Xeb, + Ynb, + Xnb, + Ysb, + Xsb, + width, + eastl, + westl, + ) + assert np.any(props is None) is not False, "Property array should not contain empty entries" + assert np.all(props is None) is not False, "Property array should not be empty" + + +def test_loop_variables(): + for i in range(len(cont)): + assert len(cont[i]) > 100, "Contour length should be greater than 100" + + +def test_fill_polygon(): + cand = np.array([[0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]]) + + iarr = np.array([[0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]]) + + cont = np.array([[[[1, 2], [1, 3], [2, 3], [2, 2]]]]) + + slate = np.array(iarr) + + polygon_vertices = np.array(list(zip(cont[0][:, 0, 1], cont[0][:, 0, 0]))) + + if ( + cand[ + np.max(cont[0][:, 0, 0]) + 1, + cont[0][np.where(cont[0][:, 0, 0] == np.max(cont[0][:, 0, 0]))[0][0], 0, 1], + ] + > 0 + ) and ( + iarr[ + np.max(cont[0][:, 0, 0]) + 1, + cont[0][np.where(cont[0][:, 0, 0] == np.max(cont[0][:, 0, 0]))[0][0], 0, 1], + ] + > 0 + ): + mahotas.polygon.fill_polygon(polygon_vertices, slate) + + print("After filling polygon:") + print(slate) + iarr[np.where(slate == 1)] = 0 + slate[:] = 0 + + assert np.array_equal( + slate, + np.array([[0, 0, 0, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]]), + ), "Test failed: slate array does not match expected result" + + assert np.array_equal( + iarr, + np.array([[0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]]), + ), "Test failed: iarr array does not match expected result" + + +def test_limb(): + arccent = csys.all_pix2world(cent[0], cent[1], 0) + if (((arccent[0] ** 2) + (arccent[1] ** 2)) > (rs**2)) or ( + np.sum(np.array(csys.all_pix2world(cont[i][0, 0, 0], cont[i][0, 0, 1], 0)) ** 2) > (rs**2) + ): + mahotas.polygon.fill_polygon(np.array(list(zip(cont[i][:, 0, 1], cont[i][:, 0, 0]))), offarr) + assert np.sum(offarr) > 0, "Offarr was not modified successfully" + else: + mahotas.polygon.fill_polygon(np.array(list(zip(cont[i][:, 0, 1], cont[i][:, 0, 0]))), slate) + slate[:] = 0 + assert np.sum(slate) > 0, "Slate was not modified successfully" + + +def test_magpol(): + assert npix[0][np.where(npix[0] == 0)] != 0, "Npix[0] should not be equal to zero at its zeros" + assert npix[1] != 0, "Npix[1] should not be equal to 0" + npixtest = [ + np.array([2, -1, 0, 3, -2, 1]), + np.array([1, -1, 0, 2, -2, 1]), + ] + wh1_expected = np.where(npixtest[1] > 0) + wh1_actual = np.where(npixtest[1] > 0) + + wh2_expected = np.where(npixtest[1] < 0) + wh2_actual = np.where(npixtest[1] < 0) + + assert np.array_equal( + wh1_actual, wh1_expected + ), f"Test failed for wh1: Expected {wh1_expected}, got {wh1_actual}" + assert np.array_equal( + wh2_actual, wh2_expected + ), f"Test failed for wh1: Expected {wh2_expected}, got {wh2_actual}" + + +assert area is not None, "Area variable not created successfully" +assert arcar is not None, "Arcar variable not created successfully" +assert chpts is not None, "Chpts variable not created successfully" +assert cent is not None, "Cent variable not created successfully" +assert arccent is not None, "Arccent variable not created successfully" +assert ypos is not None, "Ypos variable not created successfully" +assert xpos is not None, "Xpos variable not created suffessfully" +assert dist is not None, "Dist variable not created successfully" +assert ang is not None, "Ang variable not created successfully" +assert trupixar is not None, "Trupixar variable not created successfully" +assert truarcar is not None, "Truarcar variable not created successfully" +assert trummar is not None, "Trummar variable not created successfully" +assert mB is not None, "mB variable not created successfully" +assert mBpos is not None, "mBpos variable not created successfully" +assert mBneg is not None, "mBneg variable not created successfully" +assert width is not None, "width variable not created successfully" +assert eastl is not None, "eastl variable not created successfully" +assert westl is not None, "westl variable not created successfully" +assert centlat is not None, "centlat variable not created successfully" +assert centlon is not None, "centlon variable not created successfully" + +assert os.path.exists("ch_summary.txt"), "Summary file not saved correctly" +assert os.path.exists("tricolor.png"), "Tricolor image not saved correctly" +assert os.path.exists("CH_mask_" + hedb["DATE"] + ".png") diff --git a/chimerapy/tests/test_chimera_legacy.py b/chimerapy/tests/test_chimera_legacy.py new file mode 100644 index 0000000..41e30a3 --- /dev/null +++ b/chimerapy/tests/test_chimera_legacy.py @@ -0,0 +1,29 @@ +import os +from pathlib import Path + +from parfive import Downloader + +from chimerapy.chimera_legacy import chimera_legacy + +INPUT_FILES = { + "aia171": "https://solarmonitor.org/data/2016/09/22/fits/saia/saia_00171_fd_20160922_103010.fts.gz", + "aia193": "https://solarmonitor.org/data/2016/09/22/fits/saia/saia_00193_fd_20160922_103041.fts.gz", + "aia211": "https://solarmonitor.org/data/2016/09/22/fits/saia/saia_00211_fd_20160922_103046.fts.gz", + "hmi_mag": "https://solarmonitor.org/data/2016/09/22/fits/shmi/shmi_maglc_fd_20160922_094640.fts.gz", +} + + +def test_chimera(tmp_path): + Downloader.simple_download(INPUT_FILES.values(), path=tmp_path) + os.chdir(tmp_path) + chimera_legacy() + + test_summary_file = Path(__file__).parent / "test_ch_summary.txt" + with test_summary_file.open("r") as f: + test_summary_text = f.read() + + ch_summary_file = tmp_path / "ch_summary.txt" + with ch_summary_file.open("r") as f: + ch_summary_text = f.read() + + assert ch_summary_text == test_summary_text