diff --git a/chimerapy/chimera_copy.py b/chimerapy/chimera_copy.py index 7232d68..1e4bc71 100644 --- a/chimerapy/chimera_copy.py +++ b/chimerapy/chimera_copy.py @@ -1,957 +1,160 @@ -"""Package for Coronal Hole Identification Algorithm""" +import copy + -import astropy.units as u -import cv2 -import mahotas -import matplotlib.pyplot as plt import numpy as np +from matplotlib import pyplot as plt +from matplotlib.colors import LogNorm import sunpy import sunpy.map -from astropy import wcs -from astropy.modeling.models import Gaussian2D -from astropy.wcs import WCS -from skimage.util import img_as_ubyte -from sunpy.coordinates import (HeliographicStonyhurst, - propagate_with_solar_surface) -from sunpy.map import Map, all_coordinates_from_map -from sunpy.coordinates import frames - - -INPUT_FILES = { - "aia171": "http://jsoc.stanford.edu/data/aia/synoptic/2024/01/31/H1000/AIA20240131_1000_0171.fits", - "aia193": "http://jsoc.stanford.edu/data/aia/synoptic/2024/01/31/H1000/AIA20240131_1000_0193.fits", - "aia211": "http://jsoc.stanford.edu/data/aia/synoptic/2024/01/31/H1000/AIA20240131_1000_0211.fits", - "hmi_mag": "http://jsoc.stanford.edu/data/hmi/fits/2024/01/31/hmi.M_720s.20240131_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.plot -im193.plot -im211.plot -imhmi.plot - -def reproject_diff_rot(target_wcs: wcs.wcs.WCS, input_map: sunpy.map.Map): - """ - Rescale the input aia image dimensions. - - Parameters - ---------- - proj_to: 'sunpy.map.Map' - input_map: 'sunpy.map.Map - - Returns - ------- - array: 'np.array' - - """ - with frames.Helioprojective.assume_spherical_screen(target_wcs.observer_coordinate): - with propagate_with_solar_surface(): - 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): - """ - Removes negative values from each map by setting each equal to zero - - Parameters - ---------- - map1: 'sunpy.map.Map' - map2: 'sunpy.map.Map' - map3: 'sunpy.map.Map' - - Returns - ------- - map1: 'sunpy.map.Map' - map2: 'sunpy.map.Map' - map3: 'sunpy.map.Map' - - """ - map1.data[np.where(map1.data <= 0)] = 0 - map2.data[np.where(map2.data <= 0)] = 0 - map3.data[np.where(map3.data <= 0)] = 0 - - return map1, map2, map3 - - -im171, im193, im211 = filter(im171, im193, im211) - - - - -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(amap: sunpy.map.Map): - """ - Defines conversion values between pixels and arcsec - - Parameters - ---------- - amap: 'sunpy.map.Map' - - Returns - - """ - 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(amap: sunpy.map.Map): - """ - Converts maps to the Heliographic Stonyhurst coordinate system - - Parameters - ---------- - amap: 'sunpy.map.Map' - - Returns - ------- - hpc: 'astropy.coordinates.sky_coordinate.SkyCoord' - hg: 'astropy.coordinates.sky_coordinate.SkyCoord' - csys: 'astropy.wcs.wcs.WCS' - - - """ - hpc = all_coordinates_from_map(amap) - hg = hpc.transform_to(sunpy.coordinates.frames.HeliographicStonyhurst) - # Filter the header to contain only ASCII characters and exclude specified keys - filtered_header = {key: value for key, value in amap.meta.items() if key not in ['keycomments', 'comment']} - csys = wcs.WCS(dict(filtered_header)) - return hpc, hg, csys - - -hpc, hg, csys = 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:1024, 0:1024] -width = 2000 * u.arcsec -garr = Gaussian2D( - 1, - im171.reference_pixel.x.to_value(u.pix), - im171.reference_pixel.y.to_value(u.pix), - width.value / im171.scale[0].value, - width.value / im171.scale[1].value, -)(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""" - - -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 - ---------- - map1: 'sunpy.map.Map' - map2: 'sunpy.map.Map' - map3: 'sunpy.map.Map' - - Returns - ------- - 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) - - -class Bounds: - """Class 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) - - -# 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. - - Parameters - ---------- - t0: 'np.array' - t1: 'np.array' - t2: 'np.array'' - - Returns - ------- - 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 - 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 - 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 - 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) - else: - print("Must input valid logarithmic arrays") - return t0, t1, t2 - - -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 -): - """ - Creates 3 segmented bitmasks - - Parameters - ------- - tm1: 'np.array' - tm2: 'np.array' - tm3: 'np.array' - map1: 'sunpy.map.Map' - map2: 'sunpy.map.Map' - map3: 'sunpy.map.Map' - - Returns - ------- - bmmix: 'np.array' - bmhot: 'np.array' - bmcool: 'np.array' - - """ - with np.errstate(divide="ignore", invalid="ignore"): - bmmix[np.where(tm3 / tm1 >= ((np.mean(map1.data) * 0.6357) / (np.mean(map3.data))))] = 1 - bmhot[np.where(tm1 + tm2 < (0.7 * (np.mean(map2.data) + np.mean(map3.data))))] = 1 - bmcool[np.where(tm3 / tm2 >= ((np.mean(map2.data) * 1.5102) / (np.mean(map2.data))))] = 1 - return bmmix, bmhot, bmcool - - -bmmix, bmhot, bmcool = create_mask(t0, t1, t2, im171, im193, im211) - -# 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 - - Parameters - ---------- - can: 'np.array' - cir: 'np.array' - xgir: 'np.array' - ygir: 'np.array' - - Returns - ------- - 'np.array' - - """ - # 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 - cand = can * cir - return r, w, cir, cand - - -r, w, cir_off, cand_off = misid(cand, circ, xgrid, ygrid, (s[1] / 2.0) - 100) - - -def on_off(cir: np.array, can: np.array): - """ - Seperates on-disk and off-limb coronal holes - - Parameters - ---------- - cir: 'np.array' - can: 'np.array' - - Returns - ------- - 'np.array' - - """ - cir[:] = 0 - r = (rs.value / dattoarc) - 10 - inside = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 <= r**2) - cir[inside] = 1.0 - r = (rs.value / dattoarc) + 40 - outside = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 >= r**2) - cir[outside] = 1.0 - can = can * cir - return can - - -#cand = on_off(circ, cand) - - -def contour_data(cand: np.array): - """ - Contours the identified datapoints - - Parameters - ---------- - cand: 'np.array' - - Returns - ------- - cand: 'np.array' - cont: 'tuple' - heir: 'np.array' - - """ - cand = np.array(cand, dtype=np.uint8) - cont, heir = cv2.findContours(cand, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) - return cand, cont, heir - - -cand, cont, heir = contour_data(cand) - - -def sort(cont: tuple): - """ - Sorts the contours by size - - Parameters - ---------- - cont: 'tuple' - - Returns - ------- - reord: 'list' - tmp: 'list' - cont: 'list' - sizes: 'list' - - """ - sizes = np.array([]) - 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) - return cont, sizes, reord, tmp - - -cont, sizes, reord, tmp = sort(cont) - - -# =====cycles through contours========= - - -def extent(amap: sunpy.map.Map, cont: tuple, xpos: int, ypos: int): - """ - Finds coronal hole extent in latitude and longitude - - Parameters - ---------- - amap: 'sunpy.map.Map' - cont: 'tuple' - xpos: 'int' - ypos: 'int' - - 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' - - """ +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) - 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() - - # =====CH centroid in lat/lon======= + # 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 - centlat = hg.lat[int(ypos), int(xpos)] - centlon = hg.lon[int(ypos), int(xpos)] - return maxlat, maxlon, minlat, minlon, centlat, centlon + 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), + ) -def coords(i, csys, cont): - """ - Finds coordinates of CH boundaries in world coordinates + # 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) - Parameters - ---------- - i: 'int' - csys: 'astropy.wcs.wcs.WCS' - cont: 'list' + #Finding the corresponding minimum y-index + min_y_index = np.min(non_zero_cool[1]) - 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' - """ + #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") - #exctracting coordinates for contours - contour = cont[i] - x_coords = contour[:, 0, 0] - y_coords = contour[:, 0, 1] - - #finding max and min of x coordinates - max_x = np.argmax(x_coords) - min_x = np.argmin(x_coords) - - #finding manx and min of y coordinates - max_y = np.argmax(y_coords) - min_y = np.argmin(y_coords) - - #Pixel coordinates to world coordinates (in arcsec) - Ywb, Xwb = (csys.pixel_to_world(x_coords[max_x], y_coords[max_x])).Tx, (csys.pixel_to_world(x_coords[max_x], y_coords[max_x])).Ty - Yeb, Xeb = (csys.pixel_to_world(x_coords[min_x], y_coords[min_x])).Tx, (csys.pixel_to_world(x_coords[min_x], y_coords[min_x])).Ty - Ynb, Xnb = (csys.pixel_to_world(y_coords[max_y], x_coords[max_y])).Tx, (csys.pixel_to_world(y_coords[max_y], x_coords[max_y])).Ty - Ysb, Xsb = (csys.pixel_to_world(y_coords[min_y], x_coords[min_y])).Tx, (csys.pixel_to_world(y_coords[min_y], x_coords[min_y])).Ty - - return Ywb, Xwb, Yeb, Xeb, Ynb, Xnb, Ysb, Xsb - - -def ins_prop( - imhmi, - 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 - ---------- - imhmi: '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(imhmi.data[pos[:, 0], pos[:, 1]][np.where(imhmi.data[pos[:, 0], pos[:, 1]] > 0)]) - props[21, ident + 1] = f"{tbpos:.1e}" - tbneg = np.sum(imhmi.data[pos[:, 0], pos[:, 1]][np.where(imhmi.data[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])) - #exctracting coordinates for contours - contour = cont[i] - lengths = [] - # Iterate through each element in cont and calculate its length - for elem in cont: - lengths.append(len(elem)) - x_coords = contour[:, 0, 0] - y_coords = contour[:, 0, 1] - max_x = np.max(x_coords) - max_y = np.max(y_coords) - """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)) + # 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, ) - arcar = area * (dattoarc**2) - if arcar > 1000: - """finds centroid""" - - chpts = len(cont[i]) - cent = [np.mean(x_coords), np.mean(y_coords)] - - """remove quiet sun regions encompassed by coronal holes""" - if (cand[max_x + 1, contour[np.where(x_coords == max_x)[0][0], 0, 1]] > 0) and (iarr[max_x + 1, cont[i][np.where(x_coords == max_x)[0][0], 0, 1]] > 0): - mahotas.polygon.fill_polygon(np.array(list(zip(y_coords, x_coords))), slate) - iarr[np.where(slate == 1)] = 0 - slate[:] = 0 - - else: - """Create a simple centre point if coronal hole region 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.value**2)) or ( - np.sum(np.array(csys.all_pix2world(x_coords, y_coords, 0)) ** 2) > (rs.value**2) - ): - mahotas.polygon.fill_polygon(np.array(list(zip(y_coords, x_coords))), offarr) - else: - """classifies on disk coronal holes""" - - mahotas.polygon.fill_polygon(np.array(list(zip(y_coords, x_coords))), slate) - poslin = np.where(slate == 1) - print(poslin) - slate[:] = 0 - """create an array for magnetic polarity""" - pos_x = np.array((poslin[0] - (s[0] / 2)) * convermul + (s[1] / 2), dtype=np.uint) - pos_y = np.array((poslin[1] - (s[0] / 2)) * convermul + (s[1] / 2), dtype=np.uint) - pos = np.column_stack((pos_x, pos_y)) - npix = list( - np.histogram( - imhmi.data[pos_x, pos_y], - bins=np.arange( - np.round(np.min(imhmi.data[pos_x, pos_y])) - 0.5, - np.round(np.max(imhmi.data[pos_x, pos_y])) + 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) - - """Filters magnetic cutoff values by 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(imhmi.data[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""" - - maxxlon, minxlon, centlat, centlon = extent(i, ypos, xpos, hg, cont) - - """caluclate the mean magnetic field""" - - mB = np.mean(imhmi.data[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, 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( - imhmi, - 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 + #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) - """ - if cmin or cmax: - arr = np.clip(arr, cmin, cmax) - return (b - a) * ((arr - np.min(arr)) / (np.max(arr) - np.min(arr))) + a + #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] -def plot_tricolor(): - """ - Plots a tricolor mask of image data + 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") - Returns - ------- - plot: 'matplotlib.image.AxesImage' + plt.show() - """ - - tricolorarray = np.zeros((1024, 1024, 3)) - - data_a = img_as_ubyte(rescale01(np.log10(im171.data), cmin=1.2, cmax=3.9)) - data_b = img_as_ubyte(rescale01(np.log10(im193.data), cmin=1.4, cmax=3.0)) - data_c = img_as_ubyte(rescale01(np.log10(imhmi.data), 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") - plt.contour(xgrid, ygrid, slate, colors="white", linewidths=0.5) - plt.savefig("tricolor.png") - plt.close() +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 -def plot_mask(slate=slate): - """ - Plots the contour mask +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) - Parameters - ---------- - slate: 'np.array' - Returns - ------- - plot: 'matplotlib.image.AxesImage' +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) - 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) +tri_color_img = make_lupton_rgb(cs_171, cs_193, cs_211, Q=10, stretch=50) +comb_mask = mask211_171 * mask211_193 * mask171_193 - circ[:] = 0 - r = rs / dattoarc - w = np.where((xgrid - center[0]) ** 2 + (ygrid - center[1]) ** 2 <= r.value**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) +fig, axes = plt.subplot_mosaic([["tri", "comb_mask"]], layout="constrained", figsize=(6, 3)) - plt.savefig("CH_mask_" + im193.meta["date-obs"] + ".png", transparent=True) +axes["tri"].imshow(tri_color_img, origin="lower") +axes["comb_mask"].imshow(comb_mask, origin="lower") +plt.show() -plot_tricolor() -plot_mask() +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) -""" -Document detailing process summary and all functions/variables: + 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) -https://docs.google.com/document/d/1V5LkZq_AAHdTrGsnCl2hoYhm_fvyjfuzODiEHbt4ebo/edit?usp=sharing + # 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