在徐老师之前的推文 GSOD 全球逐日气象站点数据介绍与下载 中，徐老师介绍了 GSOD 数据以及如何在 R 语言中下载 GSOD 数据。不过实际上 GSOD 网站的结构非常简单，不如直接使用网络数据爬取的方法进行爬取。
GSOD 数据网站首页的链接是：https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/
选择某个年份点击进去就是各个文件的下载链接了：
因此我们的爬取思路就是先爬取首页所有年份的页面，然后从中提取数据链接再下载，最后合并下载得到的文件即可。
爬取所有年份的页面
为了提升爬取速度，这里我们使用并行方法：
library(tidyverse)
library(parallel)
|
根据自己电脑的内核数创建并行集群：
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)
|
这样我们就获取了 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))
}
|
同样，我们使用多线程的方式进行下载：
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(平均能见度 = 平均能见度 * 1.852,
平均风速 = 平均风速 * (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 短书平台获取附件：使用 R 语言爬取 1929～2023 年 GSOD 气象站点数据
评论