Occasional negative or excessive positive irradiance in GEFS

[williamhobbs]

latitudes = [34.82]
longitudes = [-89.43]

gefs_test = get_solar_forecast_ensemble(
    latitudes,
    longitudes,
    init_date='2026-03-11 00:00',
    run_length=24,
    lead_time_to_start=9,
    model='gefs',
    priority='aws',
    )

point=0
# plot all members
for member in range(1,31):
    gefs_test[(gefs_test['point']==point) &
        (gefs_test['member']==member)]['ghi'].plot(alpha=0.5, color='C0')

# plot member 12 in red
gefs_test[(gefs_test['point']==point) &
        (gefs_test['member']==12)]['ghi'].plot(color='red')

gives

Similar for 2025-11-29

And multiple issues with 2025-08-12

[williamhobbs]

But then re-running 2025-08-12, everything looks ok:

latitudes = [34.82]
longitudes = [-89.43]

gefs_test = get_solar_forecast_ensemble(
    latitudes,
    longitudes,
    init_date='2025-08-12 00:00',
    run_length=21,
    lead_time_to_start=9,
    model='gefs',
    priority='aws',
    )

point=0
# plot all members
for member in range(1,31):
    gefs_test[(gefs_test['point']==point) &
        (gefs_test['member']==member)]['ghi'].plot(alpha=0.5, color='C0')

but not for 2025-11-29

so, it seems there are two issues:

• negative irradiance
• some issue with files not always downloading/processing correctly (maybe related to check number of ensemble members with IFS and AIFS #28? Incomplete or erroneous but still usable grib files?)

[williamhobbs]

This could be related to #26.

[williamhobbs]

I now think this is just an issue in GEFS, likely related to precision/rounding. Code with a simple example/illustration checking this is below.

Trying to find a place to ask about this via NOAA-GSL/GSL-Regional-Model-Forum#55.

import hefty
import pandas as pd
import matplotlib.pyplot as plt
from herbie.core import Herbie
import pvlib
import xarray as xr

# hefty sometimes returns negative GHI from GEFS. See https://github.com/williamhobbs/hefty/issues/29. 

# Start with a single example:

latitudes = [34.82]
longitudes = [-89.43]

gefs_test_single_mem = hefty.solar.get_solar_forecast(
    latitudes,
    longitudes,
    init_date='2025-11-29 00:00',
    run_length=18,
    lead_time_to_start=9,
    model='gefs',
    member=4,
    priority='aws',
    )

point=0
gefs_test_single_mem[(gefs_test_single_mem['point']==point)]['ghi'].plot()
plt.legend()
plt.show()

# Now replicate what happens in `hefty.solar.get_solar_forecast` to troubleshoot steps along the way:

latitude = [34.82]
longitude = [-89.43]

date = '2025-11-29 00:00'
run_length = 18
lead_time_to_start = 9
fxx_max = run_length + lead_time_to_start
fxx_range = range(lead_time_to_start, fxx_max + 1, 3)
model = 'gefs'
member = 4
product = 'atmos.25'  # 0.25 deg, 'pgrb2.0p25'
search_str = 'DSWRF|:TMP:2 m above|[UV]GRD:10 m above'
# search_str = 'DSWRF'
priority = 'aws'

i = []
for fxx in fxx_range:
    ds = Herbie(
        date,
        model=model,
        product=product,
        fxx=fxx,
        member=member,
        priority=priority
        ).xarray(search_str)
    
    ds = xr.merge(ds, compat='override')

    # use pick_points for single point or list of points
    i.append(
        ds.herbie.pick_points(
            pd.DataFrame(
                {
                    "latitude": latitude,
                    "longitude": longitude,
                }
            )
        )
    )

    ts = xr.concat(i, dim="valid_time")  # concatenate

    # df_temp = ts.to_dataframe()[['sdswrf', 't2m', 'si10']]
    df_temp = ts.to_dataframe()[['sdswrf']]
    # add timezone
    df_temp = df_temp.tz_localize('UTC', level='valid_time')

num_sites = len(latitude)

# work through sites
dfs = {}  # empty list of dataframes
for j in range(num_sites):
    df = df_temp[df_temp.index.get_level_values('point') == j]
    df = df.droplevel('point')

    loc = pvlib.location.Location(
        latitude=latitude[j],
        longitude=longitude[j],
        tz=df.index.tz
        )

    if model in {'gfs', 'gefs'}:
        # for gfs and gefs ghi: we have to "unmix" the rolling average
        # irradiance that resets every 6 hours
        mixed = df[['sdswrf']].copy()
        mixed['hour'] = mixed.index.hour
        mixed['hour'] = mixed.index.hour
        mixed['hour_of_mixed_period'] = ((mixed['hour'] - 1) % 6) + 1
        mixed['sdswrf_prev'] = mixed['sdswrf'].shift(
            periods=1,
            fill_value=0
            )
        mixed['int_len'] = mixed.index.diff().total_seconds().values / 3600

        # set the first interval length:
        if model == 'gfs' and lead_time_to_start >= 120:
            mixed.loc[mixed.index[0], 'int_len'] = 1
        else:
            mixed.loc[mixed.index[0], 'int_len'] = 3
        unmixed = ((mixed['hour_of_mixed_period'] * mixed['sdswrf']
                    - (mixed['hour_of_mixed_period'] - mixed['int_len'])
                    * mixed['sdswrf_prev']) / mixed['int_len'])
        df['ghi'] = unmixed

print(df_temp)
print(df)

gives this output:

                                 sdswrf
valid_time                point        
2025-11-29 09:00:00+00:00 0         0.0
2025-11-29 12:00:00+00:00 0         0.0
2025-11-29 15:00:00+00:00 0        50.0
2025-11-29 18:00:00+00:00 0        80.0
2025-11-29 21:00:00+00:00 0        50.0
2025-11-30 00:00:00+00:00 0        20.0
2025-11-30 03:00:00+00:00 0         0.0
                           sdswrf    ghi
valid_time                              
2025-11-29 09:00:00+00:00     0.0    0.0
2025-11-29 12:00:00+00:00     0.0    0.0
2025-11-29 15:00:00+00:00    50.0   50.0
2025-11-29 18:00:00+00:00    80.0  110.0
2025-11-29 21:00:00+00:00    50.0   50.0
2025-11-30 00:00:00+00:00    20.0  -10.0
2025-11-30 03:00:00+00:00     0.0    0.0

See https://earthscience.stackexchange.com/questions/24235/understanding-gfs-parameter-averaging-times, https://www.nco.ncep.noaa.gov/pmb/products/gens/gep01.t00z.pgrb2a.0p50.f003.shtml, https://solarforecastarbiter-core.readthedocs.io/en/1.0.4/generated/solarforecastarbiter.reference_forecasts.forecast.unmix_intervals.html

The forecast valid hour 21 (2025-11-29 21:00:00+00:00) is the average for hours 18 - 21, 50 W/m^2. Call this A1.

The forecast valid hour 24 (2025-11-30 00:00:00+00:00) is the average for hours 18 - 24, 20 W/m^2. Call this A2.

We want the average for each interval, 18 - 21 (call this I1) and 21 - 24 (call this I2).

I1 = A1

I2 = (2 * A2) - I1

The average for hours 21 - 24, I2, is (2 * 20) - 50 = -10.

[williamhobbs]

NOAA NCEP EMC confirmed this is an issue, likely related to GEFS precision, and they are looking into a possible solution.

I mentioned this in a semi-relevant thread on the UCAR WRF&MPAS forum: https://forum.mmm.ucar.edu/threads/resolved-issues-with-gefs-ics-causing-wrf-to-crash-in-first-model-time-step.12770/#post-61309

And on LinkedIn: https://www.linkedin.com/posts/will-hobbs-93215023_if-you-use-gefs-for-solar-energy-applications-activity-7457476251046522881-N2WF/

And here's another illustration:

from herbie import Herbie
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

date = '2025-11-29 06:00'
model = 'gefs'
product = 'atmos.25'
member = 'c00'  # could be 'c00', 'p04', or others...
search_str = 'DSWRF'
data_var = 'sdswrf'

# irradiance at forecast hour 21, which is the average from hour 18 to 21
ds1 = Herbie(date, model=model, product=product, fxx=21, member=member).xarray(search_str)

# irradiance at forecast hour 24, which is the average from hour 18 to 24
ds2 = Herbie(date, model=model, product=product, fxx=24, member=member).xarray(search_str)

# invert the bottom half of the reversed RdBu cmap to get more contrast below zero
base_cmap = plt.get_cmap('RdBu_r')
colors = base_cmap(np.linspace(0, 1, 256))
midpoint = 256 // 2
colors[:midpoint] = colors[:midpoint][::-1]
inverted_bottom_cmap = ListedColormap(colors)

# calculate irradiance from hour 21 to 24 ("unmix" the rolling averages)
i2 = (2 * ds2[data_var]) - ds1[data_var]

# make a plot
i2.plot(figsize=(10,4), cmap=inverted_bottom_cmap)
plt.show()

# convert to 1d array of values
i2_vals = i2.values.reshape(-1)

# split into positive and negative
# drop zeros to make the histogram easier to read
i2_pos = np.where(i2_vals > 0, i2_vals, np.nan)
i2_neg = np.where(i2_vals < 0, i2_vals, np.nan)

# make a histogram
plt.figure(figsize=(10,4))
plt.hist(i2_pos, bins=range(-30,1100,10), color='C3', label='positive')
plt.hist(i2_neg, bins=range(-30,1100,10), color='C0', label='negative')
plt.legend()
plt.xlabel('DSWRF, W/m-2 (excluding zero)')
plt.ylabel('count')
plt.show()