1 引言
2 加载包
3 数据
4 生成台风螺旋数据
5 地图数据
6 制作静态图
7 制作视频
8 使用 magick 包合成 .mp4 格式视频
9 结语
10 附录：通过
source调用的代码
1 引言
最近一个月，多个台风陆续登陆……这些台风的路径是什么样的？它们到底有多猛烈？见以下视频（和图片）。
需要注意的是，下文中通过 source 调用的代码见文末附录。
2 加载包
library(readxl)
library(dplyr)
library(lubridate)
library(scales)
library(prismatic)
library(ggplot2)
library(ggview)
library(sf)
library(geocn)
library(parallel)
library(furrr)
library(magick)3 数据
数据来自中国天气台风网。
# 处理经纬度的函数
source("Lon_Lat_Decimal.R")
# 读入多个数据
Typhoon_Yagi <- read_excel("Typhoon_Yagi.xlsx")
Typhoon_Yagi$Lon <- sapply(Typhoon_Yagi$Lon, Lon_Lat_Decimal)
Typhoon_Yagi$Lat <- sapply(Typhoon_Yagi$Lat, Lon_Lat_Decimal)
Typhoon_Yagi$Max_Speed <- as.numeric(gsub("m/s", "", Typhoon_Yagi$Max_Speed))
Typhoon_Yagi$Name <- "摩羯" # "Yagi"
nrow(Typhoon_Yagi)
head(Typhoon_Yagi)
Typhoon_Bebinca <- read_excel("Typhoon_Bebinca.xlsx")
Typhoon_Bebinca$Lon <- sapply(Typhoon_Bebinca$Lon, Lon_Lat_Decimal)
Typhoon_Bebinca$Lat <- sapply(Typhoon_Bebinca$Lat, Lon_Lat_Decimal)
Typhoon_Bebinca$Max_Speed <- as.numeric(gsub("m/s", "", Typhoon_Bebinca$Max_Speed))
Typhoon_Bebinca$Name <- "贝碧嘉" # "Bebinca"
nrow(Typhoon_Bebinca)
head(Typhoon_Bebinca)
Typhoon_Pulasan <- read_excel("Typhoon_Pulasan.xlsx")
Typhoon_Pulasan$Lon <- sapply(Typhoon_Pulasan$Lon, Lon_Lat_Decimal)
Typhoon_Pulasan$Lat <- sapply(Typhoon_Pulasan$Lat, Lon_Lat_Decimal)
Typhoon_Pulasan$Max_Speed <- as.numeric(gsub("m/s", "", Typhoon_Pulasan$Max_Speed))
Typhoon_Pulasan$Name <- "普拉桑" # "Pulasan"
nrow(Typhoon_Pulasan)
head(Typhoon_Pulasan)
# 合并数据 ----
Typhoons <- rbind(Typhoon_Yagi, Typhoon_Bebinca, Typhoon_Pulasan)
nrow(Typhoons)3.1 数据处理
最大风速缩放。
# 最大风速缩放
Typhoons$Max_Speed_Scaled <- Typhoons$Max_Speed |>
rescale(to = c(5, 20)) |>
round()4 生成台风螺旋数据
使用带有随机波动螺旋数据代表台风。
# 生成台风螺旋的函数
source("Make_Typhoon_Spiral.R")
# 生成台风螺旋的数据
Typhoons_Data_List <- lapply(
1:nrow(Typhoons),
function(i) {
temp <- Make_Typhoon_Spiral(
n_spiral = Typhoons$Max_Speed_Scaled[i],
npt_spiral = 25,
max_r = Typhoons$Max_Speed[i] / 50,
center_x = Typhoons$Lon[i],
center_y = Typhoons$Lat[i]
)
cbind(
temp,
date_time = Typhoons$Date_Time[i],
max_wind_speed = Typhoons$Max_Speed[i],
name = Typhoons$Name[i]
)
}
)
typhoon_spirals <- bind_rows(Typhoons_Data_List)
typhoon_spirals$date_time <- ymd_hm(typhoon_spirals$date_time)
nrow(typhoon_spirals)4.1 台风数据处理
转为 sf 对象。
# 转为 sf 对象
typhoon_spirals_sf <- typhoon_spirals |>
st_as_sf(
coords = c("x", "y"),
# 投影方式
crs = st_crs(4326)
)
# 每个时间点，台风的中心点
typhoon_spirals_sf_center <- typhoon_spirals_sf |>
group_by(
# 按 date_time 和 name 分组
date_time, name
) |>
summarise(
# 将几何对象合并为一个整体，然后计算几何中心
geometry = st_centroid(st_union(geometry)),
max_wind_speed = max(max_wind_speed),
# 防止保留分组信息，简化输出
.groups = "drop"
)5 地图数据
使用适合于我国地图的 geocn 包。
# 世界地图
world_map <- load_world_continent()
# 我国地图
china_map <- load_cn_province()
china_line_map <- load_cn_tenline()
# 筛选出上海市、海南省
shanghai_hainan_map <- filter(
china_map,
省 %in% c("上海市", "海南省")
)6 制作静态图
6.1 基础 ggplot2 静态图
# 统一设置主题
source("Preset_Theme.R")
# 轴的范围
explim <- c(1, 1) * 0
(xlims <- c(100, 150) + explim)
(ylims <- c(0, 40) + explim)
# 静态图
p0 <- ggplot() +
geom_sf(
data = world_map,
color = "#606063",
fill = "#34353C",
linewidth = 0.3
) +
geom_sf(
data = typhoon_spirals_sf,
aes(group = name, color = max_wind_speed),
shape = 16,
size = 0.3,
alpha = 0.3
) +
geom_sf_label(
data = typhoon_spirals_sf_center |>
group_by(
name
) |>
summarise(
geometry = st_centroid(st_union(geometry)),
max_wind_speed = max(max_wind_speed),
.groups = "drop"
),
aes(
label = name,
fill = max_wind_speed
),
size = 4,
color = c("grey100", "grey100", "grey20"),
hjust = c(0.5, 1, -1.5),
vjust = c(0.5, -3, 1),
alpha = 0.9,
label.r = unit(0, "pt"),
family = "Hiragino Sans GB",
show.legend = FALSE
) +
geom_sf_text(
data = shanghai_hainan_map,
aes(label = 省),
size = 4,
color = "grey100",
alpha = 0.8,
family = "Hiragino Sans GB"
) +
scale_color_distiller(
name = "最大风速",
palette = "RdYlBu",
direction = -1,
breaks = seq(0, 60, 20),
labels = paste(seq(0, 60, 20), "m/s")
) +
scale_fill_distiller(
name = "最大风速",
palette = "RdYlBu",
direction = -1,
breaks = seq(0, 60, 20),
labels = paste(seq(0, 60, 20), "m/s")
) +
labs(
x = NULL,
y = NULL,
caption = "台风追踪 2024 | 图 @ecologyR | 数据 typhoon.weather.com.cn"
) +
coord_sf(
xlim = xlims,
ylim = ylims
)
# 预览
p0
# 存图
ggsave(
"Typhoons_Static.png",
p0,
width = 6,
height = 5.5,
device = "png",
dpi = 600,
bg = "#34353C"
)7 制作视频
7.1 准备工作
# 获取所有时间点，并升序排序
unique_dates <- sort(unique(typhoon_spirals$date_time))
n_dates <- length(unique_dates)
# 创建保存图片的文件夹
folder <- "Saved_PNGs"
if (!dir.exists(folder)) {
dir.create(folder)
}
# 设定所有要保存的 .png 的名称
png_names <- paste0(folder, "/typhoon-", 1:n_dates, ".png")
# 转换为 sf 对象
typhoon_spirals_sf <- typhoon_spirals |>
st_as_sf(
coords = c("x", "y"),
crs = st_crs(4326)
)
# 准备每个时间点的数据
typhoon_data <- lapply(
unique_dates,
function(current_time) {
# 每个个螺旋
filtered_data <- filter(
typhoon_spirals_sf,
date_time == current_time
)
# 中心点
center_data <- filter(
typhoon_spirals_sf_center,
date_time == current_time
)
# 返回列表
list(filtered = filtered_data, center = center_data)
}
)7.2 并行计算
在每个时间点存一张图。
# 设置并行策略
plan(
multisession,
workers = detectCores() - 1
)
# 开始并行计算（并计时）
system.time({
future_map(1:n_dates, function(i) {
# 各时间点数据
typhoon_spirals_filtered <- typhoon_data[[i]]$filtered
typhoon_spirals_center_filtered <- typhoon_data[[i]]$center
# 统一设置主题
source("Preset_Theme.R")
# 绘制图形
p <- ggplot() +
geom_sf(
data = world_map,
color = "#606063",
fill = "#34353C",
linewidth = 0.3
) +
geom_sf(
data = typhoon_spirals_filtered,
aes(group = name, color = max_wind_speed),
shape = 16,
size = 0.3,
alpha = 1
) +
geom_sf_text(
data = typhoon_spirals_center_filtered,
aes(label = name, color = max_wind_speed),
size = 4,
hjust = 1,
vjust = 1,
alpha = 0.5,
family = "Hiragino Sans GB"
) +
geom_sf_text(
data = shanghai_hainan_map,
aes(label = 省),
size = 4,
color = "grey100",
alpha = 0.6,
family = "Hiragino Sans GB"
) +
scale_color_distiller(
name = "最大风速",
palette = "RdYlBu",
direction = -1,
limits = c(10, 70),
breaks = seq(0, 60, 20),
labels = paste(seq(0, 60, 20), "m/s")
) +
labs(
x = NULL,
y = NULL,
caption = paste(
"台风追踪 2024 | 时间",
format(unique_dates[i], "%Y-%m-%d %H"),
"| 数据 https://typhoon.weather.com.cn/"
)
) +
coord_sf(
xlim = c(100, 150),
ylim = c(0, 40)
)
# 保存当前时间点的图片
ggsave(
filename = png_names[i],
plot = p,
width = 6,
height = 5.5,
device = "png",
bg = "#34353C"
)
})
})
# 完成后，停止多核处理
plan(sequential)8 使用 magick 包合成 .mp4 格式视频
# 读回所有图片
all_frames <- image_read(png_names)
# 制作 .mp4 视频（并计时）
system.time({
image_write_video(
all_frames,
"typhoon_animation.mp4",
framerate = 30
)
})9 结语
微信朋友圈 集赞 10 个 以上，可直接发送完整代码及数据。记录之。
10 附录：source 调用的代码
10.1 Preset_Theme.R
预设主题函数。
# 主题设置
Preset_Theme <- function() {
theme_minimal(
base_family = "Hiragino Sans GB",
base_size = 12
) +
theme(
plot.background = element_rect(
color = "#34353C",
fill = "#34353C"
),
panel.background = element_rect(
color = "#606063",
fill = "#34353C",
linewidth = 1
),
plot.caption = element_text(
size = 10,
color = "grey100"
),
panel.grid = element_blank(),
axis.text = element_blank(),
legend.position = "inside",
legend.position.inside = c(0, 1),
legend.justification.inside = c(0, 1),
legend.title = element_text(
size = 10,
color = "grey100"
),
legend.text = element_text(
size = 10,
color = "grey100"
),
legend.key.width = unit(30, "pt"),
legend.key.height = unit(15, "pt"),
legend.direction = "horizontal",
legend.text.position = "top",
legend.ticks = element_line(color = "grey100"),
legend.ticks.length = unit(15, "pt")
)
}
# 设定预设主题
theme_set(Preset_Theme())10.2 Lon_Lat_Decimal.R
定义一个函数，用于将度分格式，转换为十进制。
Lon_Lat_Decimal <- function(degree_str) {
# 提取方向
direction <- substr(degree_str, 1, 1)
# 提取度数部分，并去除度数符号
degree <- as.numeric(gsub("[^0-9.]", "", sub("^[EWSN]:", "", degree_str)))
# 根据方向调整符号
if (direction == "W" || direction == "S") {
degree <- -degree
}
# 返回
return(degree)
}10.3 Lon_Lat_Decimal.R
生成螺旋状箭头数据，添加顺时针和逆时针旋转控制。
Make_Typhoon_Spiral <- function(n_spiral = 100, npt_spiral = 50, max_r = 1, center_x = 0, center_y = 0, clockwise = FALSE) {
data <- data.frame()
for (i in 1:n_spiral) {
# 根据 clockwise 参数决定螺旋方向
if (clockwise) {
angles <- seq(0, 2 * pi, length.out = npt_spiral) + runif(1, 0, 2 * pi)
# 箭头方向变化
u <- -sin(angles)
v <- cos(angles)
} else {
angles <- -seq(0, 2 * pi, length.out = npt_spiral) + runif(1, 0, 2 * pi)
# 箭头方向变化
u <- sin(angles)
v <- -cos(angles)
}
radii <- seq(0, max_r, length.out = npt_spiral) + runif(1, 0, max_r / 0.5)
x <- radii * cos(angles) + center_x # 将中心点移到 center_x
y <- radii * sin(angles) + center_y # 将中心点移到 center_y
temp_data <- data.frame(
x = x,
y = y,
u = u,
v = v,
group = i
)
data <- rbind(data, temp_data)
}
# 返回
return(data)
}