piecewise-fit-gat-series.R

### For plots at https://twitter.com/priscian/status/1184499661656219651

## Install package "climeseries":
Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS = "true") # https://github.com/r-lib/remotes#environment-variables
devtools::install_github("priscian/climeseries")

## First plot: piecewise linear fit of major GAT series at significant changepoints
library(climeseries)
d <- get_climate_data(download = FALSE, baseline = TRUE)
series <- c("GISTEMP Global", "NCEI Global", "HadCRUT4 Global", "Cowtan & Way Krig. Global", "BEST Global (Air Ice Temp.)", "JMA Global", "RSS TLT 4.0 -70.0/82.5", "UAH TLT 6.0 Global", "ERA-Interim 2m Global", "ERA5 2m Global", "NCEP Surface Air Global")
offset <- 0.0
for (s in series) {
  d[[s]] <- d[[s]] + offset
  offset <- offset + 0.25
}
plot_climate_data(d, series = series, 1970, yearly = TRUE, lwd = 1, segmented = TRUE, plot.segmented... = list(lwd = 1), mark_segments = "points", points.segmented... = list(lwd = 2), ylim = NULL, save_png = FALSE)

## Second plot: piecewise linear fit of major GAT series, with influence of exogenous factors removed, at significant changepoints
## N.B. Give this one a couple of minutes to run!
library(climeseries)
series <- c("GISTEMP Global", "NCEI Global", "HadCRUT4 Global", "Cowtan & Way Krig. Global", "BEST Global (Air Ice Temp.)", "JMA Global", "RSS TLT 4.0 -70.0/82.5", "UAH TLT 6.0 Global", "ERA-Interim 2m Global", "ERA5 2m Global", "NCEP Surface Air Global")
start <- 1970; end <- NULL
g <- recenter_anomalies(remove_exogenous_influences(series = series, start = start, end = end, max_lag = 12), baseline = TRUE)
h <- g
series_adj <- paste(series, "(adj.)")
offset <- 0.0
for (s in series_adj) {
  h[[s]] <- h[[s]] + offset
  offset <- offset + 0.25
}
main <- "Adjusted for ENSO, Volcanic, and Solar Influences"
plot_climate_data(h, series = series_adj, end = NULL, ma = 12, yearly = TRUE, main = main, lwd = 1, baseline = FALSE, segmented = TRUE, plot.segmented... = list(lwd = 1), mark_segments = "points", points.segmented... = list(lwd = 2), trend = FALSE, ylim = NULL, save_png = FALSE)