Download and Plot SGPS Proton Flux Data

Download, read, and plot proton flux from SGPS L2 data.

Import modules

import os
import requests
import netCDF4 as nc
import numpy as np
import cftime
import matplotlib.pyplot as plt
from matplotlib import colors, gridspec
import matplotlib.dates as mdates

Download, read, and get relevant data from files

files = ["sci_sgps-l2-avg5m_g19_d20260118_v3-0-3.nc",
         "sci_sgps-l2-avg5m_g19_d20260119_v3-0-3.nc",
         "sci_sgps-l2-avg5m_g19_d20260120_v3-0-3.nc"]

url_path = "https://data.ngdc.noaa.gov/platforms/solar-space-observing-satellites/goes/goes19/l2/data/sgps-l2-avg5m/2026/01/"

proton_diff_flux_west = []
proton_diff_flux_east = []
proton_int_flux = []
time = []

for filename in files:
    # Download `filename` if it does not exist locally
    if not os.path.exists(filename):
        with open(filename, "wb") as f:
            r = requests.get(url_path + filename)
            f.write(r.content)

    # Load data from file
    sgps_data = nc.Dataset(filename)

    # Get flux data
    proton_diff_flux_west.append(sgps_data['AvgDiffProtonFlux'][:, 0, :])
    proton_diff_flux_east.append(sgps_data['AvgDiffProtonFlux'][:, 1, :])
    proton_int_flux.append(sgps_data['AvgIntProtonFlux'][:])

    # Get timestamps
    # Two different time variable names may have been used
    try:
        time.append(sgps_data['L2_SciData_TimeStamp'][:])
    except IndexError:
        time.append(sgps_data['time'][:])

# Convert list of arrays into single array
proton_diff_flux_west = np.ma.concatenate(proton_diff_flux_west)
proton_diff_flux_east = np.ma.concatenate(proton_diff_flux_east)
proton_int_flux = np.ma.concatenate(proton_int_flux)
time = np.ma.concatenate(time)

# replace zeros with nans
proton_diff_flux_west = np.where(proton_diff_flux_west < 1.e-12, np.nan, proton_diff_flux_west)
proton_diff_flux_east = np.where(proton_diff_flux_east < 1.e-12, np.nan, proton_diff_flux_east)
proton_int_flux = np.where(proton_int_flux < 1.e-12, np.nan, proton_int_flux)

# Convert J2000 time to python datetime
time = cftime.num2pydate(time[:], sgps_data['time'].units)

Plot SGPS flux

plt.figure(1, figsize=[12, 8], layout='constrained')
gridspec.GridSpec(5, 1)

# Plot differential flux
NUM_DIFF_CHANNELS = 13

chan_colors = [[.8, 0., 0.], colors.to_rgba('orangered')[0:3], colors.to_rgba('orange')[0:3],
               [.95, .95, .0], colors.to_rgba('greenyellow')[0:3], colors.to_rgba('yellowgreen')[0:3],
               colors.to_rgba('green')[0:3], [0., .78, .7], [0., .6, .88], [0., 0., .9],
               colors.to_rgba('violet')[0:3], [0.8, .1, 1.0], [0.49411765, 0., 0.70980392]]

chan_labels = ['P1 (1.0-1.9 MeV)', 'P2A (1.9-2.3 MeV)', 'P2B (2.3-3.4 MeV)', 'P3 (3.4-6.5 MeV)',
               'P4 (6.5-12 MeV)', 'P5 (12-25 MeV)', 'P6 (25-40 MeV)', 'P7 (40-80 MeV)',
               'P8A (83-99 MeV)', 'P8B (99-118 MeV)', 'P8C (118-150 MeV)', 'P9 (150-275 MeV)',
               'P10 (275-500 MeV)']

# SGPS-X (west)
ax1 = plt.subplot2grid((5, 1), (0, 0), colspan=1, rowspan=2)
for i in range(NUM_DIFF_CHANNELS):
    plt.plot(time, proton_diff_flux_west[:, i], color=chan_colors[i], label=chan_labels[i])
textstr = 'G19 SGPS-X (west looking field-of-view)'
ax1.text(0.042, 0.97, textstr, transform=ax1.transAxes, fontsize=14, verticalalignment='top')
ax1.xaxis.set_minor_locator(mdates.HourLocator())
ax1.xaxis.set_major_locator(mdates.DayLocator())
ax1.tick_params(which = 'minor', length=3)
ax1.tick_params(which = 'major', length=5)
ax1.set_xlim(time[0], time[-1])
ax1.tick_params(labelbottom=False)
plt.yscale('log')
plt.ylim([1.e-8, 1.e3])
plt.ylabel('protons/cm$^2$-s-sr-keV')
ax1.legend(bbox_to_anchor=(1.28, 1.2), loc='upper right', prop={'size': 12})

# SGPS+X (east)
ax2 = plt.subplot2grid((5, 1), (2, 0), colspan=1, rowspan=2)
for i in range(NUM_DIFF_CHANNELS):
    plt.plot(time, proton_diff_flux_east[:, i], color=chan_colors[i], label=chan_labels[i])
textstr = 'G19 SGPS+X (east looking field-of-view)'
ax2.text(0.042, 0.97, textstr, transform=ax2.transAxes, fontsize=14, verticalalignment='top')
ax2.xaxis.set_minor_locator(mdates.HourLocator())
ax2.xaxis.set_major_locator(mdates.DayLocator())
ax2.tick_params(which = 'minor', length=3)
ax2.tick_params(which = 'major', length=5)
ax2.set_xlim(time[0], time[-1])
ax2.tick_params(labelbottom=False)
plt.yscale('log')
plt.ylim([1.e-8, 1.e3])
plt.ylabel('protons/cm$^2$-s-sr-keV')

# Plot integral flux
ax3 = plt.subplot2grid((5, 1), (4, 0), colspan=1, rowspan=1)
plt.plot(time, proton_int_flux[:, 1], color='grey', label='SGPS+X P11 (>500 MeV)')
plt.plot(time, proton_int_flux[:, 0], color='k', label='SGPS-X P11 (>500 MeV)')
textstr = 'G19 SGPS P11'
ax3.text(0.042, 0.94, textstr, transform=ax3.transAxes, fontsize=14, verticalalignment='top')
ax3.xaxis.set_minor_locator(mdates.HourLocator())
ax3.xaxis.set_major_locator(mdates.DayLocator())
ax3.xaxis.set_major_formatter(mdates.DateFormatter('%m/%d'))
ax3.tick_params(which = 'minor', length=3)
ax3.tick_params(which = 'major', length=5)
ax3.set_xlim(time[0], time[-1])
plt.xlabel('2026 (UT hours)')
plt.ylabel('protons/cm$^2$-s-sr')
plt.yscale('log')
ylims = ax3.get_ylim()
plt.ylim([ylims[0], 1.4*ylims[1]])

ax3.legend(bbox_to_anchor=(1.28, 1), prop={'size': 10})

plt.show()