diff --git a/CHIMERA_V2.py b/CHIMERA_V2.py deleted file mode 100644 index 48b0111..0000000 --- a/CHIMERA_V2.py +++ /dev/null @@ -1,686 +0,0 @@ -# import required libraries -import glob -import sys - -import astropy.units as u -import cv2 -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 - -plt.style.use(astropy_mpl_style) - -# loading in the images as fits files - -im171 = glob.glob("171.fts") -im193 = glob.glob("193.fts") -im211 = glob.glob("211.fts") -imhmi = glob.glob("hmi.fts") - -# ensure that all images are present - -if im171 == [] or im193 == [] or im211 == [] or imhmi == []: - print("Not all required files present") - sys.exit() - -# Two functions that rescale the aia and hmi images from any original size to any final size - -# didn't normalize by exposure time for hmi because it was equal to 0 - - -def rescale_aia(image: np.array, orig_res: int, desired_res: int): - hed = fits.getheader(image[0], 0) - dat = fits.getdata(image[0], 0) / (hed["EXPTIME"]) - 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, dat) - if len(dn(np.arange(0, desired_res), np.arange(0, desired_res))) != desired_res: - print("Incorrect image resolution") - sys.exit() - else: - 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, dat) - if len(dn(np.arange(0, desired_res), np.arange(0, desired_res))) != desired_res: - print("Incorrect image resolution") - sys.exit() - else: - return dn(np.arange(0, desired_res), np.arange(0, desired_res)) - - -def rescale_hmi(image: np.array, orig_res: int, desired_res: int): - hdu_number = 0 - hed = fits.getheader(image[0], hdu_number) - dat = fits.getdata(image[0], ext=0) - 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, dat) - if len(dn(np.arange(0, desired_res), np.arange(0, desired_res))) != desired_res: - print("Incorrect image resolution") - sys.exit() - else: - 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, dat) - if len(dn(np.arange(0, desired_res), np.arange(0, desired_res))) != desired_res: - print("Incorrect image resolution") - sys.exit() - else: - return dn(np.arange(0, desired_res), np.arange(0, desired_res)) - - -# defining data and headers which are used in later steps -hdu_number = 0 - -data = rescale_aia(im171, 1024, 4096) -datb = rescale_aia(im193, 1024, 4096) -datc = rescale_aia(im211, 1024, 4096) -datm = rescale_hmi(imhmi, 1024, 4096) - - -# rescales 'cdelt1' 'cdelt2' 'cpix1' 'cipix2' if 'cdelt1' > 1 -# ensures 'ctype1' 'ctype2' are correctly defined as 'solar_x' and 'solar_y' respectively -# rotates array if 'crota1' is greater than 90 degrees -def filter(aiaa: np.array, aiab: np.array, aiac: np.array, aiam: np.array): - global heda, hedb, hedc, hedm - 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 " - 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, - ) - if hedm["crota1"] > 90: - datm = np.rot90(np.rot90(datm)) - - -filter(im171, im193, im211, imhmi) - - -# removes negative values from an array -def remove_neg(aiaa: np.array, aiab: np.array, aiac: np.array): - global data, datb, datc - data[np.where(data <= 0)] = 0 - datb[np.where(datb <= 0)] = 0 - datc[np.where(datc <= 0)] = 0 - if len(data[data < 0]) != 0: - print("data contains negative") - if len(datb[datb < 0]) != 0: - print("data contains negative") - if len(datc[datc < 0]) != 0: - print("datc contains negative") - - -remove_neg(im171, im193, im211) - -# defines the shape (length) of the array as "s" and the solar radius as "rs" -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"] - - -pix_arc(im171) - -# converts to the Heliographic Stonyhurst coordinate system - - -def to_helio(image: 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 - 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) - - -to_helio(im171) - -# setting up arrays to be used in later processing -# only difference between iarr and bmcool is integer vs. float -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) -# plt.plot(garr) -garr[w] = 1.0 - -# creates sub-arrays of props to isolate column of index 0 and column of index 1 -# what is props?? -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 s - -with np.errstate(divide="ignore"): - t0 = np.log10(datc) - t1 = np.log10(datb) - t2 = np.log10(data) - - -class Bounds: - 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): - 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 - - -threshold(t0) -threshold(t1) -threshold(t2) - - -def set_contour(tval: 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(): - 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 - - -create_mask() - - -def conjunction(): - global bmhot, bmcool, bmmix, cand - cand = bmcool * bmmix * bmhot - - -conjunction() - - -def misid(): - 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 - - -misid() - - -def on_off(): - global circ, cand - circ[:] = 0 - r = (rs / dattoarc) - 10 - inside = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 <= r**2) - circ[inside] = 1.0 - r = (rs / dattoarc) + 40 - outside = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 >= r**2) - circ[outside] = 1.0 - cand = cand * circ - - -on_off() - - -def contours(): - global cand, cont, heir - cand = np.array(cand, dtype=np.uint8) - cont, heir = cv2.findContours(cand, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) - - -contours() - - -def sort(): - global sizes, reord, tmp, cont - 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) - - -sort() - - -# =====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 - 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 - - wh1 = np.where(npix[1] > 0) - wh2 = np.where(npix[1] < 0) - - # =====magnetic cut offs dependant 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 lattitude 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)] - - # ====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======= - - 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 - -# =====sets ident back to max value of iarr====== - -ident = ident - 1 -np.savetxt("ch_summary.txt", props, fmt="%s") - - -from skimage.util import img_as_ubyte - - -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(): - 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 = ) - cs = plt.contour(xgrid, ygrid, slate, colors="white", linewidths=0.5) - plt.savefig("tricolor.png") - plt.close() - - -def plot_mask(slate=slate): - 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") - cs = plt.contour(xgrid, ygrid, slate, colors="black", linewidths=0.5) - cs = plt.contour(xgrid, ygrid, circ, colors="black", linewidths=1.0) - - plt.savefig("CH_mask_" + hedb["DATE"] + ".png", transparent=True) - # plt.close() - - -# ====stores all CH properties in a text file===== - -plot_tricolor() -plot_mask() - -# ====EOF==== diff --git a/CHIMERA_V3.py b/CHIMERA_V3.py deleted file mode 100644 index 48b0111..0000000 --- a/CHIMERA_V3.py +++ /dev/null @@ -1,686 +0,0 @@ -# import required libraries -import glob -import sys - -import astropy.units as u -import cv2 -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 - -plt.style.use(astropy_mpl_style) - -# loading in the images as fits files - -im171 = glob.glob("171.fts") -im193 = glob.glob("193.fts") -im211 = glob.glob("211.fts") -imhmi = glob.glob("hmi.fts") - -# ensure that all images are present - -if im171 == [] or im193 == [] or im211 == [] or imhmi == []: - print("Not all required files present") - sys.exit() - -# Two functions that rescale the aia and hmi images from any original size to any final size - -# didn't normalize by exposure time for hmi because it was equal to 0 - - -def rescale_aia(image: np.array, orig_res: int, desired_res: int): - hed = fits.getheader(image[0], 0) - dat = fits.getdata(image[0], 0) / (hed["EXPTIME"]) - 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, dat) - if len(dn(np.arange(0, desired_res), np.arange(0, desired_res))) != desired_res: - print("Incorrect image resolution") - sys.exit() - else: - 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, dat) - if len(dn(np.arange(0, desired_res), np.arange(0, desired_res))) != desired_res: - print("Incorrect image resolution") - sys.exit() - else: - return dn(np.arange(0, desired_res), np.arange(0, desired_res)) - - -def rescale_hmi(image: np.array, orig_res: int, desired_res: int): - hdu_number = 0 - hed = fits.getheader(image[0], hdu_number) - dat = fits.getdata(image[0], ext=0) - 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, dat) - if len(dn(np.arange(0, desired_res), np.arange(0, desired_res))) != desired_res: - print("Incorrect image resolution") - sys.exit() - else: - 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, dat) - if len(dn(np.arange(0, desired_res), np.arange(0, desired_res))) != desired_res: - print("Incorrect image resolution") - sys.exit() - else: - return dn(np.arange(0, desired_res), np.arange(0, desired_res)) - - -# defining data and headers which are used in later steps -hdu_number = 0 - -data = rescale_aia(im171, 1024, 4096) -datb = rescale_aia(im193, 1024, 4096) -datc = rescale_aia(im211, 1024, 4096) -datm = rescale_hmi(imhmi, 1024, 4096) - - -# rescales 'cdelt1' 'cdelt2' 'cpix1' 'cipix2' if 'cdelt1' > 1 -# ensures 'ctype1' 'ctype2' are correctly defined as 'solar_x' and 'solar_y' respectively -# rotates array if 'crota1' is greater than 90 degrees -def filter(aiaa: np.array, aiab: np.array, aiac: np.array, aiam: np.array): - global heda, hedb, hedc, hedm - 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 " - 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, - ) - if hedm["crota1"] > 90: - datm = np.rot90(np.rot90(datm)) - - -filter(im171, im193, im211, imhmi) - - -# removes negative values from an array -def remove_neg(aiaa: np.array, aiab: np.array, aiac: np.array): - global data, datb, datc - data[np.where(data <= 0)] = 0 - datb[np.where(datb <= 0)] = 0 - datc[np.where(datc <= 0)] = 0 - if len(data[data < 0]) != 0: - print("data contains negative") - if len(datb[datb < 0]) != 0: - print("data contains negative") - if len(datc[datc < 0]) != 0: - print("datc contains negative") - - -remove_neg(im171, im193, im211) - -# defines the shape (length) of the array as "s" and the solar radius as "rs" -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"] - - -pix_arc(im171) - -# converts to the Heliographic Stonyhurst coordinate system - - -def to_helio(image: 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 - 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) - - -to_helio(im171) - -# setting up arrays to be used in later processing -# only difference between iarr and bmcool is integer vs. float -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) -# plt.plot(garr) -garr[w] = 1.0 - -# creates sub-arrays of props to isolate column of index 0 and column of index 1 -# what is props?? -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 s - -with np.errstate(divide="ignore"): - t0 = np.log10(datc) - t1 = np.log10(datb) - t2 = np.log10(data) - - -class Bounds: - 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): - 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 - - -threshold(t0) -threshold(t1) -threshold(t2) - - -def set_contour(tval: 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(): - 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 - - -create_mask() - - -def conjunction(): - global bmhot, bmcool, bmmix, cand - cand = bmcool * bmmix * bmhot - - -conjunction() - - -def misid(): - 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 - - -misid() - - -def on_off(): - global circ, cand - circ[:] = 0 - r = (rs / dattoarc) - 10 - inside = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 <= r**2) - circ[inside] = 1.0 - r = (rs / dattoarc) + 40 - outside = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 >= r**2) - circ[outside] = 1.0 - cand = cand * circ - - -on_off() - - -def contours(): - global cand, cont, heir - cand = np.array(cand, dtype=np.uint8) - cont, heir = cv2.findContours(cand, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) - - -contours() - - -def sort(): - global sizes, reord, tmp, cont - 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) - - -sort() - - -# =====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 - 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 - - wh1 = np.where(npix[1] > 0) - wh2 = np.where(npix[1] < 0) - - # =====magnetic cut offs dependant 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 lattitude 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)] - - # ====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======= - - 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 - -# =====sets ident back to max value of iarr====== - -ident = ident - 1 -np.savetxt("ch_summary.txt", props, fmt="%s") - - -from skimage.util import img_as_ubyte - - -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(): - 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 = ) - cs = plt.contour(xgrid, ygrid, slate, colors="white", linewidths=0.5) - plt.savefig("tricolor.png") - plt.close() - - -def plot_mask(slate=slate): - 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") - cs = plt.contour(xgrid, ygrid, slate, colors="black", linewidths=0.5) - cs = plt.contour(xgrid, ygrid, circ, colors="black", linewidths=1.0) - - plt.savefig("CH_mask_" + hedb["DATE"] + ".png", transparent=True) - # plt.close() - - -# ====stores all CH properties in a text file===== - -plot_tricolor() -plot_mask() - -# ====EOF==== diff --git a/chimerapy/chimera.py b/chimerapy/chimera.py index d5a80fe..098bb9e 100644 --- a/chimerapy/chimera.py +++ b/chimerapy/chimera.py @@ -209,7 +209,7 @@ def to_helio(image: np.array): hg = hpc.transform_to(sunpy.coordinates.frames.HeliographicStonyhurst) global csys csys = wcs.WCS(hedb) - + to_helio(im171) diff --git a/chimerapy/chimera_copy.py b/chimerapy/chimera_copy.py index 68b0966..32ee6d2 100644 --- a/chimerapy/chimera_copy.py +++ b/chimerapy/chimera_copy.py @@ -1,50 +1,34 @@ """Package for Coronal Hole Identification Algorithm""" -import glob -import sys import astropy.units as u import cv2 import mahotas import matplotlib.pyplot as plt import numpy as np -import scipy -from sunpy.map import Map -from sunpy.coordinates import frames -from sunpy.coordinates import propagate_with_solar_surface -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 -from astropy.wcs import WCS -from sunpy.map import all_coordinates_from_map -from sunpy.coordinates import HeliographicStonyhurst - - -#--noverify - -plt.style.use(astropy_mpl_style) - -"""Defining the paths for example files used to run program locally""" - -file_path = "./" - -INPUT_FILES = {"aia171": 'http://jsoc.stanford.edu/data/aia/synoptic/2016/09/22/H1000/AIA20160922_1030_0171.fits', -"aia193": 'http://jsoc.stanford.edu/data/aia/synoptic/2016/09/22/H1000/AIA20160922_1030_0193.fits', -"aia211": 'http://jsoc.stanford.edu/data/aia/synoptic/2016/09/22/H1000/AIA20160922_1030_0211.fits', -"hmi_mag": 'http://jsoc.stanford.edu/data/hmi/fits/2016/09/22/hmi.M_720s.20160922_010000_TAI.fits', +from sunpy.coordinates import (HeliographicStonyhurst, + propagate_with_solar_surface) +from sunpy.map import Map, all_coordinates_from_map + +INPUT_FILES = { + "aia171": "http://jsoc.stanford.edu/data/aia/synoptic/2016/09/22/H1000/AIA20160922_1030_0171.fits", + "aia193": "http://jsoc.stanford.edu/data/aia/synoptic/2016/09/22/H1000/AIA20160922_1030_0193.fits", + "aia211": "http://jsoc.stanford.edu/data/aia/synoptic/2016/09/22/H1000/AIA20160922_1030_0211.fits", + "hmi_mag": "http://jsoc.stanford.edu/data/hmi/fits/2016/09/22/hmi.M_720s.20160922_010000_TAI.fits", } -im171 = Map(INPUT_FILES['aia171']) -im193 = Map(INPUT_FILES['aia193']) -im211 = Map(INPUT_FILES['aia211']) -imhmi = Map(INPUT_FILES['hmi_mag']) +im171 = Map(INPUT_FILES["aia171"]) +im193 = Map(INPUT_FILES["aia193"]) +im211 = Map(INPUT_FILES["aia211"]) +imhmi = Map(INPUT_FILES["hmi_mag"]) -def rescale(proj_to: sunpy.map.Map, input_map: sunpy.map.Map): +def reproject_diff_rot(target_wcs: wcs.wcs.WCS, input_map: sunpy.map.Map): """ Rescale the input aia image dimensions. @@ -56,26 +40,27 @@ def rescale(proj_to: sunpy.map.Map, input_map: sunpy.map.Map): Returns ------- array: 'np.array' - - """ + + """ with propagate_with_solar_surface(): - map1 = proj_to.reproject_to(input_map.wcs) - new_x_scale = map1.scale[0].to(u.arcsec / u.pixel).value - new_y_scale = map1.scale[1].to(u.arcsec / u.pixel).value - map1.meta['cdelt1'] = new_x_scale - map1.meta['cdelt2'] = new_y_scale - map1.meta['cunit1'] = 'arcsec' - map1.meta['cunit2'] = 'arcsec' - return map1 - -im171 = rescale(im171, im171) -im193 = rescale(im171, im193) -im211 = rescale(im171, im211) -imhmi = rescale(im171, imhmi) + amap = input_map.reproject_to(target_wcs.wcs) + new_x_scale = amap.scale[0].to(u.arcsec / u.pixel).value + new_y_scale = amap.scale[1].to(u.arcsec / u.pixel).value + amap.meta["cdelt1"] = new_x_scale + amap.meta["cdelt2"] = new_y_scale + amap.meta["cunit1"] = "arcsec" + amap.meta["cunit2"] = "arcsec" + return amap + + +im193 = reproject_diff_rot(im171, im193) +im211 = reproject_diff_rot(im171, im211) +imhmi = reproject_diff_rot(im171, imhmi) + def filter(map1: np.array, map2: np.array, map3: np.array): """ - Defines headers and filters aia arrays to meet header requirements. Removes 0 values from each array. + Removes negative values from each map by setting each equal to zero Parameters ---------- @@ -96,59 +81,83 @@ def filter(map1: np.array, map2: np.array, map3: np.array): return map1, map2, map3 + im171, im193, im211 = filter(im171, im193, im211) -"""defines the shape of the arrays as "s" and "rs" as the solar radius""" -s = np.shape(im171.data) -#do we want the solar radius in arcsec or pixels? -rs = im171.rsun_obs -rs_pixels = im171.rsun_obs/im171.scale[0] +def shape(map1: sunpy.map.Map, map2: sunpy.map.Map, map3: sunpy.map.Map): + """ + defines the shape of the arrays as "s" and "rs" as the solar radius + + Parameters + ---------- + map1: 'sunpy.map.Map' + map2: 'sunpy.map.Map' + map3: 'sunpy.map.Map' + Returns + ------- + s: 'tuple' + rs: 'astropy.units.quantity.Quantity' + rs_pixels: 'astropy.units.quantity.Quantity' + + """ + + im171, im193, im211 = filter(map1, map2, map3) + # defines the shape of the arrays as "s" and "rs" as the solar radius + s = np.shape(im171.data) + rs = im171.rsun_obs + print(rs) + rs_pixels = im171.rsun_obs / im171.scale[0] + return s, rs, rs_pixels +s, rs, rs_pixels = shape(im171, im193, im211) -def pix_arc(map: sunpy.map.Map): - ''' - Calculates the conversion value of pixel to arcsec + +def pix_arc(amap: sunpy.map.Map): + """ + Defines conversion values between pixels and arcsec Parameters ---------- - map1: 'sunpy.map.Map' - map2: 'sunpy.map.Map' + amap: 'sunpy.map.Map' Returns - ''' - dattoarc = map.scale[0].value - conver = ((s[0]) / 2) * dattoarc / map.meta["cdelt1"] - (s[1] / 2) - convermul = dattoarc / map.meta["cdelt1"] + """ + dattoarc = amap.scale[0].value + s = amap.dimensions + conver = (s.x/2) * amap.scale[0].value / amap.meta['cdelt1'], (s.y / 2) + convermul = dattoarc / amap.meta["cdelt1"] return dattoarc, conver, convermul dattoarc, conver, convermul = pix_arc(im171) +print(conver) -def to_helio(map: sunpy.map.Map): + +def to_helio(amap: sunpy.map.Map): """ - Converts arrays to the Heliographic Stonyhurst coordinate system + Converts maps to the Heliographic Stonyhurst coordinate system Parameters ---------- - map: 'sunpy.map.Map' + amap: 'sunpy.map.Map' Returns ------- - x, y: 'astropy.units.quantity.Quantity' hpc: 'astropy.coordinates.sky_coordinate.SkyCoord' hg: 'astropy.coordinates.sky_coordinate.SkyCoord' csys: 'astropy.wcs.wcs.WCS' """ - hpc = all_coordinates_from_map(map) + hpc = all_coordinates_from_map(amap) hg = hpc.transform_to(sunpy.coordinates.frames.HeliographicStonyhurst) - csys = wcs.WCS(dict(map.meta)) + csys = wcs.WCS(dict(amap.meta)) return hpc, hg, csys + hpc, hg, csys = to_helio(im171) """Setting up arrays to be used in later processing""" @@ -165,8 +174,14 @@ def to_helio(map: sunpy.map.Map): 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:1024, 0:1024] -pix_size = (2000 * u.arcsec).value -garr = Gaussian2D(1, im171.reference_pixel.x.value, im171.reference_pixel.y.value, pix_size / im171.scale[0].value, pix_size / im171.scale[0].value)(x, y) +width = (2000 * u.arcsec) +garr = Gaussian2D( + 1, + im171.reference_pixel.x.to_value(u.pix), + im171.reference_pixel.y.to_value(u.pix), + width / im171.scale[0], + width / im171.scale[1], +)(x, y) garr[w] = 1.0 """creates sub-arrays of props to isolate column of index 0 and column of index 1""" @@ -229,8 +244,10 @@ def to_helio(map: sunpy.map.Map): ) """define threshold values in log space""" + + def log_dat(map1: sunpy.map.Map, map2: sunpy.map.Map, map3: sunpy.map.Map): - ''' + """ Takes the log base-10 of all sunpy map data Parameters @@ -244,13 +261,14 @@ def log_dat(map1: sunpy.map.Map, map2: sunpy.map.Map, map3: sunpy.map.Map): t0: 'np.array' t1: 'np.array' t2: 'np.array' - ''' + """ with np.errstate(divide="ignore"): t0 = np.log10(map1.data) t1 = np.log10(map2.data) t2 = np.log10(map3.data) return t0, t1, t2 + t0, t1, t2 = log_dat(im171, im193, im211) @@ -276,7 +294,8 @@ def new_s(self, new_slope): t1b = Bounds(1.4, 3.0, 255) t2b = Bounds(1.2, 3.9, 255) -#set to also take in boundaries + +# set to also take in boundaries def set_contour(t0: np.array, t1: np.array, t2: np.array): """ Threshold arrays based on desired boundaries and sets contours. @@ -295,15 +314,15 @@ def set_contour(t0: np.array, t1: np.array, t2: np.array): """ if t0 is not None and t1 is not None and t2 is not None: - #set the threshold and contours for t0 + # set the threshold and contours for t0 t0[np.where(t0 < t0b.upper)] = t0b.upper t0[np.where(t0 > t0b.lower)] = t0b.lower t0 = np.array(((t0 - t0b.upper) / (t0b.lower - t0b.upper)) * t0b.slope, dtype=np.float32) - #set the threshold and contours for t1 + # set the threshold and contours for t1 t1[np.where(t1 < t1b.upper)] = t1b.upper t1[np.where(t1 > t1b.lower)] = t2b.lower t1 = np.array(((t1 - t1b.upper) / (t1b.lower - t1b.upper)) * t1b.slope, dtype=np.float32) - #set the threshold and contours for t2 + # set the threshold and contours for t2 t2[np.where(t2 < t2b.upper)] = t2b.upper t2[np.where(t2 > t2b.lower)] = t2b.lower t2 = np.array(((t2 - t2b.upper) / (t2b.lower - t2b.upper)) * t2b.slope, dtype=np.float32) @@ -314,7 +333,10 @@ def set_contour(t0: np.array, t1: np.array, t2: np.array): t0, t1, t2 = set_contour(t0, t1, t2) -def create_mask(tm1: np.array, tm2: np.array, tm3: np.array, map1: sunpy.map.Map, map2: sunpy.map.Map, map3: sunpy.map.Map): + +def create_mask( + tm1: np.array, tm2: np.array, tm3: np.array, map1: sunpy.map.Map, map2: sunpy.map.Map, map3: sunpy.map.Map +): """ Creates 3 segmented bitmasks @@ -326,7 +348,7 @@ def create_mask(tm1: np.array, tm2: np.array, tm3: np.array, map1: sunpy.map.Map map1: 'sunpy.map.Map' map2: 'sunpy.map.Map' map3: 'sunpy.map.Map' - + Returns ------- bmmix: 'np.array' @@ -343,9 +365,10 @@ def create_mask(tm1: np.array, tm2: np.array, tm3: np.array, map1: sunpy.map.Map bmmix, bmhot, bmcool = create_mask(t0, t1, t2, im171, im193, im211) -#conjunction of 3 bitmasks +# conjunction of 3 bitmasks cand = bmcool * bmmix * bmhot + def misid(can: np.array, cir: np.array, xgir: np.array, ygir: np.array, thresh_rad: int): """ Removes off-detector mis-identification @@ -362,7 +385,7 @@ def misid(can: np.array, cir: np.array, xgir: np.array, ygir: np.array, thresh_r 'np.array' """ - #make r a function argument, give name and unit + # make r a function argument, give name and unit r = thresh_rad w = np.where((xgir - center[0]) ** 2 + (ygir - center[1]) ** 2 <= thresh_rad**2) cir[w] = 1.0 @@ -395,17 +418,16 @@ def on_off(cir: np.array, can: np.array): outside = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 >= r**2) cir[outside] = 1.0 can = can * cir - plt.figure() - plt.imshow(cand, cmap='viridis') - plt.show return can + cand = on_off(circ, cand) + def contour_data(cand: np.array): """ Contours the identified datapoints - + Parameters ---------- cand: 'np.array' @@ -419,15 +441,16 @@ def contour_data(cand: np.array): """ cand = np.array(cand, dtype=np.uint8) cont, heir = cv2.findContours(cand, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) - #I think cont might be the x-y coordinates in pixels? return cand, cont, heir + cand, cont, heir = contour_data(cand) + def sort(cont: tuple): """ Sorts the contours by size - + Parameters ---------- cont: 'tuple' @@ -457,14 +480,16 @@ def sort(cont: tuple): # =====cycles through contours========= -def extent(map: sunpy.map.Map, cont: tuple): +def extent(amap: sunpy.map.Map, cont: tuple, xpos: int, ypos: int): """ Finds coronal hole extent in latitude and longitude Parameters ---------- - map: 'sunpy.map.Map' + amap: 'sunpy.map.Map' cont: 'tuple' + xpos: 'int' + ypos: 'int' Returns ------- @@ -475,8 +500,8 @@ def extent(map: sunpy.map.Map, cont: tuple): """ - coord_hpc = map.world2pix(cont) - maxlat = coord_hpc.transform_to(HeliographicStonyhurst).lat.max() + coord_hpc = amap.world2pix(cont) + maxlat = coord_hpc.transform_to(HeliographicStonyhurst).lat.max() maxlon = coord_hpc.transform_to(HeliographicStonyhurst).lon.max() minlat = coord_hpc.transform_to(HeliographicStonyhurst).lat.min() minlon = coord_hpc.transform_to(HeliographicStonyhurst).lat.min() @@ -488,10 +513,9 @@ def extent(map: sunpy.map.Map, cont: tuple): return maxlat, maxlon, minlat, minlon, centlat, centlon - def coords(i, csys, cont): """ - Finds coordinates of CH boundaries + Finds coordinates of CH boundaries in world coordinates Parameters ---------- @@ -684,7 +708,7 @@ def ins_prop( slate[:] = 0 else: - """Create a simple centre point if coronal hole regions is not quiet""" + """Create a simple centre point if coronal hole region is not quiet""" arccent = csys.all_pix2world(cent[0], cent[1], 0) @@ -863,7 +887,7 @@ def plot_tricolor(): """ - tricolorarray = np.zeros((4096, 4096, 3)) + tricolorarray = np.zeros((1024, 1024, 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)) @@ -921,4 +945,13 @@ def plot_mask(slate=slate): plot_mask() if __name__ == "__main__": - import_functions(INPUT_FILES['aia171'], INPUT_FILES['aia193'], INPUT_FILES['aia211'], INPUT_FILES['hmi_mag']) \ No newline at end of file + import_functions( + INPUT_FILES["aia171"], INPUT_FILES["aia193"], INPUT_FILES["aia211"], INPUT_FILES["hmi_mag"] + ) + +''' +Document detailing process summary and all functions/variables: + +https://docs.google.com/document/d/1V5LkZq_AAHdTrGsnCl2hoYhm_fvyjfuzODiEHbt4ebo/edit?usp=sharing + +''' diff --git a/chimerapy/tests/test_chimera.py b/chimerapy/tests/test_chimera.py index c3da28b..714c112 100644 --- a/chimerapy/tests/test_chimera.py +++ b/chimerapy/tests/test_chimera.py @@ -1,4 +1,3 @@ -import glob import os import warnings @@ -7,64 +6,24 @@ from astropy.io import fits from astropy.modeling.models import Gaussian2D -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, -) from chimerapy.chimera import * - -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", +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", + "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 img_present(): assert im171 != [] or im193 != [] or im211 != [] or imhmi != [], "Not all required files present" @@ -365,7 +324,7 @@ def test_fill_polygon(): 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]]]]) + cont = np.array([[[[1, 2], [1, 3], [2, 3], [2, 2]]]]) slate = np.array(iarr) @@ -384,14 +343,12 @@ def test_fill_polygon(): ] > 0 ): - mahotas.polygon.fill_polygon(polygon_vertices, slate) - print("After filling polygon:") print(slate) iarr[np.where(slate == 1)] = 0 - slate[:] = 0 + slate[:] = 0 assert np.array_equal( slate, @@ -421,8 +378,8 @@ 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]), + 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)