pyGHCND is a Python library for collecting and plotting weather data from NOAA's National Centers for Environmental Information (NCEI) Global Historical Climatology Network daily (GHCNd) datasets. These data are accessed using NCEI's Climate Data Online (CDO) API and filtered based on the corresponding attribute flags. Some statistical analyses on the data trends are then performed. The raw and statistical data tables (Pandas DataFrames) are saved to disk. Several plotting functions are provided to facilitate visualizing trends in these data over time.
- python >= 3.10
- requests >= 2.20
- pandas >= 1.5
- pyarrow >= 8.0
- scipy >= 1.10
- matplotlib >= 3.7
- openpyxl >= 3.0
- tqdm >= 4.60
- Caveats and Notes on GHCNd Data
- Simple Usage Example
- Basic Data Collection
- Raw Data Tables
- Statistical Analyses Table
- Plotting Function Overview
- Yearly Temperature Plot
- Recent Temperature Deviation Plot
- Precipitation Plot
- Daily Temperature Trend Plot
- Daily Trends Overview Plot
- Sorting Daily Trends
The GHCNd data are raw weather measurement data, and as such, they may suffer
from deviations caused by a variety of factors, such as changes in measurement
location/time/equipment/etc. The UCAR website has a nice whitepaper on the
GHCNd
data
that highlights some of these issues. Except for filtering measurements that
are flagged as being problematic, pyGHCND does not currently apply any
adjustments to compensate for these additional factors. With that being said,
the statistical analyses performed by this module need to be interpreted with
some caution.
The GHCNd network consists of a large number of weather stations that vary in
their yearly coverage and data quality. The US Historical Climate
Network
is a subset of the GHCNd stations that are deemed to contain data of higher
quality. The stations that are part of this network are flagged with 'HCN' in
the GHCNd station
list, for
example. These stations are good starting points for data analysis using
pyGHCND.
The following code block is a complete example of the data collection, analysis, and plotting using pyGHCND. The weather data in this example comes from the Missoula International Airport. With one exception, this code be saved into a .py file and executed from the command line, for example. The only additional step is requesting an API token from the Climate Data Online website and saving it into a file called 'token.txt'. A description of each line of code is found in the sections below.
from pyghcnd import GHCND, mplvis
import matplotlib.pyplot as plt
plt.style.use('ggplot')
with open('token.txt') as token_file:
token = token_file.read().strip()
temps = GHCND('USW00024153', token, data_folder='MissoulaMSO')
temps.update_data(status=True)
temps.stats.to_excel('MSO_stats.xlsx')
mplvis.plot_temp(temps)
mplvis.plot_temp_diffs(temps)
mplvis.plot_prcp(temps, ptype='SNPR')
mplvis.plot_daily_temp(temps, day=1, month=1)
mplvis.plot_daily_trends(temps)
for cat in ('TMIN', 'TMAX'):
slope_sorted = temps.stats[cat]\
.sort_values(by='-log_p*abs_slope', ascending=False)
first_slope = slope_sorted.iloc[0]
month, day = first_slope.name
mplvis.plot_daily_temp(temps, month=month, day=day, temp_type=cat)
The core element of pyGHCND is a simple object class called GHCND.
from pyghcnd import GHCND
To create a GHCND object instance, two pieces of data are necessary: an API
token and a GHCNd station ID. The API token is freely available from the
Climate Data Online website. The
desired station ID can be found using the NCEI mapping
tool, for example. The full
station ID will be something like "USW00024153", which corresponds to
the Missoula, Montana, airport. The station ID with the preceding "GHCND:"
-- e.g. "GHCND:USW00024153" -- is also acceptable.
The GHCND class has the following initialization call signature.
GHCND(stationid, token, store_type='parquet', data_folder='.')
The stationid and token are the strings representing the station ID and
API token as described above. The store_type keyword argument determines the
storage format for the raw and statistical data. The default ('parquet')
uses Parquet binary files, which were developed
specifically for column-oriented data (e.g. DataFrames). The data_folder
keyword argument sets the folder where data and plots will be stored. The
default is the current directory ('.'). Note: the stored data
files and plots have a fixed naming scheme, so downloading multiple datasets
into the same folder will overwrite previous data sets.
In the following example, a GHCND object instance for the historical weather
data at the Missoula International
Airport
will be created. This assumes that the API token has been saved into a
text file called 'token.txt' in the current directory (i.e. not the
data_folder directory). Note: this code will make an API call, so an
external internet connection is necessary.
with open('token.txt') as token_file:
token = token_file.read().strip()
temps = GHCND('USW00024153', token, data_folder='MissoulaMSO')
On initialization, the GHCND object performs several operations. First, an API
call is made to collect station metadata, such as the start and end dates for
the measurements at that location -- these are saved as the start_date and
end_date attributes, respectively. (The raw metadata is stored as the
private _station_info dictionary.) It then checks for the existence of the
given data folder. If it is not found, the folder is created. However, if the
folder is present, a check is performed to determine if previously
downloaded/processed data is present, and if so, these data are loaded into
the object instance (more info below). This makes it possible to manipulate
previously downloaded data or append newly downloaded information without
having to re-download the entire dataset.
The update_data method is used to download new data using the CDO API. This
method accepts only a single Boolean keyword argument: status. The default
for this argument is False, which means that the data will be downloaded
silently. If this argument is True, then a download status bar will be
printed to the screen, using the tqdm
package. The example below will download the
data with a progress bar.
temps.update_data(status=True)
The download process can be slow. The CDO API restrict a single data request to a total of 1000 data points, and only 5 requests can be made per second. Each year may consist of several thousand data points, so if data is being collected for long time intervals, the total number of requests will be substantial. For example, the Missoula airport has daily weather data from 1948 to the present. To protect against broken connections, a temporary file (pickle format) of raw download data is periodically saved in the data folder during the download process, and in the event that the download must be restarted, the temporary file will first be loaded as a starting point.
After completion of the data download, one private DataFrame attribute --
_raw_full -- will be created and saved to disk. A selection of this
DataFrame as rendered in a Jupyter notebook is shown below. This DataFrame is
stored as a "hidden" attribute because it is probably not necessary/desired to
work with this directly; however, it is instructive to view this table to
understand what the raw data looks like as it comes from the API and how it is
filtered, as described below. The API data is initially delivered in JSON
format, and after complete download, these data are reformatted into a Pandas
DataFrame indexed by the collection date. Temperatures are converted to
degrees Fahrenheit, and precipitation is converted to inches. The columns of
this DataFrame are the data types, such as max temp (TMAX), min temp (TMIN),
etc. Under these data type columns, there are several sub-columns for the
recorded value ("value") and the measurement attributes, such as the data
quality ("qflag"). The data types and attributes are described in the GHCNd
documentation. This DataFrame is saved to disk using the store_type selected
above.
This full data set is then filtered to remove measurements flagged as being
problematic -- this filtered DataFrame is called raw. An example of this
DataFrame, as rendered in a Jupyter notebook, is shown below. Bad measurements
in the _raw_full DataFrame, as determined by their quality attribute
("qflag"), are replaced with NaN in the raw table. The raw DataFrame
retains the date index, but simplifies the columns to be only the measured
values for each data type. Two new columns are also created. The first is a
'yeardiff' column, which is only used for linear regression analysis as
detailed below. The second column, called SNPR, is calculated by the
equation (PRCP + SNOW/10) and is the rough estimate of the total liquid
precipitation on a particular day assuming that 1 inch of snow equals 0.1 inch
of liquid precipitation. The raw attribute is not saved to disk. If
previously downloaded data is reloaded in a new GHCND instance, the raw
DataFrame is instead recreated from the saved _raw_full DataFrame.
A variety of statistical analyses are then conducted on the filtered raw
data and saved into a stats DataFrame attribute. An example of this
DataFrame, as rendered in a Jupyter notebook, is shown below. Several of these
analyses are self-evident: lowest value=='min', highest value=='max',
average='mean', standard deviation='std'. However, the maximum and minimum
daily temperatures (TMAX and TMIN, respectively) are also subjected to a
linear regression analysis. In other words, these are the trends in the
highest and lowest temperatures for any given day in the year. The slope and
intercept of these trends are the 'slope' and 'icept' columns, respectively.
The 'p_slope' column is the p-value for the calculated slope. The two columns
'-log_p*slope' are the slopes weighted by the negative base-10 logarithm of
the p value. These values are convenient for ranking these slope values, as
will be described in the plotting section. This DataFrame is also saved to
disk using the store_type selected above.
These DataFrames attributes retain all of the methods for this object class as
well. For example, the statistics data can be saved as an Excel file using the
to_excel method.
temps.stats.to_excel('MSO_stats.xlsx')
Several plotting functions are also provided for visualizing the raw and
statistical weather data. Currently, these functions are only available for
Matplotlib and are found in the mplvis module, which can be imported
directly from pyghcnd.
from pyghcnd import GHCND, mplvis
A general call to the plotting functions is shown below -- additional positional and keyword arguments are discussed with the plotting examples below.
mplvis.general_function_call(ghcnd, ... , show=True, save=True, dpi=300)
The initial positional argument ghcnd must be a GHCND object instance
containing the data to be plotted. The final three keyword arguments --
show, save, and dpi -- control the plotting output. The show boolean
argument turns on/off visual display of the final figure using the
matplotlib.pyplot.show function. From the command line, this command will
open an interactive GUI window; in a Jupyter notebook, this will display the
plot after the code execution cell. The save boolean keyword argument turn on/off
saving of the final plot in the data_folder directory using the
matplotlib.pyplot.savefig function at the defined dpi value. The dpi
must be a positive integer. If both save and show are False, then no
output will be generated from the plotting function call.
Several examples of the plotting function output and argument descriptions are provided below. These examples are for the Missoula International Airport, and the weather data were updated in mid-March of 2023.
The plot_temp function is used to plot the statistical data as well as the
high/low temps for a given year, as shown above. The solid red and blue lines
are the 2023 high and low temps, respectively, and the solid black line is the
average daily temperature. The thin black vertical line indicates the current
date. The light blue/red shaded area shows the total range between the
largest and smallest of measured temperatures for that particular day. The
slightly darker shaded region shows the range defined by the average temp +/-
one standard deviation. I.e. 68% of the measured values fall in that range.
The dotted and dashed lines are the average temperatures estimated from the
linear regression analysis of the daily temperatures. For example, in the plot
shown, both the low and high daily temperatures have increased considerably
from 1948 (dotted line) to the present (dashed line) during January, March,
and the summer/early fall months. (But from this plot, it is not possible to
infer statistical significance from these changes.)
The call signature for this function is shown below.
mplvis.plot_temp(ghcnd, use_year=None, trends=True, smooth=15,
show=True, save=True, dpi=300)
The use_year keyword argument is used to select a year for plotting the
high/low temps as an overlay on the statistical values. Its default of None
causes it to select the current year. The trends Boolean keyword argument
plots additional lines representing the high/low temperatures as extrapolated
from the linear regression equation fitted to the daily temperature values.
The smooth keyword argument is used to apply a simple
moving average smoothing function to some of the plotted data, which can
improve the appearance of the plots. The value can either be a positive
integer (default: 15) that represents the number of points to use in the
smoothing average, or this parameter can be None or False, which plots the
unsmoothed data.
The plot_temp_diffs function is used to plot the deviations of the most
recent high/low temperatures from the average values. The color coding and
plotting elements are equivalent to the Yearly Temperature
Plot above. In this case, the center point (y = 0)
represents the average temperature for that particular day. Red/blue lines
greater than 0 are days that are warmer than average, and vice versa for
points below the center of the y-axis. In the example plot, the temperatures
for most days appear to be about average; however, there are a few
particularly cold streaks, with the cold snap in late February being records
lows for the daily high temperatures.
The call signature of this function is shown below.
mplvis.plot_temp_diffs(ghcnd, ndays=90, show=True, save=True, dpi=300)
The ndays keyword argument is a positive integer that determines the number
of days to use in the plot. The default (90) means that the deviations for
the final 90 days of data will be plotted.
The plot_prcp function is used to plot various precipitation data, such as
the cumulative precipitation for the year, daily precipitation values, and a
distribution of yearly total precipitation. The figure consists of three
subplots. The upper left plot show the cumulative precipitation (in inches)
for the year. The solid black line is the average cumulative precipitation for
the year based on all the data, and the solid red line is the cumulative data
for the year defined by the use_year keyword argument, described below. The
thin, semi-transparent blue lines are cumulative precipitation for every year
for which enough data is present (see n_missing keyword argument below). The
lower left bar plot shows daily precipitation values. Black bars are the
average daily precipitation, and red bars are the measured precipitation for
the year determined by the use_year keyword argument. The darker blue bars
are the average daily precipitation plus one standard deviation, and the
lightest blue bars are the maximum precipitation for a given day. The
histogram on the right shows a distribution of total precipitation for all the
years of station data. The solid vertical black line shows the average total
precipitation for the year, and the solid red vertical line is the total
precipitation for the year defined by the use_year keyword argument.
The call signature for this function is as shown below.
mplvis.plot_prcp(ghcnd, use_year=None, ptype='PRCP', n_missing=10,
show=True, save=True, dpi=300)
The use_year keyword argument a positive integer used to set the main
display year. If this argument is left with its default value of None, the
current year will be used. The ptype keyword argument is a string used to
select the type of data to plot. The acceptable values are 'PRCP'
(precipitation, default), 'SNOW' (snow), and 'SNPR' (sum of precipitation
and snow/10). The example shown above is using the 'SNPR' data. The
n_missing keyword is a positive integer that is used to filter out yearly
data for which the number of missing data are equal to or greater than this
value.
The plot_daily_temp function is used to plot high and low temperature trends
for a particular day of the year. The blue and red lines in the upper plots
are the measured low and high temperatures, respectively, for every year in
the stations records. The black line is a linear regression fit line for these
data. This line will either be a solid or dashed line, depending on the
p-value of the regression slope. If the p-value is greater than the p
keyword argument value, then the slope is not significant, and the trend line
will be dashed. Alternatively, statistically significant slopes will be
plotted as solid lines. For the example plot, it appears that these
temperatures do not show a statistically significant change over the lifetime
of the station data. A method for plotting the most significant trends will be
discussed below. The lower plots show the residuals, or trend minus actual
temperatures, for the low and high temperature values. Ideally, the residuals
should be randomly distributed about zero.
The call signature for this function is shown below.
mplvis.plot_daily_temp(ghcnd, month, day, temp_type='both', p=0.05,
show=True, save=True, dpi=300):
The two positional arguments month and day are positive integers that
correspond to the month (1-12) and day (1-31) of the year for which to plot
the trends. The temp_type keyword argument is a string used to select the
type of data to plot. The default 'both' creates a plot for both the high
and low daily temps. This can also be 'TMIN' or 'TMAX', if a plot of only
the low or high temperature data, respectively, is desired. The p keyword
argument is a floating point number that corresponds to the p-value cutoff for
determining statistical significance. The default is 0.05, which corresponds
to a 95% confidence level.
The plot_daily_trends function is used to plot all of the daily trend
information for every day on a single plot, which allows for comparisons of all
trends throughout the year. The left-hand
scatter plots show the daily trend slopes by date for the daily low and high
temps. Statistically significant trends, as defined by p-value<p keyword
argument value, are highlighted in black. The two plots on the right are
showing the same data but in histogram form. In this case, a considerable
number of days appear to be getting warmer in both the day and night.
The call signature for this function is shown below.
mplvis.plot_daily_trends(ghcnd, p=0.05, show=True, save=True, dpi=300)
The p keyword argument is a floating point number
that corresponds to a p-value cutoff for highlighting statistically
significant trends, as described below. The default 0.05 corresponds to a
95% confidence level cutoff.
The slope values can be sorted directly; however, this does not take into
account the statistical significance of this fit parameter. Two additional
columns of data -- -log_p*slope and -log_p*abs_slope -- have been provided
to sort trends based on the p-value of the slope. The weighting in this case
is based on the -log10(p-value), which will be a larger positive value for
smaller (i.e. more significant) p-values. These weighting values are then
either multiplied by the slopes, -log_p*slope, which retains the sign of the
data. These can be used to find either the most positive or most negative
slope. The weights are also multiplied by the absolute value of the of the
slope data, -log_p*abs_slope, this can be used to find the most significant
slopes regardless of whether the temperatures are increasing or decreasing.
Sorting the daily temperature trends can be done using the builtin DataFrame sorting methods, as shown in Jupyter notebook code cell below.
Using a simple loop and the plot_daily_temp function, this sorting code can
be used to plot the largest daily trend for both the high and low
temperatures, as shown below.
for cat in ('TMIN', 'TMAX'):
slope_sorted = temps.stats[cat]\
.sort_values(by='-log_p*abs_slope', ascending=False)
first_slope = slope_sorted.iloc[0]
month, day = first_slope.name
mplvis.plot_daily_temp(temps, month=month, day=day, temp_type=cat)
