Air temperature doesn't look right in GFS reference forecast

I was looking at air temperature forecasts from GFS for Dec 23-24 and noticed that the forecasted values were surprisingly high. Opening the NetCDF file in NOAA's WCT, I'm seeing different values that are much lower (and more believable). I'm wondering if SFA is processing air temperature incorrectly.

As an example, I'm using the GFS forecast initialized at 2022-12-22 06:00Z and a valid time of 2022-12-24 00:00Z, for a location near Atlanta, GA (33.75, -84.30). I used the workshop reference_forecasts.ipynb as a starting point. Here's the NetCDF file zipped: sfa_test.zip.

Opening the .nc file in NOAA WCT, I get -8 deg C (265.1 K). But running through sfa, I get about -2 C. For reference, the actual temperature ended up being about -11.5 C (https://www.wunderground.com/history/daily/us/ga/atlanta/KATL/date/2022-12-24).

I tried a few other locations and got similar results. Maybe I'm doing something wrong, but it seems like the forecasted temperature never gets low enough, so something seems off.

Viewing the file in NOAA's WCT:

SFA results (code is included below):

[EDIT: THIS CODE CONTAINS A TYPO!!]

import datetime
from functools import partial
from pathlib import Path

import numpy as np
import pandas as pd

from bokeh.core.properties import value
from bokeh.io import output_notebook
from bokeh.layouts import gridplot
from bokeh.models import ColumnDataSource
from bokeh.plotting import figure, show
from bokeh import palettes
from bokeh.palettes import Category10_10 as PALETTE
TOOLS = "pan,box_zoom,xwheel_zoom,reset,save"
output_notebook()

from solarforecastarbiter import datamodel
from solarforecastarbiter.io import nwp
from solarforecastarbiter.io.api import APISession, request_cli_access_token

# location of netcdf files 
base_path = 'C:/Users/willh/Downloads/sfa_test' # path to folder used for "solararbiter fetchnwp"
# define file loading function that knows where to find the files
load_forecast = partial(nwp.load_forecast, base_path=base_path)

def plot_nwp(ghi, dni, dhi, air_temperature, wind_speed, ac_power):
    """Make a plot of the variables from processed NWP forecasts"""
    # data prep
    if ac_power is None:
        ac_power = pd.Series(np.nan, index=ghi.index)
    data = dict(zip(('GHI', 'DNI', 'DHI', 'Air temperature', 'Wind speed', 'AC power', 'index'),
                    (ghi, dni, dhi, air_temperature, wind_speed, ac_power, ghi.index)))
    source = ColumnDataSource(data)

    figs = []
    fig_kwargs = dict(tools=TOOLS, x_axis_type="datetime", plot_height=200)
    palette = iter(PALETTE)

    # irradiance figure
    fig1 = figure(**fig_kwargs)
    for name in ('GHI', 'DNI', 'DHI'):
        fig1.line(x='index', y=name, source=source, legend_label=name, color=next(palette), line_width=2)
    fig1.yaxis.axis_label = "Irradiance (W/m^2)"
    fig1.legend.location = "top_left"
    figs.append(fig1)

    # air temperature figure
    fig_kwargs['x_range'] = fig1.x_range  # link x-zoom
    fig2 = figure(**fig_kwargs)
    fig2.line(x='index', y='Air temperature', source=source, legend_label="Air temperature", color=next(palette), line_width=2)
    fig2.yaxis.axis_label = "Air temperature (C)"
    fig2.legend.location = "top_left"
    figs.append(fig2)

    # wind speed figure
    fig3 = figure(**fig_kwargs)
    fig3.line(x='index', y='Wind speed', source=source, legend_label="Wind speed", color=next(palette), line_width=2)
    fig3.yaxis.axis_label = "Wind speed (m/s)"
    fig3.legend.location = "top_left"
    figs.append(fig3)

    # ac power figure
    if not ac_power.dropna().empty:
        fig4 = figure(**fig_kwargs)
        fig4.line(x='index', y='AC power', source=source, legend_label="AC power", color=next(palette), line_width=2)
        fig4.yaxis.axis_label = "AC power (MW)"
        fig4.legend.location = "top_left"
        figs.append(fig4)

    figs[-1].xaxis.axis_label = 'Time (UTC)'

    grid = gridplot(figs, ncols=1, plot_width=800)
    return grid

# define coordinates (Atlanta, GA)
latitide = 33.75 
longitude = -84.30
elevation = 225 

# define a Site object
site = datamodel.Site(
    name='Atlanta, GA',
    latitude=latitude,
    longitude=longitude,
    elevation=elevation,
    timezone='America/New York'
)
# Complete GFS is available 5 hours after its initialization time
run_time = pd.Timestamp('20221222T1100Z')
issue_time = pd.Timestamp('20221222T0700Z')

# define a Forecast object
forecast = datamodel.Forecast(
    name='GFS 7day Weather Fcast',
    site=site,
    variable='ghi',
    interval_label='instant',
    interval_value_type='interval_mean',
    interval_length='1h',
    issue_time_of_day=datetime.time(hour=7),
    run_length=pd.Timedelta('24h')*7,
    lead_time_to_start=pd.Timedelta('0h')
)

# select the model
model = models.gfs_quarter_deg_to_hourly_mean

# tell the model where to find the NWP data
model_wrapped = partial(model, load_forecast=load_forecast)

# load and process the NWP data. ac_power is None because site is not a power plant.
ghi, dni, dhi, air_temperature, wind_speed, ac_power = main.run_nwp(forecast, model_wrapped, run_time, issue_time)

grid = plot_nwp(ghi, dni, dhi, air_temperature, wind_speed, ac_power)
show(grid)

SolarArbiter / solarforecastarbiter-core

Air temperature doesn't look right in GFS reference forecast #801