Plots the Ljungqvist 2010 temperature reconstruction series along with the instrumental HadCRUT4 30N-90N (120-mo. moving average) temp series.

Ljungqvist 2010 Reconstruction+HadCRUT4_0004.6-2017.9_ma120_baseline1961-1990.R

library(climeseries) # devtools::install_github("priscian/climeseries")
library(xlsx)

ljungqvist <- xlsx::read.xlsx2(system.file("extdata/paleo/ljungqvist2010.xls", package="climeseries"), sheetName="Reconstruction", startRow=11, check.names=FALSE, stringsAsFactors=FALSE)
matches <- str_match(ljungqvist$Decade, "(\\d+)\u2013(\\d+)")
ljungqvist$yr_part <- apply(matches[, 2:3], 1, function(x) mean(as.numeric(x)))

yearRange <- range(as.numeric(matches[, 2:3]))
allYears <- seq(yearRange[1], yearRange[2])

e <- get_climate_data(download=FALSE, baseline=FALSE)
e <- base::merge(expand.grid(month=1:12, year=allYears), e, by=c("year", "month"), all=TRUE)
e$yr_part <- e$year + (2 * e$month - 1)/24

e <- tbl_dt(e)
l <- tbl_dt(dataframe(
  yr_part = ljungqvist$yr_part,
  `Ljungqvist 2010 Reconstruction` = as.numeric(ljungqvist$`Reconstructed temperature`),
  `Ljungqvist 2010 Reconstruction_uncertainty` = as.numeric(ljungqvist$`Upper error bar`) - as.numeric(ljungqvist$`Lower error bar`)))
data.table::setkey(e, yr_part)
g <- dplyr::left_join(e, e[, .(year, month, yr_part), with=TRUE][l, roll="nearest"][, yr_part := NULL], by=c("year", "month"))

## HadCRUT4: https://crudata.uea.ac.uk/cru/data/temperature/HadCRUT.4.6.0.0.median.nc
dataDir <- getOption("climeseries_data_dir")
flit <- "HadCRUT.4.6.0.0.median.nc"
download.file("https://crudata.uea.ac.uk/cru/data/temperature/HadCRUT.4.6.0.0.median.nc", paste(dataDir, flit, sep="/"), mode="wb", quiet=TRUE)
n <- nc_open(paste(dataDir, flit, sep="/")) # 'print(n)' or just 'n' for details.
a <- ncvar_get(n, "temperature_anomaly")
## Structure of 'a' is temperature_anomaly[longitude, latitude, time], 72 × 36 × Inf (monthly since Jan. 1850)
lat <- ncvar_get(n, "latitude")
long <- ncvar_get(n, "longitude")
times <- ncvar_get(n, "time")
tunits <- ncatt_get(n,"time", "units")
nc_close(n)

tunits
# $value
# [1] "days since 1850-1-1 00:00:00"
dtimes <- as.Date(times, origin="1850-01-01")
## This data set should have the same length as the "time" dimension of 'a':
h <- data.frame(year=year(dtimes), month=month(dtimes), check.rows=FALSE, check.names=TRUE, fix.empty.names=TRUE, stringsAsFactors=FALSE)
# h$yr_part <- h$year + (h$month - 0.5)/12

subLat <- c(30, 90); subLong <- c(-180, 180)

flit <- apply(a, 3,
  function(y)
  {
    x <- t(y)
    w <- cos(matrix(rep(lat, ncol(x)), ncol=ncol(x), byrow=FALSE) * (pi / 180)) # Latitude weights.

    ## Use only subgrid for calculations.
    keepSubGrid <- list(lat=lat >= subLat[1] & lat <= subLat[2], long=long >= subLong[1] & long <= subLong[2])
    x1 <- x[keepSubGrid$lat, keepSubGrid$long, drop=FALSE]
    w1 <- w[keepSubGrid$lat, keepSubGrid$long, drop=FALSE]

    nlat <- length(lat[keepSubGrid$lat])
    temp <- NULL
    for (i in seq(1L, nrow(x1), by=nlat)) {
      xi <- data.matrix(x1[i:(i + nlat - 1L), ])
      tempi <- stats::weighted.mean(xi, w1, na.rm=TRUE)

      temp <- c(temp, tempi)
    }

    temp
  })

instSeries <- "HadCRUT4 30N-90N"
instSeriesMa <- instSeries %_% " (120-mo. MA)"
h[[instSeries]] <- flit
h[[instSeriesMa]] <- MA(h[[instSeries]], 120)
h <- tbl_dt(h)

g <- dplyr::full_join(g, h, by=c("year", "month"))
## The following is a hack, but it should be very close:
g[[instSeries %_% "_uncertainty"]] <- g$`HadCRUT4 NH_uncertainty`
g[[instSeriesMa %_% "_uncertainty"]] <- MA(g$`HadCRUT4 NH_uncertainty`, 120)

series <- c("Ljungqvist 2010 Reconstruction", "HadCRUT4 30N-90N (120-mo. MA)")
plot_climate_data(as.data.frame(g), series, conf_int=TRUE, get_x_axis_ticks...=list(min=-3000, by=100), col=c("green", "red"), baseline=1961:1990)

### Calculate difference in yearly means for deniers.

## See code for 'get_yearly_difference()' at https://github.com/priscian/climeseries/blob/master/R/helper.R#L930

series <- c("HadCRUT4 30N-90N (120-mo. MA)")
ytd <- get_yearly_difference(series, 2000, loess=FALSE, loess...=list(span=0.4), data=g)

# Difference in °C from 2000–2017
#                                [,1]
# HadCRUT4 30N-90N (120-mo. MA) 0.526