中国气象局热带气旋资料中心（https://tcdata.typhoon.org.cn/zjljsjj.html）提供了历年中国热带气旋最佳路径数据集。据此我们可以合成每次台风的路径以及判断每次台风经过的省市区县。
原始数据介绍
原始数据是这样的，以 CH2024BST.txt 数据为例：
66666 2401 38 0001 2401 0 6 EWINIAR 20250301
2024052400 1 83 1283 1004 13
2024052406 1 90 1273 1004 13
......
2024060200 9 448 1738 992 20
2024060206 9 462 1792 992 20
66666 2402 23 0002 2402 0 3 MALIKSI 20250301
2024053006 1 181 1120 1002 13
2024053012 1 183 1120 1000 15
......
2024060203 0 246 1168 1005 10
2024060206 0 243 1171 1005 10
......
|
比较长的行是头记录，含义如下：
比较短的是热带气旋的最佳路径数据记录，含义如下：
特别需要注意的是这个时间是世界时，也就是 UTC，如果想转换成北京时间（UTC+8）还需要加上 8 个小时。
想要判断每次台风是否经过某个区县。一种简单的方法就是直接判断这些路径坐标是否在该区县内。但是这样可能会忽略中间路径经过的曲线。因此，我们首先应该根据这些坐标点合成热带气旋的路径数据。
热带气旋路径生成
附件中的 data 文件夹就是从上述网址下载的 txt 文件合集。读取合并：
library(tidyverse)
library(sf)
fs::dir_ls("data") %>%
lapply(read_lines) %>%
c() %>%
unlist() %>%
as_tibble() %>%
mutate(len = str_length(value)) -> filedf
filedf
|
由于比较长的行是头记录，比较短的是路径数据，所以我们可以根据长度来给每个气旋路径生成一个标志符：
filedf %>%
mutate(z = case_when(len >= 40 ~ row_number())) %>%
fill(z) %>%
filter(len < 40) %>%
select(-len) %>%
separate(value, into = paste0("X", 1:6)) -> filedf2
filedf2
|
处理时间和经纬度：
filedf2 %>%
mutate(X1 = lubridate::ymd_h(X1) + hours(8)) %>%
select(id = z, date = X1, lat = X3, lon = X4) %>%
mutate(date = lubridate::date(date)) %>%
type_convert() %>%
mutate(lat = lat / 10, lon = lon / 10) -> df
df
|
然后就是路径合成了，如果气旋有多个路径点，那么就可以合成路线：
df %>%
group_by(id, date) %>%
mutate(n = n()) %>%
ungroup() %>%
filter(n > 1) -> df1
df1
|
如果气旋持续多个日期，那么从第一天到第二天的路径很难判断到底是第一天经过的某个区县，还是第二天经过的。所以这里我选择了分日期生成路径线条：
df1 %>%
nest(-c(id, date)) -> nestdf
nestdf
nestdf$data[[1]] %>%
select(lon, lat) %>%
as.matrix() %>%
st_linestring() %>%
st_sfc(crs = 4326)
|
所有的线条：
nestdf$data[[1]] %>%
select(lon, lat) %>%
as.matrix() %>%
st_linestring() %>%
st_sfc(crs = 4326)
nestdf %>%
mutate(geometry = map(data, function(x){
x %>%
select(lon, lat) %>%
as.matrix() %>%
st_linestring() %>%
st_sfc(crs = 4326)
})) -> nestdf2
nestdf2
|
再转换成 sf 对象：
nestdf2 %>%
select(-data) %>%
unnest(geometry) %>%
st_sf(crs = 4326) -> dfsf
dfsf
|
保存：
dfsf %>%
write_rds("台风分日期路径（linestring）.rds")
|
如果路径坐标只有一个，就不必合成了，直接保存即可：
df %>%
group_by(id, date) %>%
mutate(n = n()) %>%
ungroup() %>%
filter(n == 1) %>%
select(-n) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) -> dfb
dfb %>%
write_rds("台风分日期路径（point）.rds")
|
最终合成的台风路径是这样的：
判断每天每个省市区县是否有台风
这里我们以区县为例。计算上述的线条和区县矢量数据的相交情况即可：
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
read_sf("2021行政区划/县.shp") %>%
select(-contains("类型")) %>%
st_transform(mycrs) -> county
dfsf %>%
st_transform(mycrs) -> dfsf
dfsf %>%
mutate(len = st_length(.),
len = as.numeric(len)) %>%
filter(len > 0) %>%
st_intersection(county) -> resdf1
resdf1
dfsf %>%
mutate(len = st_length(.),
len = as.numeric(len)) %>%
filter(len == 0) %>%
st_cast("POINT") -> dfsf2
dfsf2 %>%
st_intersection(county) -> resdf2
bind_rows(
resdf1 %>%
st_drop_geometry() %>%
select(date, 省, 省代码, 市, 市代码, 县, 县代码),
resdf2 %>%
st_drop_geometry() %>%
select(date, 省, 省代码, 市, 市代码, 县, 县代码)
) -> resdf
resdf %>%
distinct() -> resdf
resdf
|
如果想要生成平衡面板，可以先准备一个平衡面板：
(ymd("1949-01-01") + 0:(ymd("2024-12-31") - ymd("1949-01-01"))) %>%
as_tibble() %>%
crossing(
county %>%
st_drop_geometry()
) %>%
rename(date = value) -> crossdf
crossdf
|
然后再和前面的数据连接即可：
crossdf %>%
left_join(
resdf %>%
mutate(是否有台风 = 1)
) -> resdfn
resdfn
|
最后保存数据：
resdfn %>%
mutate(是否有台风 = if_else(is.na(是否有台风), 0, 是否有台风)) %>%
rename(日期 = date) %>%
haven::write_dta("1949～2024年各区县每天是否有台风.dta")
resdfn %>%
mutate(是否有台风 = if_else(is.na(是否有台风), 0, 是否有台风)) %>%
rename(日期 = date) %>%
filter(是否有台风 == 1) %>%
haven::write_dta("1949～2024年各区县每天是否有台风（有台风的天）.dta")
|
点击这里跳转到 RStata 短书平台获取附件：使用 R 语言处理台风路径数据及判断每天每个区县是否有台风
评论