140+年的气候变化建模与预测

本文是学习Climate change: Modeling 140+ years of temperature data with tsibble and fable的笔记。

加载包

首先是加载本文所需要的一些包：

library(tsibble)
library(tidyverse)
library(magrittr)
library(lubridate)
library(patchwork)
library(fable)

如果你的fable包安装失败了，你可能可以参考这篇文章Mac OS X R error “ld: warning: directory not found for option” - Stack Overflow，在终端运行下面两行代码解决：

Shell

1 2	curl -O http://r.research.att.com/libs/gfortran-4.8.2-darwin13.tar.bz2 sudo tar fvxz gfortran-4.8.2-darwin13.tar.bz2 -C /

准备数据集

我已经把两份数据集都下载好了：

TN_STAID000011.txt
TX_STAID000011.txt

将温度变成tsibbles并添加日平均

# 整理最低温数据：
raw <- read_csv('TN_STAID000011.txt', skip = 20)
raw %<>% select(-STAID, -SOUID) %>% 
  rename(
    date = DATE,
    min_temp = TN,
    qual = Q_TN
  ) %>% 
  mutate(
    min_temp = replace(min_temp, min_temp <= -9999, NA),
    min_temp = replace(min_temp, qual == 1 | qual == 9, NA),
    min_temp = min_temp * 0.1, 
    date = ymd(date)
  ) %>% 
  select(-qual)

min_temp <- raw %>% as_tsibble(index = date)

# 整理最高温数据
raw <- read_csv('TX_STAID000011.txt', skip = 20) %>% 
  select(-STAID, -SOUID) %>% 
  rename(
    date = DATE, 
    max_temp = TX,
    qual = Q_TX
  ) %>% 
  mutate(
    max_temp = replace(max_temp, max_temp <= -9999, NA),
    max_temp = replace(max_temp, qual == 1 | qual == 9, NA),
    max_temp = max_temp * 0.1,
    date = ymd(date)
  ) %>% 
  select(-qual)

max_temp <- raw %>% as_tsibble(index = date)

krem <- left_join(max_temp, min_temp) %>% 
  mutate(month = month(date),
         mean_temp = (max_temp + min_temp)/2)

绘图：

krem %>% 
  ggplot(aes(x = date, y = mean_temp)) + 
  geom_line(alpha = 0.8) + 
  scale_x_date(
    '年份', 
    limits = c('2001-01-01', '2019-01-01') %>% ymd(),
    date_breaks = '5 years',
    date_labels = '%Y'
  ) + 
  scale_y_continuous(
    '平均气温˚C',
    limits = c(-20, 30),
    breaks = seq(-20, 30, 10)
  ) + 
  geom_smooth(size = 1.2, color = 'orange') + 
  geom_hline(yintercept = 0, color = 'dodgerblue', 
             linetype = 2, size = 1.5) + 
  labs(
    title = '奥地利地区克雷姆斯明斯特: 2001-2018每日平均气温',
    caption = '数据来源：https://www.ecad.eu/'
  )

这幅图展示了奥地利地区克雷姆斯明斯特: 2001-2018每日平均气温，从图中可以略微看出，气温似乎有上升趋势。

可视化不同时期的月份温度范围

krem_narrow <- gather(krem, key = 'type',
                      value = 'Temp', 
                      min_temp, max_temp, mean_temp) %>% 
  as_tibble()

krem_range <- krem_narrow %>% 
  group_by(type, month) %>% 
  summarise(
    mean = mean(Temp, na.rm = T),
    max = max(Temp, na.rm = T),
    min = min(Temp, na.rm = T)
  )
month <- c('一月', '二月', '三月', '四月',
           '五月', '六月', '七月', '八月',
           '九月', '十月', '十一月', '十二月')  

krem_range %>% 
  ggplot(aes(x = month, y = mean, color = type)) +
  geom_pointrange(
    data = subset(krem_range, type == "min_temp"), 
    aes(
      ymin = min,
      ymax = max
    )
  ) + 
  geom_pointrange(
    data = subset(krem_range, type == "max_temp"), 
    aes(
      ymin = min,
      ymax = max
    ),
    position = position_nudge(0.2)
  ) + 
  scale_x_continuous("月份", breaks = 1:12, labels = month) + 
  scale_y_continuous("温度˚C", 
                     breaks = seq(-25, 35, 5)) + 
  scale_color_discrete(
    breaks = c('max_temp', 'min_temp'),
    labels = c('每日最高气温', '每日最低气温')
  ) + 
  guides(color = guide_legend(title = '温度范围')) + 
  theme(legend.position = c(0.12, 0.88)) + 
  geom_line(data = subset(krem_range, type == 'min_temp')) + 
  geom_line(data = subset(krem_range, type == 'max_temp'), position = position_nudge(0.2)) + 
  labs(title = '奥地利地区克雷姆斯明斯特: 1876-2018气温范围', caption = '数据来源：https://www.ecad.eu/')

温度异常与平均值相比

气候学家经常将温度与基准线的偏差定义为温度的异常。在这里，我们绘制出年平均异常与100年的基准线（1876——1976）。我们可以注意到自20世纪70年代以来，温度呈上升趋势。

krem_narrow %<>% mutate(
  day = yday(date)
)

krem_range <- krem_narrow %>% 
  filter(type == 'mean_temp',
         date < as.Date('1976-01-01')) %>% 
  group_by(day) %>% 
  summarise(centurymean = mean(Temp, na.rm = T))

anomaly <- krem_narrow %>% 
  filter(type == "mean_temp") %>% 
  left_join(krem_range) %>% 
  mutate(anomaly = Temp - centurymean) %>% 
  select(date, anomaly) %>% 
  mutate(year = year(date)) %>% 
  as_tsibble(index = date) %>% 
  index_by(year) %>% 
  summarise(anom = mean(anomaly, na.rm = T))

anomaly %>% 
  ggplot(aes(x = year, y = anom)) + 
  geom_line(alpha = 0.8) + 
  scale_x_continuous('年份', breaks = seq(1880, 2020, 20)) + 
  scale_y_continuous('温度异常˚C',
                     limits = c(-2, 4),
                     breaks = seq(-2, 4, 0.4)) + 
  geom_smooth(size = 1.5, color = 'red') + 
  geom_hline(yintercept = 0, 
             color = 'dodgerblue', 
             linetype = 2, 
             size = 2) + 
  labs(title = '奥地利地区克雷姆斯明斯特：1876—2018温度异常情况',
       subtitle = '基准时期：1876—1976, 红线为Loess拟合趋势',
       caption = '数据来源：https://www.ecad.eu/')

可以看出1970年以来温度异常确实有上升的趋势。

每年低于零度的天数

krem_narrow %<>% mutate(
  year = year(date),
  winter = ifelse(month > 6, 1, 0),
  year = year + winter
) %>% 
  select(-winter)

krem_high <- krem_narrow %>% 
  filter(type == 'max_temp') %>% 
  select(date, Temp, year)

krem_narrow %<>% filter(type == 'mean_temp') %>% 
  select(date, Temp, year)

num_freeze_days <- krem_narrow %>% 
  group_by(year) %>% 
  filter(Temp <= 0) %>% 
  summarise(num_freeze_days = n()) %>% 
  filter(year <= 2018) %>% 
  as_tsibble(index = year) %>% 
  fill_na()

ggplot(num_freeze_days, 
       aes(x = year, y = num_freeze_days)) + 
  geom_line(alpha = 0.8) + 
  scale_x_continuous('年份', 
                     breaks = seq(1880, 2020, 20)) + 
  scale_y_continuous('结冰的天数',
                     limits = c(0, 110),
                     breaks = seq(0, 110, 20)) + 
  geom_smooth(size = 1.5, color = 'dodgerblue') + 
  labs(title = '奥地利地区克雷姆斯明斯特：1876-2017年结冰的天数',
       subtitle = '每年气温低于0度的天数，蓝线为Loess拟合',
       caption = '数据来源：https://www.ecad.eu/')

由于气候变暖，每年气温低于0度的天数越来越少了。

每年气温高于30˚C的天数

num_hot_days <- krem_high %>% 
  group_by(year) %>% 
  filter(Temp >= 30) %>% 
  summarise(num_hot = n()) %>% 
  as_tsibble(index = year) %>% 
  fill_gaps() %>% 
  replace_na(list(num_hot = 0))

num_hot_days %>% 
  ggplot(aes(x = year,
             y = num_hot)) + 
  geom_line(alpha = 0.7) + 
  scale_x_continuous("年份", breaks = seq(1880, 2020, 20)) + 
  scale_y_continuous("炎热的天数", limits = c(0, 40), 
                     breaks = seq(0, 40, 10),
                     oob = scales::rescale_none) + 
  geom_smooth(size = 1.5, color = "red") + 
  labs(title = "奥地利克雷姆斯明斯特: 1876—2018年炎热的天数", 
       subtitle = "每年最高气温高于30˚C的天数，红线为Loess拟合趋势", 
       caption = "数据来源：https://www.ecad.eu/")

年度超出30度的温度累计值在飙升

我们通过计算一年中超过30˚C的每日最高温度的累计来构建一个有多热的度量标准。这不仅仅可以了解每年有多少炎热天气，还能了解累计的多余热量。

krem_high %>% 
  mutate(
    excess_heat = ifelse(Temp - 30 > 0, Temp-30, 0)
  ) %>% 
  group_by(year) %>% 
  summarise(
    excess_degress_total = sum(excess_heat)
  ) %>% 
  as_tsibble(index = year) %>% 
  fill_gaps() %>% 
  replace_na(list(excess_degress_total = 0)) %>% 
  filter(year <= 2018) %>% 
  ggplot(aes(
    x = year,
    y = excess_degress_total
  )) + 
  geom_line(alpha = 0.7) + 
  geom_smooth(size = 1.5, color = "red") + 
  labs(title = "奥地利克雷姆斯明斯特: 1876—2018年每年的累计超额热量", 
       subtitle = "每年最高气温高于30˚C的超出值的总和，红线为Loess拟合趋势", 
       caption = "数据来源：https://www.ecad.eu/")