Global Temperature

Global land and ocean temperature anomaly from NASA GISS in degrees Celcius.

temp <- readRDS("giss.globe.rds")
summary(temp)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.7800 -0.2300 -0.0500  0.0231  0.2250  1.3300
tsp(temp)
## [1] 1880.000 2016.833   12.000

Moving Averages

Create filters and plot three moving averages. Filters include: three years (36 months), ten years (120 months), and twenty years (240 months).

# Moving average of 3 year
ma03 <- filter(x=temp, filter=rep(x=1/36, times=60), sides=2)

# Moving average of 10 years
ma10 <- filter(x=temp, filter=rep(x=1/120, times=120), sides=2)

# Moving average of 20 years
ma20 <- filter(x=temp, filter=rep(x=1/240, times=240), sides=2)

plot(temp, col="grey", ylab= "Global Temperature (C)")
lines(ma03, col="deeppink", lwd=2)
lines(ma10, col="green", lwd=2)
lines(ma20, col="blue", lwd=2)

Low frequencies retain the most value (three years), while higher frequency (twenty years) looses some of the values on the ends. We get a smoother line, but at the expense of lost data at the ends. All moving averages show a general upward trend in global temperatures from 1880 to 2016.

Hanning Filter

The Hanning filter is similar to the moving average by highlighting low frequency variability, but is smoother than the moving average. This filter also loses the ends of the dataset, because it doesn’t have enough information at the ends to calculate averages.

library(dplR)
## Warning: package 'dplR' was built under R version 3.3.2
han03 <- hanning(temp, n=36)
han10 <- hanning(temp, n=120)
han20 <- hanning(temp, n=240)

# Turn each hanning function object into a time series with the same properties as the temp time series, so they will plot properly together
han03.ts <- ts(han03, start=c(1880), frequency = 12)
han10.ts <- ts(han10, start=c(1880), frequency = 12)
han20.ts <- ts(han20, start=c(1880), frequency = 12)

par(tcl=0.5, mar=rep(3,4), mgp=c(1.1, 0.1, 0), xaxs="i")
plot(temp, col="grey", ylab= "Global Temperature (C)")
abline(h=-0.05)
lines(han03.ts, col="deeppink", lwd=2)
lines(han10.ts, col="green", lwd=2)
lines(han20.ts, col="blue", lwd=2)

Smoothing Splines

The smoothing spline function fits a smooth curve to the dataset using piecewise polynomials. The term spline comes from a physical object used in architectural drawings to bend and run a pencil along. The spline runs across the entire dataset, so in this method we don’t lose data at the ends. The spar parametric adjusts the “smoothness” of the data. Larger numbers smooth the data more than smaller numbers.

Below I plot three different “spar” values to adjust the smoothness of the temperature data.

spl.2 <- smooth.spline(temp,spar=0.2)
spl.4 <- smooth.spline(temp,spar=0.4)
spl.8 <- smooth.spline(temp,spar=0.8)

par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(temp,type="l",xlab="Year",ylab="Temp (C)",col="grey")
abline(h=-0.05)
lines(spl.2, col = "deeppink", lwd = 2)
lines(spl.4, col = "green", lwd = 2)
lines(spl.8, col = "blue", lwd = 2)
axis(3);axis(4)

The nyrs argument fits a spline that many “years” to the data.

spl32 <- ffcsaps(temp,nyrs=32)
spl64 <- ffcsaps(temp,nyrs=64)
spl128 <- ffcsaps(temp,nyrs=128)

# Make each spline object have time series properties
spl32.ts <- ts(spl32, start=c(1880), frequency = 12)
spl64.ts <- ts(spl64, start=c(1880), frequency = 12)
spl128.ts <- ts(spl128, start=c(1880), frequency = 12)

par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(temp,type="l",xlab="Year",ylab="Temp (C)",col="grey")
abline(h=-0.05)
lines(spl32.ts, col = "deeppink", lwd = 2)
lines(spl64.ts, col = "green", lwd = 2)
lines(spl128.ts, col = "blue", lwd = 2)
axis(3);axis(4)

Lowess

The lowess uses a local linear polynomial fit to smooth the data. Adjusting the span f, you change the proportion of points used in the smoothing, which means that a smaller f value makes a stiffer curve fit. Larger f values make a smoother curve.

I start with a 36 months (3 year) span and plot this filter. Then, I take the residuals from the filter to remove low-frequency variability.

# Make an object with the length of the temp dataset
n <- length(temp)

# Make an object containing 1 through the lengh, which represents the number of years.
yrs <- 1:n

# 36 months divided by the total number of observations in the dataset to be used in the lowess filter
f36 <- 36/n
f36.lo <- lowess(temp, f= f36)

# Get the residuals of the lowess filter by dividing the temp time series by the values of the lowess filter
residual.36 <- temp/f36.lo$y

par(mfcol=c(2,1),mar=c(0.1,3,3,3),mgp=c(1.2,0.2,0),tcl=0.5)

plot(temp,type="l",xlab="Year",ylab="Temp (C)",col="grey", xaxt='n')
abline(h=0.05)
lines(f36.lo$x, f36.lo$y, col = "deeppink", lwd = 2)
axis(3)

# Plot the residuals below
par(mar=c(3,3,0.1,3))
plot(f36.lo$x, residual.36, type="l", ylab='Residuals', xlab='Years', col='deeppink')

Now, I’ll look at a 120 month (10 year) filter.

n <- length(temp)
yrs <- 1:n

f120 <- 120/n
f120.lo <- lowess(temp, f= f120)
residual.120 <- temp/f120.lo$y

par(mfcol=c(2,1),mar=c(0.1,3,3,3),mgp=c(1.2,0.2,0),tcl=0.5)

plot(temp,type="l",xlab="Year",ylab="Temp (C)",col="grey", xaxt='n')
abline(h=0.05)
lines(f120.lo$x, f120.lo$y, col = "green", lwd = 2)
axis(3)

par(mar=c(3,3,0.1,3))
plot(f120.lo$x, residual.120, type="l", ylab='Residuals', xlab='Years', col='green')

Aggregating the Sea Ice Data

This data shows the Northern hemisphere sea ice extent in million km^2 from 1979 to 2016. The moving average for 1 year and 2 years isn’t as interesting, because it has such a strong seasonal cycle.

We can use a filter to aggregate the sea ice dataset to look at the average sea ice extent for the winter months of January and February.

# Read in sea ice data
seaice <- readRDS("ice.nx.rds")
summary(seaice)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.630   9.218  12.380  11.750  14.540  16.520
# Moving average of 12 months or 1 year (ma01)
ma01 <- filter(x=seaice, filter=rep(x=1/12, times=12), sides=2)

# Moving average of 2 year (ma02)
ma02 <- filter(x=seaice, filter=rep(x=1/24, times=24), sides=2)

plot(seaice, col="grey", ylab= "Sea Ice Extent (million Km^2)")
lines(ma01, col="blue", lwd=2)
lines(ma02, col="pink", lwd=2)

yrs <- 1979:2016
# Make a 2-month filter
ice.2.filt <- filter(x= seaice, filter = rep(x=1/2, times=2), sides=2)

winter.index <- seq(from=1 , to= length(ice.2.filt), by=12)
ice.winter <- ice.2.filt[winter.index]

# Note length
length(ice.winter)
## [1] 38
# Compare
mean(seaice[1:2])
## [1] 15.99
ice.winter[1]
## [1] 15.99
par(mar=rep(4,4))
plot(yrs, ice.winter, type="b", xlab="Year", ylab= expression(paste("Sea Ice Extent (", 1e6~km^2 , ")")))
abline(lm(ice.winter ~ yrs), col="blue")

lm(ice.winter~yrs)
## 
## Call:
## lm(formula = ice.winter ~ yrs)
## 
## Coefficients:
## (Intercept)          yrs  
##   106.05373     -0.04561

The linear model and plot indicates that sea ice extent in Jan and Feb months from 1979 to 2016 has decreased at a rate of about 50,000 \(km^2\) per year.

Now aggregate for the temperature data

Use filter to aggregate the global temperature dataset to look at the average temperature (C) for the summer months of June, July, and August. This enables us to see a trend in one part of the year.

yrs <- 1880:2016

# 3 month moving average filter
ma3 <- filter(x=temp, filter=rep(x=1/3, times=3), sides=2)
# Note the NA
head(cbind(temp, ma3))
##       temp        ma3
## [1,] -0.30         NA
## [2,] -0.21 -0.2300000
## [3,] -0.18 -0.2200000
## [4,] -0.27 -0.1966667
## [5,] -0.14 -0.2333333
## [6,] -0.29 -0.2233333
plot(temp, col="grey", ylab= "Global Temperature (C)")
lines(ma3, col="red", lwd=2)

summer.index <- seq(from=7, to=length(ma3), by=12)
temp.summer <- ma3[summer.index]

# Check the length of the temp.summer
length(temp.summer)
## [1] 137
# Compare
mean(temp[6:8])
## [1] -0.2
temp.summer[1]
## [1] -0.2
par(mar=rep(4,4))
plot(yrs, temp.summer, type="b", xlab="Year", ylab= expression(degree*C))
abline(lm(temp.summer~yrs), col="red")

Spectral Analysis on Milankovitch Cycles

the Milankovitch cycles describe the changes in the way the Earth orbits the sun and how this influences the Earth’s climate over thousands of years. These changes mark the sequence of ice ages and warm periods in Earth’s history.

The Earth’s orbit varies from nearly circular and slightly eliptical (eccentricity), which influences the amount of solar radiation at different times of the year and overall climate patterns.

The dataset (jul65N.rds) looks at the July insolation at 65N in thousands of years ago.

jul65N <- readRDS("jul65N.rds")

plot(jul65N$kya,jul65N$W.per.m2,type="l",
     xlab="Thousands of years ago",
     ylab=expression(W/m^2),main="July insolation at 65N")

# Make it a time series
jul.ts <- ts(jul65N$W.per.m2, start = -5000, end=0 , frequency=1)
spectrum(jul.ts, log="n")

Let’s take a closer look at the interesting part of the periodogram with a span of 4.

jul.spec <- spectrum(jul.ts, span=4, plot=FALSE)
plot(jul.spec, log="n", type="h", xlab= expression(paste("Frequency of July Solar Irradiance (W ", m^2, ") in thousands of years")), ylab = "Spectral density", main="Milankovitch Cycles", xlim= c(0, 0.06))

Period = 1 / Frequency

The strong peak in fequency just above 0.04 indicates there is a dominant periodicy in the data every 25 thousand years. The slightly smaller peak in frequency just above 0.05, indicates that there is also a cycle evident around the 20 million year cycle. The much smaller peak in fequency at about 0.025 indicates a less apparent low frequency cycle every 40-41 thousand year.

These data indicate that there are three dominant Milankovitch cycles when the tilt of the Earth’s axis changes and influences an increase or decrease in the amount of solar irriadiance (W\(m^2\)) at 65N.

Now, we can plot these dominant periodicies in the data using the Lowess Filter to get a better visual of these periodicies.

# 
n <- length(jul.ts)
yrs <- 1:n

# Make a filter of 20 thousand years
f20 <- 20/n
f20.lo <- lowess(jul.ts, f= f20)

# Make a filter of 25 thousand years
f25 <- 25/n
f25.lo <- lowess(jul.ts, f=f25)

# Make a filter of 40 thousand years
f40 <- 40/n
f40.lo <- lowess(jul.ts, f=f40)

plot(jul.ts,type="l", main= "20 Thousand Year Filter", xlab="Thousands of years", ylab= expression(paste("Summer solar irradiance (W", m^2, ")" )), col="grey")
abline(h=0.05)
lines(f20.lo$x, f20.lo$y, col = "chocolate", lwd = 2)

plot(jul.ts,type="l", main= "25 Thousand Year Filter", xlab="Thousands of years", ylab= expression(paste("Summer solar irradiance (W", m^2, ")" )), col="grey")
abline(h=0.05)
lines(f25.lo$x, f25.lo$y, col = "red", lwd = 2)

plot(jul.ts,type="l", main= "40 Thousand Year Filter", xlab="Thousands of years", ylab= expression(paste("Summer solar irradiance (W", m^2, ")" )), col="grey")
abline(h=0.05)
lines(f40.lo$x, f40.lo$y, col = "orange", lwd = 2)

# 
n <- length(jul.ts)
yrs <- 1:n

# Make a filter of 20 thousand years
f20 <- 20/n
f20.lo <- lowess(jul.ts, f= f20)

# Make a filter of 40 thousand years
f40 <- 40/n
f40.lo <- lowess(jul.ts, f=f40)

plot(jul.ts,type="l", xlab="Thousands of years", ylab= expression(paste("Summer solar irradiance (W", m^2, ")" )), col="grey")
abline(h=0.05)
lines(f25.lo$x, f25.lo$y, col = "red", lwd = 2)
lines(f40.lo$x, f40.lo$y, col = "orange", lwd = 2)

References

Berger A. and Loutre M. F. (1991) Insolation values for the climate of the last 10 million years. Quaternary Sciences Review, Vol. 10 No. 4, pp. 297-317, 1991.