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

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

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

加载包

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

R
1
2
3
4
5
6
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并添加日平均

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
# 整理最低温数据:
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)

绘图:

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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每日平均气温,从图中可以略微看出,气温似乎有上升趋势。

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

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
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年代以来,温度呈上升趋势。

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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年以来温度异常确实有上升的趋势。

每年低于零度的天数

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
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的天数

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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的每日最高温度的累计来构建一个有多热的度量标准。这不仅仅可以了解每年有多少炎热天气,还能了解累计的多余热量。

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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/")

建模

由于tsibble包的改进,fable包最近出了一些问题,所以只能展示略去建模部分的学习了。

问题可以参考fable包的issues-100:Unable to forecast using fable ets #100

# R

评论

程振兴

程振兴 @czxa.top
截止今天,我已经在本博客上写了607.9k个字了!

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×