其他
使用 R 语言爬取 1929～2023 年 GSOD 气象站点数据
该课程准备之后有机会再进行直播讲解！感兴趣的小伙伴也可以先自行学习，有问题可以私信我提问！
在徐老师之前的推文 GSOD 全球逐日气象站点数据介绍与下载 中，徐老师介绍了 GSOD 数据以及如何在 R 语言中下载 GSOD 数据。不过实际上 GSOD 网站的结构非常简单，不如直接使用网络数据爬取的方法进行爬取。
GSOD 数据网站首页的链接是：https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/
选择某个年份点击进去就是各个文件的下载链接了：
因此我们的爬取思路就是先爬取首页所有年份的页面，然后从中提取数据链接再下载，最后合并下载得到的文件即可。
爬取所有年份的页面
为了提升爬取速度，这里我们使用并行方法：
library(tidyverse)
library(parallel)
根据自己电脑的内核数创建并行集群：
# 可以运行 detectCores() 查看自己电脑的 CPU 内核数
makeCluster(36) -> cl
# 分发 %>%
clusterExport(cl, "%>%")
下载 1929～2023 年的所有首页数据：
dir.create("title")
parLapply(cl, 1929:2023, fun = function(x){
if(!file.exists(paste0("title/", x, ".html"))){
download.file(paste0("https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/", x, "/"), paste0("title/", x, ".html"))
}
})
这里使用了这样的结构：
if(!file.exists(paste0("title/", x, ".html"))){
....
}
这样就可以通过多次反复运行这段代码来确保所有的文件都下载成功，并且自动跳过已经下载成功的。
然而这里还存在一种难以筛选的情况，也就是仅仅下载了文件的开头（下载不完全），使用下面的代码筛选没完成成功的再下载就行：
# 重新下载那些没有完全下载的
parLapply(cl, 1929:2023, fun = function(x){
readLines(paste0("title/", x, ".html")) -> text
if(!stringr::str_detect(text[length(text)], "script")){
download.file(paste0("https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/", x, "/"), paste0("title/", x, ".html"))
}
})
这里我是逐次读取了每个文件，判断该文件的最后一行是否有“script”字符（如果有的话，说明该文件下载完全了）。
反复运行上面的代码就可以成功下载所有年份的首页 html 文件了。
从 首页 html 里面提取数据下载链接
我们先处理一个文件：
library(rvest)
library(tidyverse)
filelink <- "title/2020.html"
read_html(filelink) -> html
str_extract(filelink, "\\d{4}") -> year
html %>%
html_elements("a") %>%
html_attr("href") %>%
as_tibble() %>%
dplyr::filter(str_detect(value, "\\.csv")) %>%
mutate(value = paste0("https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/", year, "/", value)) %>%
mutate(年份 = year)
#> # A tibble: 12,299 × 2
#> value 年份
#> <chr> <chr>
#> 1 https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/2020/0… 2020
#> 2 https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/2020/0… 2020
#> 3 https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/2020/0… 2020
#> 4 https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/2020/0… 2020
#> 5 https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/2020/0… 2020
#> 6 https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/2020/0… 2020
#> 7 https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/2020/0… 2020
#> 8 https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/2020/0… 2020
#> 9 https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/2020/0… 2020
#> 10 https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/2020/0… 2020
#> # ℹ 12,289 more rows
这样我们就获取了 2020 年所有数据文件的链接，循环所有的文件即可获取全部年份的：
lapply(fs::dir_ls("title"), function(x){
read_html(x) -> html
str_extract(x, "\\d{4}") -> year
html %>%
html_elements("a") %>%
html_attr("href") %>%
as_tibble() %>%
dplyr::filter(str_detect(value, "\\.csv")) %>%
mutate(value = paste0("https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/", year, "/", value)) %>%
mutate(年份 = year)
}) %>%
bind_rows() -> filelist
下载所有的数据文件
下面就可以下载所有的文件了，首先我们创建一个文件夹保存分年份的数据：
dir.create("data")
for (y in 1929:2023) {
dir.create(paste0("data/", y))
}
同样，我们使用多线程的方式进行下载：
# 分发 filelist 对象
clusterExport(cl, "filelist")
parLapply(cl, 1:nrow(filelist), fun = function(x){
stringr::str_remove_all(filelist$value[x], "https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/") %>%
stringr::str_replace_all("/", "_") -> y
if(!file.exists(paste0("data/", filelist$年份[x], "/", y)) | file.size(paste0("data/", filelist$年份[x], "/", y)) < 10){
download.file(filelist$value[x],
destfile = paste0("data/", filelist$年份[x], "/", y))
}
}) -> res
这里的 if 我用了两个条件进行判断是否下载，一个是文件不存在，另外一个是文件大小小于 10kb（太小的文件很可能是个空文件或者下载不完全的）。
反复运行上面的代码即可下载全部的数据文件。
文件较多，耗时较久。所以反复运行上面的代码可能会很费事，建议每次运行完之后从 filelist 里面先剔除运行成功的：
filelist %>%
mutate(link = value) %>%
mutate(value = str_remove_all(value, "https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/"),
value = str_replace_all(value, "/", "_")) -> filelist2
fs::dir_ls("data", recurse = T, regexp = "[.]csv$", type = "file") -> filefin
filefin %>%
as_tibble() %>%
rename(file = value) %>%
mutate(value = basename(file)) -> filefin2
filelist2 %>%
anti_join(filefin2) -> filelist3
然后再运行没成功的：
clusterExport(cl, "filelist3")
parLapply(cl, 1:nrow(filelist3), fun = function(x){
stringr::str_remove_all(filelist3$link[x], "https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/") %>%
stringr::str_replace_all("/", "_") -> y
if(!file.exists(paste0("data/", filelist3$年份[x], "/", y)) | file.size(paste0("data/", filelist3$年份[x], "/", y)) < 10){
download.file(filelist3$link[x],
destfile = paste0("data/", filelist3$年份[x], "/", y))
}
}) -> res
合并全部的数据
最后我们再合并全部的数据生成一个大文件（总数据会多大几十 GB，如果电脑内存不大，可以每次只合并一个年份的）：
# 合并所有的文件
fs::dir_ls("data", recurse = T,
regexp = "[.]csv$", type = "file") -> filefin
parLapply(cl, filefin, fun = function(x){
readr::read_csv(x) %>%
dplyr::mutate_all(as.character)
}) %>%
bind_rows() -> res
# 临时保存
res %>%
write_rds("res.rds")
然后再把变量重命名成中文变量：
res %>%
rename(气象站代码 = STATION,
日期 = DATE,
纬度 = LATITUDE,
经度 = LONGITUDE,
气象站高程 = ELEVATION,
气象站名称 = NAME,
平均气温 = TEMP,
平均气温属性 = TEMP_ATTRIBUTES,
平均露点 = DEWP,
平均露点属性 = DEWP_ATTRIBUTES,
平均海平面压强 = SLP,
平均海平面压强属性 = SLP_ATTRIBUTES,
平均观测站压强 = STP,
平均观测站压强属性 = STP_ATTRIBUTES,
平均能见度 = VISIB,
平均能见度属性 = VISIB_ATTRIBUTES,
平均风速 = WDSP,
平均风速属性 = WDSP_ATTRIBUTES,
最大持续风速 = MXSPD,
最大持续风速属性 = GUST,
最高气温 = MAX,
最高气温属性 = MAX_ATTRIBUTES,
最低气温 = MIN,
最低气温属性 = MIN_ATTRIBUTES,
降水量 = PRCP,
降水量属性 = PRCP_ATTRIBUTES,
积雪深度 = SNDP,
指示器 = FRSHTT) -> res2
res2 %>%
type_convert() -> res2
res2 %>%
mutate(平均观测站压强 = as.numeric(平均观测站压强)) -> res2
有些变量中使用特殊值作为缺失值，我们把这些值替换成缺失值：
res2 %>%
mutate(平均气温 = (平均气温 - 32) / 1.8,
平均气温 = case_when(平均气温 < 200 ~ 平均气温)) %>%
mutate(平均露点 = (平均露点 - 32) / 1.8,
平均露点 = case_when(平均露点 < 200 ~ 平均露点)) %>%
mutate(平均海平面压强 = case_when(平均海平面压强 < 8000 ~ 平均海平面压强)) %>%
mutate(平均观测站压强 = case_when(平均观测站压强 < 8000 ~ 平均观测站压强)) %>%
mutate(平均能见度 = case_when(平均能见度 < 900 ~ 平均能见度)) %>%
mutate(平均风速 = case_when(平均风速 < 900 ~ 平均风速)) %>%
mutate(最大持续风速 = case_when(最大持续风速 < 900 ~ 最大持续风速)) %>%
mutate(最大持续风速属性 = case_when(最大持续风速属性 < 900 ~ 最大持续风速属性)) %>%
mutate(最高气温 = case_when(最高气温 < 8000 ~ 最高气温)) %>%
mutate(最低气温 = case_when(最低气温 < 8000 ~ 最低气温)) %>%
mutate(降水量 = case_when(降水量 < 90 ~ 降水量)) %>%
mutate(积雪深度 = case_when(积雪深度 < 990 ~ 积雪深度)) %>%
mutate(平均能见度 = 平均能见度 / 0.62137,
平均风速 = 平均风速 * (1852/3600),
最大持续风速 = 最大持续风速 * (1852/3600),
最高气温 = (最高气温 - 32) / 1.8,
最低气温 = (最低气温 - 32) / 1.8,
降水量 = 降水量 * 25.4,
积雪深度 = 积雪深度 * 25.4) %>%
rename(最大阵风 = 最大持续风速属性) %>%
mutate(最大阵风 = 最大阵风 * (1852/3600)) -> res2
res2 %>%
mutate(年份 = lubridate::year(日期)) -> res2
为了方便使用，可以把数据再分年份保存：
res2 %>%
mutate(年份 = lubridate::year(日期)) -> res2
dir.create("GSOD气象数据")
unique(res2$年份) -> yearlist
for (i in yearlist) {
res2 %>%
dplyr::filter(年份 == i) %>%
write_csv(paste0("GSOD气象数据/", i, ".csv"))
}
这样我们就完成了这个数据的爬取。
课程信息
该课程准备之后有机会再进行直播讲解！感兴趣的小伙伴也可以先自行学习，有问题可以私信我提问！
直播地址：腾讯会议(需要报名 RStata 培训班参加) 讲义材料：需要报名 RStata 培训班，详情可阅读：一起来学习 R 语言和 Stata 啦！学习过程中遇到的问题也可以随时提问！
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：
https://rstata.duanshu.com/#/brief/course/d4d95d88409e47cb91fbf6d165a8444c