地温热力图

地温热力图

本文是学习zonination/ground_temp的笔记,讲述了如何绘制一幅日期-地深-地温热力图。

首先准备绘图需要的数据集,数据量很大。
groundtemp.csv
这个数据集是这样的,每一行是一次观测值,每隔五个小时观测一次(所以绘图代码里面 有一个5*60),J.005 J.010 J.016 J.023 J.033 J.046 J.060 J.079 J.103 J.134变量分别代表不同深度的地温。这幅热力图是这样构造的:热力图由许许多多的小矩形构成,横轴是日期,所以为了让矩形块连续,矩形块的宽度应为5个小时,而高度为地深差值。理清了这个思路就可以开始处理数据了:

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
46
47
48
49
50
51
52
53
54
55
56
57
58
library(ggplot2)
library(reshape2)
library(viridis)
temp <- read.csv("groundtemp.csv")
temp <- subset(temp, !is.na(Kolumn5))
temp <- temp[, 1:19]
temp <- melt(temp, id = 1:9)

# 日期中文化
temp$Tid <- gsub(pattern = "Tue", replacement = "四", x = temp$Tid)
temp$Tid <- gsub(pattern = "Wed", replacement = "三", x = temp$Tid)
temp$Tid <- gsub(pattern = "Thu", replacement = "二", x = temp$Tid)
temp$Tid <- gsub(pattern = "Fri", replacement = "五", x = temp$Tid)
temp$Tid <- gsub(pattern = "Sat", replacement = "六", x = temp$Tid)
temp$Tid <- gsub(pattern = "Sun", replacement = "日", x = temp$Tid)
temp$Tid <- gsub(pattern = "Mon", replacement = "一", x = temp$Tid)

temp$Tid <- gsub(pattern = "Jan", replacement = "一月", x = temp$Tid)
temp$Tid <- gsub(pattern = "Feb", replacement = "二月", x = temp$Tid)
temp$Tid <- gsub(pattern = "Mar", replacement = "三月", x = temp$Tid)
temp$Tid <- gsub(pattern = "Apr", replacement = "四月", x = temp$Tid)
temp$Tid <- gsub(pattern = "May", replacement = "五月", x = temp$Tid)
temp$Tid <- gsub(pattern = "Jun", replacement = "六月", x = temp$Tid)
temp$Tid <- gsub(pattern = "Jul", replacement = "七月", x = temp$Tid)
temp$Tid <- gsub(pattern = "Aug", replacement = "八月", x = temp$Tid)
temp$Tid <- gsub(pattern = "Sep", replacement = "九月", x = temp$Tid)
temp$Tid <- gsub(pattern = "Oct", replacement = "十月", x = temp$Tid)
temp$Tid <- gsub(pattern = "Nov", replacement = "十一月", x = temp$Tid)
temp$Tid <- gsub(pattern = "Dec", replacement = "十二月", x = temp$Tid)


temp$Tid <- strptime(temp$Tid, "%a %b %d %H:%M:%S %Y")
# %a: 系统(我的是中文系统)的星期缩写,例如一~六和日
# %b: 系统的月份缩写,例如一月~十二月
# %d: 日,例如01-31
# %H: 小时: 00
# %M: 分钟: 00
# %S: 秒: 00
# %Y: 四位数的年份 0000
temp$variable <- as.numeric(substr(temp$variable, 3, 5))
temp <- subset(temp, !is.na(value))

temp$height <- NA
a <- c(0, unique(temp$variable))
b <- 0
for(n in 1:length(a)){
b[n-1] <- a[n] - a[n-1]
}
b <- c(0, b)
for(n in 1:nrow(temp)){
temp$height[n] <- b[a == temp$variable[n]]
print(paste(signif(n*100/nrow(temp), 4), "% 完成!", sep = ""))
}
# 这个生成的高度变量是个相对高度(地深差),是等下绘制的矩形块的高

# 因为上面的循环需要运行很长时间,所以把 运行结果赶紧保存起来
library(io)
qwrite(temp, "temp.rds")

如果你不想浪费时间运行上面的循环,可以直接用我运行好的:
temp.rds

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
temp <- qread("temp.rds")

theme_set(theme_bw(base_size = 20, base_family = "STSong") + theme(plot.title = element_text(hjust = 0.5)) + theme(panel.grid = element_blank()))
rm(n)
temp$Tid <- as.Date(temp$Tid)
ggplot(temp) +
geom_rect(aes(xmin = Tid,
ymin = variable - height,
ymax = variable,
xmax = Tid + 5*60,
fill = value)) +
scale_fill_viridis("温度") +
scale_y_discrete(breaks = seq(0, 140, 10), limits = c(135, 0)) +
labs(title = "图:64.7N, 20.9E处的地温热力图",
x = "日期",
y = "深度(cm)")

该图片的绘制比较耗时,耐心等待。

# R

评论

Your browser is out-of-date!

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

×