Home / Tutorials

Weather Data Analysis in R

This tutorial covers the visualization and basic descriptive analysis of historic weather data in R.

R Libraries

The following R libraries will be used by code in this tutorial and you should install them if you haven't already.

install.packages("xts")
install.packages("tsbox")
install.packages("forecast")

Acquire the Data

The National Oceanic and Atmospheric Administration (NOAA) is a US federal government agency formed within the Commerce Department in 1970 by combining the U.S. Coast and Geodetic Survey (founded under President Jefferson in 1807 as the Survey of the Coast), the Weather Bureau (founded in 1870), and the US Commission of Fish and Fisheries (founded in 1871) (NOAA 2021).

NOAA's mission is, " To understand and predict changes in climate, weather, oceans, and coasts, to share that knowledge and information with others, and to conserve and manage coastal and marine ecosystems and resources" (NOAA 2021).

As part of that mission, NOAA makes a wide variety of earth observation data available to the public, including historical weather observations collected by NOAA's National Weather Service.

The video below shows how to download historical weather observation data from NOAA's Climate Data Online web portal.

1. Under Select Weather Observation Type/Dataset, select Daily Summaries.
2. Leave Select Date Range as the default because the available data will depend on which weather observation station you select on the next screen.
3. Under Search For, leave the default of Stations.
4. Under Enter a Search Term, enter the name of your county and state.
5. When the map pops up, select More Search Options, select Air Temperature and Precipitation as well as a 20-year date range to the present.
6. Run the search and click on station icons on the map to find a station that a Period of Record with adequate date coverage. Twenty years or more will be most useful for trend analysis. Airports are usually a safe bet.
7. In the selected station, click, ADD TO CART.
8. Click on the Cart in the top right corner
9. On the Select Cart Options page, under Select the Output Format, choose Custom GHCN-Daily CSV.
10. Under Select the Date Range, select the desired date range. Note that this dialog will limit you to the dates for which data is available for this particular station and will also limit you to 25 years. For trend analysis, you will likely want to select all available dates.
11. Click CONTINUE to proceed to the Custom Options: Daily Summaries page.
12. Under Station Detail & Data Flag Options, select Station Name and Geographic Location.
13. Under Select data types for custom output choose the selected options. For this example, we use all Precipitation, Air Temperature, and Wind variables and CONTINUE.
14. On the Review Order page, Enter e-mail address and SUBMIT ORDER.
15. Wait for an e-mail with a link to your data.
16. Download the linked CSV to your local drive.

Process the Data

Load the Data

Note that the default name for this file will be your order number and you should probably rename it to something more meaningful so you can keep track of which file is which.

You can load the .csv data file with the read.csv() function and list the available variables with the names() function.

csv = read.csv("2650569-daily-summary.csv", as.is=T)

names(csv)

	 [1] "STATION"   "NAME"      "LATITUDE"  "LONGITUDE" "ELEVATION" "DATE"     
	 [7] "AWND"      "FMTM"      "PGTM"      "PRCP"      "SNOW"      "SNWD"     
	[13] "TAVG"      "TMAX"      "TMIN"      "WDF2"      "WDF5"      "WSF2"     
	[19] "WSF5"

range(csv$DATE) 

	[1] "1997-06-01" "2021-07-10"

The variables of particular interest for this analysis are the following.

• TMAX: Maximum temperature in degrees Fahrenheit
• TMIN: Minimum temperature in degrees Fahrenheit
• PRCP: Precipitiation in inches

Convert to a Time Series

A time series is "a set of data collected sequentially usually at fixed intervals of time" (Merriam-Webster 2021). Time series data has characteristics that make it very useful for analyzing phenomenon that change over time.

There are a wide variety of different libraries and classes are available for manipulating and analyzing time series data in R, and deciding which one to you can be confusing (Hyndman 2021).

This tutorial will use the xts library, which provides an extension of the zoo and ts classes for representing time series data in R.

This tutorial will also use the tsbox library for conversion between time series representations in order to take advantage of capabilities in the different libraries.

• The xts() function creates a time series object. The order.by parameter uses the dates from the data, which are converted to R Date objects by the as.Date() function.
• The ts_regular() function gives the time series a regular (daily) interval by adding NA values for missing dates. A regular interval will be needed later for decomposition analysis.
• The na.fill() function fills in those missing dates by duplicitating (extending) values from previous days.
• The window() function clips off the starting and ending dates so the number of years covered is a multiple of four. This will be needed later when the data needs to be aggregated into monthly periods.

library(xts)
library(tsbox)

historical = xts(csv[,c("TMAX","TMIN","PRCP")], order.by=as.Date(csv$DATE))

historical = ts_regular(historical)

historical = na.fill(historical, "extend")

historical = window(historical, start=as.Date("2000-01-01"), end=as.Date("2020-12-31"))

A plot() of TMIN and TMAX show the annual cycle of temperatures as well as extreme temperature events that spike above the general curve.

The ts_ts() function is used to plot ts objects rather than xts objects so that more control is available over the formatting of the charts.

plot(ts_ts(historical$TMAX), col="darkred", bty="n", las=1, fg=NA, 
	ylim=c(-20, 120), ylab="Temperature (F)")

lines(ts_ts(historical$TMIN), col="navy")

grid(nx=NA, ny=NULL, lty=1, col="gray")

legend("topright", fill=c("darkred", "navy"), cex=0.7,
	legend=c("TMAX", "TMIN"), bg="white")

The plot() of daily precipitation shows no clear seasonal pattern, although the presence of a limited number of high precipitation days over the 23 years of this data is notable.

barplot(historical$PRCP, border=NA, col="darkgreen", ylim=c(0, 2),
	space=0, bty="n", las=1, fg=NA, ylab="Daily Rainfall (inches)")

grid(nx=NA, ny=NULL, lty=1)

Analyze the Data

Summary Descriptive Statistics

Running the summary() function will give basic descriptive statistics for the time series fields. You can subset the time series with logical conditions to find the dates for extreme conditions.

summary(historical)

	     Index                 TMAX             TMIN            PRCP     
	 Min.   :2000-01-01   Min.   : -4.00   Min.   :-18.0   Min.   :0.00  
	 1st Qu.:2005-04-01   1st Qu.: 45.00   1st Qu.: 28.0   1st Qu.:0.00  
	 Median :2010-07-02   Median : 66.00   Median : 43.0   Median :0.00  
	 Mean   :2010-07-02   Mean   : 62.69   Mean   : 41.9   Mean   :0.10  
	 3rd Qu.:2015-10-01   3rd Qu.: 81.00   3rd Qu.: 58.0   3rd Qu.:0.03  
	 Max.   :2020-12-31   Max.   :102.00   Max.   : 78.0   Max.   :3.76  

historical[historical$TMIN == min(historical$TMIN)]

		   TMAX TMIN PRCP
	2009-01-16    9  -18    0

historical[historical$TMAX == max(historical$TMAX)]

		   TMAX TMIN PRCP
	2001-07-20  102   55 0.02

historical[historical$PRCP == max(historical$PRCP)]

		   TMAX TMIN PRCP
	2014-05-21   91   66 3.76

Seasonal Statistics

Time series representation also facilitates summary statistics of recurring periods, like months.

• The format() function with the %m format descriptor isolates the month portion of the date.
• The split() function breaks a vector into a data frame with separate fields based on a specific factor, in this case the month. The as.numeric() conversion is used to remove the time series date index so that dates are not included in the summaries.
• The sapply() function runs the summary() function on each field in the split data frame.

historical$MONTH = format(index(historical), "%m")

months = split(as.numeric(historical$TMAX), historical$MONTH)

sapply(months, summary)

		       1        2        3        4       5        6        7        8
	Min.    -4.00000  4.00000 12.00000 33.00000 48.0000 56.00000  65.0000 63.00000
	1st Qu. 27.00000 30.00000 42.00000 57.00000 68.0000 79.00000  81.0000 80.00000
	Median  34.00000 36.00000 50.00000 65.00000 77.0000 83.00000  85.0000 83.00000
	Mean    34.26575 37.79798 51.47005 64.94524 75.3149 83.23571  84.5169 83.37327
	3rd Qu. 41.00000 45.75000 61.00000 74.00000 84.0000 88.00000  88.0000 87.00000
	Max.    67.00000 72.00000 86.00000 90.00000 97.0000 99.00000 102.0000 97.00000
		       9       10       11       12
	Min.    52.00000 40.00000 16.00000  2.00000
	1st Qu. 74.00000 57.00000 43.00000 31.00000
	Median  80.00000 66.00000 52.00000 38.00000
	Mean    79.56111 66.24424 51.85079 38.42089
	3rd Qu. 85.00000 75.00000 59.75000 47.00000
	Max.    99.00000 93.00000 81.00000 73.00000

months = split(as.numeric(historical$TMIN), historical$MONTH)

sapply(months, summary)

			1         2        3        4        5       6        7
	Min.    -18.00000 -14.00000 -9.00000 12.00000 -1.00000 -1.0000 -1.00000
	1st Qu.   8.00000  13.00000 25.00000 33.00000 46.00000 57.0000 60.00000
	Median   20.00000  21.50000 31.00000 41.00000 52.00000 62.0000 65.00000
	Mean     17.64286  20.26599 31.00998 40.89603 52.34562 61.3254 63.74731
	3rd Qu.  28.00000  29.00000 37.00000 48.00000 60.00000 66.5000 69.00000
	Max.     57.00000  52.00000 61.00000 65.00000 72.00000 74.0000 78.00000
		       8        9       10       11        12
	Min.    -1.00000 -1.00000 -1.00000 -1.00000 -10.00000
	1st Qu. 58.00000 49.00000 36.00000 26.00000  15.00000
	Median  63.00000 55.00000 42.00000 32.00000  24.00000
	Mean    61.76805 54.63095 42.88633 32.19048  22.87404
	3rd Qu. 66.00000 61.00000 49.00000 39.00000  31.00000
	Max.    75.00000 72.00000 68.00000 61.00000  58.00000

Time Series Decomposition

Time series decomposition separates a time series into three fundamental components that can be added together to create the original data:

• A seasonal component
• A long-term trend component
• A random component

When using historical data to anticipate what might happen in the future, we often wish to separates the seasonal and random components to see the long term historical trends.

The stl() function performs this decomposition.

• The ts_ts() function from the library converts an xts field to a ts object that can be used with stl().
• The s.window parameter defines the smoothing to isolate the seasonal component (365 days = one year).
• The t.window parameter defines the smoothing to isolate the long-term trend (~20 years).

decomposition = stl(ts_ts(historical$TMAX), s.window=365, t.window=7001)

plot(decomposition)

summary(decomposition$time.series[,"trend"])

	   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
	  62.33   62.43   62.63   62.72   63.00   63.36 

In this example, the trend indicates an increase of around one degree (F) in daily maximum temperatures over the period 2000 - 2020. In daily life, this is masked by seasonal swings (summer to winter) of around 45 degrees, and random daily variation of around 60 degrees.

Performing the same analysis on daily precipitation, we also see an small trend in daily increase in precipitation of around 0.06 inches over the past 23 years, which translates to around one additional inch of rain per year (0.06 * 365 / 23).

decomposition = stl(ts_ts(historical$PRCP), s.window=365, t.window=7001)

plot(decomposition)

summary(decomposition$time.series[,"trend"])

Aggregation

Although having daily weather observations is extremely useful for analysis, direct visualization of daily observations relies on inexact (and occasionally deceptive) visual identification of patterns and trends. Visualization of decomposition captures trends, but can be confusing to folks with little or no exprience with time series analysis.

One common technique to address this issue with time series data is to aggregate values by months or years.

• With xts time series objects, the period.apply() function can be used to perform operations on blocks of observations across the time series. Those operations are specified with the FUN= parameter and commonly include min(), max(), and mean().
• period.apply() requires a regular (even) interval between periods. Because months are uneven, a seq() of regularly spaced index values is passed to the INDEX= parameter. The INDEX value is is a regular spacing that approximates 12 months per year when leap years are considered (365 + 365 + 365 + 366 / 48 = 30.4375).
• Because the periods need an even number of observations, the window() dates should be chosen so they encompass even multiples of four years to account for leap years.

monthly.tmax = period.apply(historical$TMAX, INDEX = seq(1, nrow(historical) - 1, 30.4375), FUN = mean)

plot(ts_ts(monthly.tmax), col="darkred", ylim=c(20, 100), 
	lwd=3, bty="n", las=1, fg=NA, ylab="TMAX (F)")

grid(nx=NA, ny=NULL, lty=1)

The aggregation by mean() smooths out the daily randomness to make the overall cycles and trends clearer.