Plot SUVI L2 Flare Locations

Purpose:

Use Python to plot the SUVI L2 flare location (flloc) product in various coordinate systems.

Overview

In this example, we will be plotting both the flare location (flloc) and peak location (peak_loc) variables from the flare location product. Flare location is determined by center-of-mass calculations, whereas peak location is determined by the region with the highest irradiance.

We will plot both flloc and peak_loc for each of the following coordinate systems:

• Heliographic Stonyhurst (longitude, latitude)

• Heliographic Carrington (longitude, latitude)

• Pixel (x pixels, y pixels)

• R-Theta (solar radii, degrees)

• Helioprojective Cartesian (x arcsec, y arcsec)

• Heliocentric Cartesian (x AU, y AU, z AU)

• Heliocentric Radial (solar radii, degrees, z solar radii)

Imports

First, we will import the necessary libraries:

__authors__ = "elucas, ajarvis"

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
from datetime import datetime, timedelta
import netCDF4 as nc
import astropy
from astropy.coordinates import SkyCoord
import sunpy
import sunpy.map
import astropy.units as u
from sunpy.coordinates import frames
from sunpy.net import Fido, attrs as a
from astropy.io import fits

Helper Functions

We need to define a few helper functions for use throughout plotting:

# Fetch the SUVI SunPy map to use as a base for plotting.
def create_suvi_sunpymap(date, goes=16, wavelength=131, rng=2):
    ds0 = (date - timedelta(minutes=rng)).strftime("%Y/%m/%d %H:%M:%S")
    ds1 = (date + timedelta(minutes=rng)).strftime("%Y/%m/%d %H:%M:%S")
    q = Fido.search(a.Time(ds0, ds1), a.Instrument.suvi, a.Wavelength(wavelength*u.angstrom))
    tmp_files = Fido.fetch(q)
    # Select only files for level 2 composites from G16
    for tmp_file in tmp_files:
        if (f'g{goes}' in tmp_file) and ('l2' in tmp_file):
            data, header = fits.getdata(tmp_file, ext=1), fits.getheader(tmp_file, ext=1)
            suvi_map = sunpy.map.Map(data, header)
            return suvi_map, data, header
    print('No SUVI map available')
    return None

# Convert time given in the .nc file to datetime.
def convert_time(time_nc):
    date_2000 = datetime(2000, 1, 1, 12, 0)
    date = date_2000 + timedelta(seconds=time_nc)
    return date

# Define legend for this product.
def legend_handles():
    markers = ['o', '*']
    fillstyle = ['none', 'full']
    labels = ['flloc', 'peak_loc']
    f = lambda m, fill: plt.plot([], [], marker=m, color='k', ls='none', fillstyle=fill)[0]
    handles = [f(markers[i], fillstyle[i]) for i in range(len(markers))]
    return handles, labels

Retrieving the SUVI SunPy Map

Let’s use some of the helper functions above to retrieve a SUVI image.

We’ll grab an image of the sun in the 131Å wavelength for datetime 2024-01-23 03:44.

Example data file: dr_suvi-l2-flloc_g16_s20240123T034400Z_e20240123T034800Z_v1-0-6.nc

date = datetime(2024, 1, 23, 3, 44)
goes = 16
wavelength = 131
suvi_flloc_path = f'../../data/suvi/dr_suvi-l2-flloc_g16_s20240123T034400Z_e20240123T034800Z_v1-0-6.nc'
flloc_nc = nc.Dataset(suvi_flloc_path)

Let’s look at the variables contained in the flare location files.

for k in flloc_nc.variables.keys():
    print(k)

time
wavelength
degraded_status
num_flloc
brght_area
srs_status
xrs_status
euv_status
flloc_pix
flloc_hg
flloc_car
flloc_rtheta
flloc_hpc
flloc_hcc
flloc_hcr
tot_flux
peak_flux
peak_loc_pix
peak_loc_hg
peak_loc_car
peak_loc_rtheta
peak_loc_hpc
peak_loc_hcc
peak_loc_hcr

For our examples, we will be looking at the flloc and peak_loc variables. We can see the general structure by inspecting the variables in any coordinate system, for example, in Heliographic Stonyhurst.

print(flloc_nc['flloc_hg'])
print('')
print(flloc_nc['peak_loc_hg'])

<class 'netCDF4._netCDF4.Variable'>
float32 flloc_hg(time, feature_number, wavelength, location)
    long_name: Flare location in Stonyhurst/heliographic coordinates (lon, lat)
    comments: Values provided for on-disk flares only.
    units: degrees, degrees
    _FillValue: -9999.0
unlimited dimensions: time, feature_number, wavelength
current shape = (1, 2, 6, 2)
filling on

<class 'netCDF4._netCDF4.Variable'>
float32 peak_loc_hg(time, feature_number, wavelength, location)
    long_name: Peak bright region flux location in Stonyhurst/heliographic coordinates (lon, lat)
    comments: Values provided for on-disk flares only.
    units: degrees, degrees
    _FillValue: -9999.0
unlimited dimensions: time, feature_number, wavelength
current shape = (1, 2, 6, 2)
filling on

The structure for both in this example is (1, 2, 6, 2):

1: The first dimension of the structure is the time dimension, which will always have 1 value for this product.

2: The second dimension is the number of separate flares within this image.

6: The third dimension is each wavelength, in the order [094, 131, 171, 195, 284, 304].

2: The fourth dimension is the coordinate pair (3 for HCC and HCR), with units as described above.

We need to extract the time (found in the ‘time’ variable), and convert to python datetime in order to fetch a corresponding SUVI SunPy map at the closest time possible.

# Get file time
flloc_time = convert_time(flloc_nc['time'][:][0])
# Access SUVI SunPy map for that time and defined wavelength
suvi_map, suvi_data, suvi_header = create_suvi_sunpymap(flloc_time, goes=goes, wavelength=wavelength)

/usr/share/miniconda/envs/goesr-spwx-examples/lib/python3.12/site-packages/sunpy/net/vso/vso.py:222: SunpyUserWarning: VSO-D400 Bad Request - Invalid wavelength, wavetype of filter specification
  warn_user(resp["error"])

Files Downloaded:   0%|          | 0/8 [00:00<?, ?file/s]





OR_SUVI-L1b-Fe131_G18_s20240230342227_e20240230342237_c20240230342434.fits.gz:   0%|          | 0.00/1.52M [00:00<?, ?B/s]




dr_suvi-l2-ci131_g16_s20240123T034400Z_e20240123T034800Z_v1-0-2.fits:   0%|          | 0.00/1.73M [00:00<?, ?B/s]

OR_SUVI-L1b-Fe131_G16_s20240230345119_e20240230345129_c20240230345302.fits.gz:   0%|          | 0.00/1.48M [00:00<?, ?B/s]



OR_SUVI-L1b-Fe131_G16_s20240230345319_e20240230345319_c20240230345503.fits.gz:   0%|          | 0.00/2.07M [00:00<?, ?B/s]


OR_SUVI-L1b-Fe131_G16_s20240230345219_e20240230345219_c20240230345401.fits.gz:   0%|          | 0.00/1.86M [00:00<?, ?B/s]





OR_SUVI-L1b-Fe131_G18_s20240230342227_e20240230342237_c20240230342434.fits.gz:   0%|          | 1.02k/1.52M [00:00<03:34, 7.11kB/s]




dr_suvi-l2-ci131_g16_s20240123T034400Z_e20240123T034800Z_v1-0-2.fits:   0%|          | 1.02k/1.73M [00:00<04:03, 7.11kB/s]

OR_SUVI-L1b-Fe131_G16_s20240230345119_e20240230345129_c20240230345302.fits.gz:   0%|          | 1.02k/1.48M [00:00<03:28, 7.08kB/s]



OR_SUVI-L1b-Fe131_G16_s20240230345319_e20240230345319_c20240230345503.fits.gz:   0%|          | 1.02k/2.07M [00:00<04:48, 7.19kB/s]


OR_SUVI-L1b-Fe131_G16_s20240230345219_e20240230345219_c20240230345401.fits.gz:   0%|          | 1.02k/1.86M [00:00<04:28, 6.92kB/s]





OR_SUVI-L1b-Fe131_G18_s20240230342227_e20240230342237_c20240230342434.fits.gz:  18%|█▊        | 280k/1.52M [00:00<00:00, 1.38MB/s]




dr_suvi-l2-ci131_g16_s20240123T034400Z_e20240123T034800Z_v1-0-2.fits:  23%|██▎       | 392k/1.73M [00:00<00:00, 1.77MB/s]

OR_SUVI-L1b-Fe131_G16_s20240230345119_e20240230345129_c20240230345302.fits.gz:  22%|██▏       | 328k/1.48M [00:00<00:00, 1.48MB/s]



OR_SUVI-L1b-Fe131_G16_s20240230345319_e20240230345319_c20240230345503.fits.gz:  13%|█▎        | 272k/2.07M [00:00<00:01, 1.22MB/s]


OR_SUVI-L1b-Fe131_G16_s20240230345219_e20240230345219_c20240230345401.fits.gz:  22%|██▏       | 400k/1.86M [00:00<00:00, 1.80MB/s]





OR_SUVI-L1b-Fe131_G18_s20240230342227_e20240230342237_c20240230342434.fits.gz:  91%|█████████ | 1.38M/1.52M [00:00<00:00, 5.35MB/s]






Files Downloaded:  12%|█▎        | 1/8 [00:00<00:03,  1.80file/s]




dr_suvi-l2-ci131_g16_s20240123T034400Z_e20240123T034800Z_v1-0-2.fits:  90%|████████▉ | 1.55M/1.73M [00:00<00:00, 5.53MB/s]

OR_SUVI-L1b-Fe131_G16_s20240230345119_e20240230345129_c20240230345302.fits.gz:  87%|████████▋ | 1.29M/1.48M [00:00<00:00, 4.59MB/s]










OR_SUVI-L1b-Fe131_G16_s20240230345319_e20240230345319_c20240230345503.fits.gz:  69%|██████▊   | 1.42M/2.07M [00:00<00:00, 5.17MB/s]


OR_SUVI-L1b-Fe131_G16_s20240230345219_e20240230345219_c20240230345401.fits.gz:  88%|████████▊ | 1.63M/1.86M [00:00<00:00, 5.84MB/s]












OR_SUVI-L1b-Fe131_G18_s20240230342328_e20240230342328_c20240230342529.fits.gz:   0%|          | 0.00/1.67M [00:00<?, ?B/s]




dr_suvi-l2-ci131_g18_s20240123T034400Z_e20240123T034800Z_v1-0-2.fits:   0%|          | 0.00/1.72M [00:00<?, ?B/s]

OR_SUVI-L1b-Fe131_G18_s20240230342428_e20240230342428_c20240230343023.fits.gz:   0%|          | 0.00/2.05M [00:00<?, ?B/s]





OR_SUVI-L1b-Fe131_G18_s20240230342328_e20240230342328_c20240230342529.fits.gz:   0%|          | 1.02k/1.67M [00:00<03:54, 7.11kB/s]




dr_suvi-l2-ci131_g18_s20240123T034400Z_e20240123T034800Z_v1-0-2.fits:   0%|          | 1.02k/1.72M [00:00<04:06, 6.96kB/s]

OR_SUVI-L1b-Fe131_G18_s20240230342428_e20240230342428_c20240230343023.fits.gz:   0%|          | 1.02k/2.05M [00:00<04:51, 7.04kB/s]




dr_suvi-l2-ci131_g18_s20240123T034400Z_e20240123T034800Z_v1-0-2.fits:  13%|█▎        | 216k/1.72M [00:00<00:01, 983kB/s]

OR_SUVI-L1b-Fe131_G18_s20240230342428_e20240230342428_c20240230343023.fits.gz:  20%|██        | 417k/2.05M [00:00<00:00, 1.88MB/s]





OR_SUVI-L1b-Fe131_G18_s20240230342328_e20240230342328_c20240230342529.fits.gz:  25%|██▌       | 424k/1.67M [00:00<00:00, 1.79MB/s]






Files Downloaded:  75%|███████▌  | 6/8 [00:01<00:00,  6.22file/s]




dr_suvi-l2-ci131_g18_s20240123T034400Z_e20240123T034800Z_v1-0-2.fits:  68%|██████▊   | 1.16M/1.72M [00:00<00:00, 4.25MB/s]






OR_SUVI-L1b-Fe131_G18_s20240230342428_e20240230342428_c20240230343023.fits.gz:  90%|████████▉ | 1.84M/2.05M [00:00<00:00, 6.62MB/s]


Files Downloaded: 100%|██████████| 8/8 [00:01<00:00,  7.10file/s]

Plotting the Flare Locations

Now that the SUVI SunPy map has been retrieved, we can overplot the coordinates from the .nc file onto that image.

Note: Heliographic Stonyhurst and Heliographic Carrington are only defined on-disk, so any features (or pieces of features) that extend off-disk will not be plotted.

Heliographic Stonyhurst Coordinates

# Plot the SUVI SunPy composite image
fig_hg = plt.figure(figsize=(10, 10))
fig_hg.tight_layout(h_pad=-50)
fig_hg.suptitle(f'Heliographic Stonyhurst Coordinates, GOES-{goes}')
ax_hg = fig_hg.add_subplot(projection=suvi_map)
suvi_map.plot(axes=ax_hg, clip_interval=(0.1, 99.9) * u.percent)
num_points = 100

# Ignore INFO message about 'Missing metadata for solar radius: assuming the standard radius of the photosphere.'
# RSUN_REF will be included in future data products.

# Extract flare location from netCDF
flloc_hg = flloc_nc['flloc_hg'][:][0]
peak_hg = flloc_nc['peak_loc_hg'][:][0]

wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)
colors = cm.Reds_r(np.linspace(0.2, 0.6, len(flloc_hg)))

# Flare loc
for flare, c in zip(flloc_hg, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)

    x = flare[wavelength_i][0]
    y = flare[wavelength_i][1]

    # Convert flare location into SkyCoord object in coordinate frame
    flloc_hgs_sc = SkyCoord(x * u.degree, y * u.degree, obstime=flloc_time, observer=suvi_map.observer_coordinate,
                            frame=frames.HeliographicStonyhurst)

    # Plot contours
    stonyhurst_frame = frames.HeliographicStonyhurst(obstime=suvi_map.date)
    constant_lon = SkyCoord(flloc_hgs_sc.lon, np.linspace(-90, 90, num_points) * u.deg,
                            frame=stonyhurst_frame)
    constant_lat = SkyCoord(np.linspace(-90, 90, num_points) * u.deg, flloc_hgs_sc.lat,
                            frame=stonyhurst_frame)

    ax_hg.plot_coord(constant_lon, color=c, alpha=0.7)
    ax_hg.plot_coord(constant_lat, color=c, alpha=0.7)

    # Transform SkyCoord object to the SUVI SunPy map's coordinate frame
    flloc_hgs_sc = flloc_hgs_sc.transform_to(suvi_map.coordinate_frame)

    # Plot flare
    ax_hg.plot_coord(flloc_hgs_sc, color=c, marker='o', fillstyle='none', markersize=8, alpha=0.7)

# Peak loc
for flare, c in zip(peak_hg, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)

    x = flare[wavelength_i][0]
    y = flare[wavelength_i][1]

    # Convert peak location into SkyCoord object in coordinate frame
    flloc_hgs_sc = SkyCoord(x * u.degree, y * u.degree, obstime=flloc_time, observer=suvi_map.observer_coordinate,
                            frame=frames.HeliographicStonyhurst)

    # Transform SkyCoord object to the SUVI SunPy map's coordinate frame
    flloc_hgs_sc = flloc_hgs_sc.transform_to(suvi_map.coordinate_frame)

    # Plot flare
    ax_hg.plot_coord(flloc_hgs_sc, color=c, marker='*', markersize=4, alpha=0.85)

# Add legend for features
handles, labels = legend_handles()
plt.legend(handles, labels, loc=1, framealpha=1)

plt.show()

INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]

Heliographic Carrington Coordinates

# Plot the SUVI SunPy composite image
fig_car = plt.figure(figsize=(10, 10))
fig_car.tight_layout(h_pad=-50)
fig_car.suptitle(f'Heliographic Carrington Coordinates, GOES-{goes}')
ax_car = fig_car.add_subplot(projection=suvi_map)
suvi_map.plot(axes=ax_car, clip_interval=(0.1, 99.9) * u.percent)
num_points = 100

# Extract flare location from netCDF
flloc_car = flloc_nc['flloc_car'][:][0]
peak_car = flloc_nc['peak_loc_car'][:][0]
wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)
colors = cm.YlOrBr_r(np.linspace(0.5, 0.7, len(flloc_car)))

# Flare loc
for flare, c in zip(flloc_car, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)

    x = flare[wavelength_i][0]
    y = flare[wavelength_i][1]

    # Convert flare location into SkyCoord object in coordinate frame
    flloc_car_sc = SkyCoord(x * u.degree, y * u.degree, obstime=flloc_time, observer=suvi_map.observer_coordinate,
                            frame=frames.HeliographicCarrington)

    # Find Carrington rotation for cycle in degrees
    rot_car = sunpy.coordinates.sun.carrington_rotation_number(t=flloc_time)
    rot_mod = rot_car % 1
    rot_deg = rot_mod * 360.

    # Plot contours, offset by Carrington rotation
    carrington_frame = frames.HeliographicCarrington(obstime=suvi_map.date, observer=suvi_map.observer_coordinate)
    constant_lon = SkyCoord(flloc_car_sc.lon, np.linspace(-90, 90, num_points) * u.deg,
                            frame=carrington_frame)
    constant_lat = SkyCoord(np.linspace(-90. - rot_deg, 90. - rot_deg, num_points) * u.deg, flloc_car_sc.lat,
                            frame=carrington_frame)

    ax_car.plot_coord(constant_lon, color=c, alpha=0.7)
    ax_car.plot_coord(constant_lat, color=c, alpha=0.7)

    # Transform SkyCoord object to the SUVI SunPy map's coordinate frame
    flloc_car_sc = flloc_car_sc.transform_to(suvi_map.coordinate_frame)

    # Plot flare
    ax_car.plot_coord(flloc_car_sc, color=c, marker='o', fillstyle='none', markersize=8, alpha=0.7)

# Peak loc
for flare, c in zip(peak_car, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)

    x = flare[wavelength_i][0]
    y = flare[wavelength_i][1]

    # Convert peak location into SkyCoord object in coordinate frame
    flloc_car_sc = SkyCoord(x * u.degree, y * u.degree, obstime=flloc_time, observer=suvi_map.observer_coordinate,
                            frame=frames.HeliographicCarrington)

    # Transform SkyCoord object to the SUVI SunPy map's coordinate frame
    flloc_car_sc = flloc_car_sc.transform_to(suvi_map.coordinate_frame)

    # Plot flare
    ax_car.plot_coord(flloc_car_sc, color=c, marker='*', markersize=4, alpha=0.85)

# Add legend for features
handles, labels = legend_handles()
plt.legend(handles, labels, loc=1, framealpha=1)

plt.show()

Pixel Coordinates

# Plot the SUVI SunPy composite image
fig_px = plt.figure(figsize=(10, 10))
fig_px.tight_layout()
ax_px = fig_px.add_subplot(projection=suvi_map)
suvi_map.plot(axes=ax_px, clip_interval=(0.1, 99.9) * u.percent)
fig_px.suptitle(f'Pixel Coordinates, GOES-{goes}')
num_points = 100

# Extract flare location from netCDF
flloc_pix = flloc_nc['flloc_pix'][:][0]
peak_pix = flloc_nc['peak_loc_pix'][:][0]
wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)
colors = cm.Greens_r(np.linspace(0.3, 0.6, len(flloc_pix)))

# Flare loc
for flare, c in zip(flloc_pix, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)

    x = flare[wavelength_i][0]
    y = flare[wavelength_i][1]

    # Plot crosshairs
    ax_px.axvline(x=x, color=c, alpha=0.7)
    ax_px.axhline(y=y, color=c, alpha=0.7)

    # Plot flare
    ax_px.plot_coord(x, y, color=c, marker='o', fillstyle='none', markersize=8, alpha=0.7)

# Peak loc
for flare, c in zip(peak_pix, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)

    x = flare[wavelength_i][0]
    y = flare[wavelength_i][1]

    # Plot flare
    ax_px.plot_coord(x, y, color=c, marker='*', markersize=4, alpha=0.85)

# Add legend for features
handles, labels = legend_handles()
plt.legend(handles, labels, loc=1, framealpha=1)

plt.show()

R-Theta Coordinates

# Plot the SUVI SunPy map
fig_rt = plt.figure(figsize=(10, 10))
fig_rt.suptitle(f'R-Theta Coordinates, GOES-{goes}')
fig_rt.tight_layout(h_pad=-50)
ax_rt = fig_rt.add_subplot(projection=suvi_map)
suvi_map.plot(axes=ax_rt, clip_interval=(0.1, 99.9) * u.percent)
num_points = 100

# Extract flare location from netCDF
flloc_rt = flloc_nc['flloc_rtheta'][:][0]
peak_rt = flloc_nc['peak_loc_rtheta'][:][0]
wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)
colors = cm.Purples_r(np.linspace(0.3, 0.6, len(flloc_rt)))

# Flare loc
for flare, c in zip(flloc_rt, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)
    # Need to rotate psi by 90 degrees, according to SkyCoord cylindrical specifications
    rho = flare[wavelength_i][0]
    psi = flare[wavelength_i][1] + 90.

    # R-Theta trigonometry to find correct z for initializing SkyCoord object
    # z is assumed to be 0 if feature is off-disk
    solar_radius = 1.
    if -1. <= rho <= 1.:
        angle = np.arcsin(rho / solar_radius)
        z = solar_radius * np.cos(angle)
    else:
        z = 0.

    # Create SkyCoord object in heliocentric frame
    flloc_rt_sc = SkyCoord(rho=rho * u.solRad, psi=psi * u.deg, z=z * u.solRad, representation_type='cylindrical',
                           obstime=flloc_time, observer=suvi_map.observer_coordinate, frame=frames.Heliocentric)

    # Plot r/psi 'crosshairs'
    heliocentric_frame = frames.Heliocentric(obstime=suvi_map.date, observer=suvi_map.observer_coordinate)
    r_vector = SkyCoord(rho=flloc_rt_sc.rho, psi=np.linspace(-360, 0, num_points) * u.deg, z=z * u.solRad,
                        frame=heliocentric_frame, representation_type='cylindrical')
    constant_psi = SkyCoord(rho=np.linspace(0, rho, num_points) * u.solRad, psi=flloc_rt_sc.psi, z=z * u.solRad,
                            frame=heliocentric_frame, representation_type='cylindrical')

    ax_rt.plot_coord(r_vector, color=c, alpha=0.7)
    ax_rt.plot_coord(constant_psi, color=c, alpha=0.7)

    # Plot flare
    flloc_rt_sc = flloc_rt_sc.transform_to(suvi_map.coordinate_frame)
    ax_rt.plot_coord(flloc_rt_sc, color=c, marker='o', fillstyle='none', markersize=8, alpha=0.7)

# Peak loc
for flare, c in zip(peak_rt, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)
    # Need to rotate psi by 90 degrees, according to SkyCoord cylindrical specifications
    rho = flare[wavelength_i][0]
    psi = flare[wavelength_i][1] + 90.

    # R-Theta trigonometry to find correct z for initializing SkyCoord object
    # z is assumed to be 0 if feature is off-disk
    solar_radius = 1.
    if -1. <= rho <= 1.:
        angle = np.arcsin(rho / solar_radius)
        z = solar_radius * np.cos(angle)
    else:
        z = 0.

    # Create SkyCoord object in heliocentric frame
    flloc_rt_sc = SkyCoord(rho=rho * u.solRad, psi=psi * u.deg, z=z * u.solRad, representation_type='cylindrical',
                           obstime=flloc_time, observer=suvi_map.observer_coordinate, frame=frames.Heliocentric)

    # Plot flare
    flloc_rt_sc = flloc_rt_sc.transform_to(suvi_map.coordinate_frame)
    ax_rt.plot_coord(flloc_rt_sc, color=c, marker='*', markersize=4, alpha=0.85)

# Add legend for features
handles, labels = legend_handles()
plt.legend(handles, labels, loc=1, framealpha=1)

plt.show()

Helioprojective Cartesian Coordinates

# Plot the composite image
fig_hpc = plt.figure(figsize=(10, 10))
fig_hpc.tight_layout(h_pad=-50)
fig_hpc.suptitle(f'Helioprojective Cartesian Coordinates, GOES-{goes}')
ax_hpc = fig_hpc.add_subplot(projection=suvi_map)
suvi_map.plot(axes=ax_hpc, clip_interval=(0.1, 99.9) * u.percent)
num_points = 100

x_limits = ax_hpc.get_xlim()
y_limits = ax_hpc.get_ylim()

# Extract flare location from netCDF
flloc_hpc = flloc_nc['flloc_hpc'][:][0]
wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)

peak_hpc = flloc_nc['peak_loc_hpc'][:][0]

colors = cm.Blues_r(np.linspace(0., 0.3, len(flloc_hpc)))

for flare, c in zip(flloc_hpc, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)
    wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)

    x = float(flare[wavelength_i][0])
    y = float(flare[wavelength_i][1])

    # Convert flare location into skycoord object in heliographic stonyhurst frame
    flloc_hpc_sc = SkyCoord(x * u.arcsec, y * u.arcsec, obstime=flloc_time, observer=suvi_map.observer_coordinate,
                            frame=frames.Helioprojective)

    # Set params
    line_resolution = 100

    # Draw x dimension of crosshairs
    x_location = np.repeat(x, line_resolution) * u.arcsec
    y_line = np.linspace(-1600, 1600, line_resolution) * u.arcsec
    # Draw y dimension of crosshairs
    y_location = np.repeat(y, line_resolution) * u.arcsec
    x_line = np.linspace(-1600, 1600, line_resolution) * u.arcsec

    # Transform skycoord object to the suvi sunpy map's coordinate frame
    flloc_hpc_sc = flloc_hpc_sc.transform_to(suvi_map.coordinate_frame)

    ax_hpc.plot(x_line.to(u.deg), y_location.to(u.deg), color=c, transform=ax_hpc.get_transform("world"))
    ax_hpc.plot(x_location.to(u.deg), y_line.to(u.deg), color=c, transform=ax_hpc.get_transform("world"))

    # Plot flare
    ax_hpc.plot_coord(flloc_hpc_sc, color=c, marker='o', fillstyle='none', markersize=8, alpha=0.7)

for flare, c in zip(peak_hpc, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)
    wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)

    x = float(flare[wavelength_i][0])
    y = float(flare[wavelength_i][1])

    # Convert flare location into skycoord object in heliographic stonyhurst frame
    flloc_hpc_sc = SkyCoord(x * u.arcsec, y * u.arcsec, obstime=flloc_time, observer=suvi_map.observer_coordinate,
                            frame=frames.Helioprojective)

    # Transform skycoord object to the suvi sunpy map's coordinate frame
    flloc_hpc_sc = flloc_hpc_sc.transform_to(suvi_map.coordinate_frame)

    # Plot flare
    ax_hpc.plot_coord(flloc_hpc_sc, color=c, marker='*', markersize=4, alpha=0.85)

# add legend for features
handles, labels = legend_handles()
plt.legend(handles, labels, loc=1, framealpha=1)
ax_hpc.set_xlim(x_limits)
ax_hpc.set_ylim(y_limits)

plt.show()

Heliocentric Cartesian Coordinates

# Plot the composite image
fig_hcc = plt.figure(figsize=(10, 10))
fig_hcc.tight_layout(h_pad=-50)
fig_hcc.suptitle(f'Heliocentric Cartesian Coordinates, GOES-{goes}')
ax_hcc = fig_hcc.add_subplot(projection=suvi_map)
suvi_map.plot(axes=ax_hcc, clip_interval=(0.1, 99.9) * u.percent)
num_points = 100

x_limits = ax_hcc.get_xlim()
y_limits = ax_hcc.get_ylim()

# Extract flare location from netCDF
flloc_hcc = flloc_nc['flloc_hcc'][:][0]
wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)

peak_hcc = flloc_nc['peak_loc_hcc'][:][0]

colors = cm.PiYG(np.linspace(0., 0.3, len(flloc_hcc)))

for flare, c in zip(flloc_hcc, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)
    wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)

    x = float(flare[wavelength_i][0])
    y = float(flare[wavelength_i][1])
    z = float(flare[wavelength_i][2])

    # Convert flare location into skycoord object in heliographic stonyhurst frame
    flloc_hcc_sc = SkyCoord(x * u.AU, y * u.AU, z * u.AU, obstime=flloc_time, observer=suvi_map.observer_coordinate,
                            frame=frames.Heliocentric)

    # Plot flare
    ax_hcc.plot_coord(flloc_hcc_sc, color=c, marker='o', fillstyle='none', markersize=8, alpha=0.7)

    # Set params
    line_resolution = 100

    # sunpy.sun.constants.AU.value*1E-3
    # Draw x dimension of crosshairs
    x_location = (np.repeat(x, line_resolution))
    y_line = (np.linspace(-0.00764747, 0.00765689, line_resolution))
    # Draw y dimension of crosshairs
    y_location = (np.repeat(y, line_resolution))
    x_line = (np.linspace(-0.00764747, 0.00765689, line_resolution))

    # Set up the coord instance for the X-position of a feature (flare/peak/centroid, etc)
    coord_x = SkyCoord(astropy.coordinates.representation.CartesianRepresentation(x_location * u.AU, y_line * u.AU,
                                                                                  np.repeat(0,
                                                                                            line_resolution) * u.AU),
                       obstime=flloc_time, observer=suvi_map.observer_coordinate, frame=frames.Heliocentric)

    # Set up the coord instance for the Y-position of a feature (flare/peak/centroid, etc)
    coord_y = SkyCoord(astropy.coordinates.representation.CartesianRepresentation(x_line * u.AU, y_location * u.AU,
                                                                                  np.repeat(0,
                                                                                            line_resolution) * u.AU),
                       obstime=flloc_time, observer=suvi_map.observer_coordinate, frame=frames.Heliocentric)
    ax_hcc.plot_coord(coord_x, color=c)
    ax_hcc.plot_coord(coord_y, color=c)
    # Transform skycoord object to the suvi sunpy map's coordinate frame
    # flloc_hcc_sc = flloc_hcc_sc.transform_to(suvi_map.coordinate_frame)

for flare, c in zip(peak_hcc, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)
    wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)

    x = float(flare[wavelength_i][0])
    y = float(flare[wavelength_i][1])
    z = float(flare[wavelength_i][2])

    # Convert flare location into skycoord object in heliographic stonyhurst frame
    flloc_hcc_sc = SkyCoord(x * u.AU, y * u.AU, z * u.AU, obstime=flloc_time, observer=suvi_map.observer_coordinate,
                            frame=frames.Heliocentric)

    # Transform skycoord object to the suvi sunpy map's coordinate frame
    flloc_hcc_sc = flloc_hcc_sc.transform_to(suvi_map.coordinate_frame)

    # Plot flare
    ax_hcc.plot_coord(flloc_hcc_sc, color=c, marker='*', markersize=4, alpha=0.85)

# add legend for features
handles, labels = legend_handles()
plt.legend(handles, labels, loc=1, framealpha=1)

ax_hcc.set_xlim(x_limits)
ax_hcc.set_ylim(y_limits)

plt.show()

Heliocentric Radial Coordinates

# Plot the SUVI SunPy map
fig_hcr = plt.figure(figsize=(10, 10))
fig_hcr.suptitle(f'Heliocentric Radial Coordinates, GOES-{goes}')
fig_hcr.tight_layout(h_pad=-50)
ax_hcr = fig_hcr.add_subplot(projection=suvi_map)
suvi_map.plot(axes=ax_hcr, clip_interval=(0.1, 99.9) * u.percent)

# Extract flare location from netCDF
flloc_hcr = flloc_nc['flloc_hcr'][:][0]
wavelength_i = list(flloc_nc['wavelength'][:]).index(wavelength)

peak_rt = flloc_nc['peak_loc_hcr'][:][0]

colors = cm.BuPu(np.linspace(0.85, 1., len(flloc_hcr)))

for flare, c in zip(flloc_hcr, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)
    # Need to rotate theta by 90 to get psi, according to SkyCoord cylindrical specifications
    rho = flare[wavelength_i][0]
    psi = flare[wavelength_i][1] + 90.
    z = flare[wavelength_i][2]

    # Create SkyCoord object in heliocentric frame
    flloc_hcr_sc = SkyCoord(rho=rho * u.solRad, psi=psi * u.deg, z=z * u.solRad, representation_type='cylindrical',
                            obstime=flloc_time, observer=suvi_map.observer_coordinate, frame=frames.Heliocentric)

    # Plot r/psi 'crosshairs'
    heliocentric_frame = frames.Heliocentric(obstime=suvi_map.date, observer=suvi_map.observer_coordinate)
    r_vector = SkyCoord(rho=flloc_hcr_sc.rho, psi=np.linspace(-360, 0, num_points) * u.deg, z=z * u.solRad,
                        frame=heliocentric_frame, representation_type='cylindrical')
    constant_psi = SkyCoord(rho=np.linspace(0, rho, num_points) * u.solRad, psi=flloc_hcr_sc.psi, z=z * u.solRad,
                            frame=heliocentric_frame, representation_type='cylindrical')

    ax_hcr.plot_coord(r_vector, color=c, alpha=0.7)
    ax_hcr.plot_coord(constant_psi, color=c, alpha=0.7)

    # Plot flare
    flloc_hcr_sc = flloc_hcr_sc.transform_to(suvi_map.coordinate_frame)
    ax_hcr.plot_coord(flloc_hcr_sc, color=c, marker='o', fillstyle='none', markersize=8, alpha=0.7)

for flare, c in zip(peak_rt, colors):
    flare[flare == -9999.] = np.nan
    flare = flare.filled(fill_value=np.nan)
    # Need to rotate theta by 90 to get psi, according to SkyCoord cylindrical specifications
    rho = flare[wavelength_i][0]
    psi = flare[wavelength_i][1] + 90.
    z = flare[wavelength_i][2]

    # Create SkyCoord object in heliocentric frame
    flloc_hcr_sc = SkyCoord(rho=rho * u.solRad, psi=psi * u.deg, z=z * u.solRad, representation_type='cylindrical',
                            obstime=flloc_time, observer=suvi_map.observer_coordinate, frame=frames.Heliocentric)
    # Plot flare
    flloc_hcr_sc = flloc_hcr_sc.transform_to(suvi_map.coordinate_frame)
    ax_hcr.plot_coord(flloc_hcr_sc, color=c, marker='*', markersize=4, alpha=0.85)

# add legend for features
handles, labels = legend_handles()
plt.legend(handles, labels, loc=1, framealpha=1)

plt.show()