使用R语言rnoaa包批量获取佛罗里达多公园二十年气象数据求助

刚好我之前用rnoaa处理过批量站点匹配的需求，给你一套适配大量公园场景的完整操作流程，亲测好用：

首先得有所有目标公园的经纬度信息——如果还没整理好，你可以用tidygeocoder包批量通过公园名称获取坐标。先给你一个示例数据框，你替换成自己的完整公园列表就行：

library(tidyverse)
# 示例：佛罗里达州3个公园的经纬度，替换成你的研究站点列表
fl_parks <- tibble(
  park_name = c("Everglades National Park", "Biscayne National Park", "Dry Tortugas National Park"),
  lat = c(25.3213, 25.6338, 24.6281),
  lon = c(-80.9372, -80.0844, -82.8793)
)

用rnoaa的ghcnd_stations()获取全球历史气候网络（GHCN-D）的站点元数据，我们先筛选佛罗里达州范围内、且有近20年连续数据的站点，缩小范围能大幅提高后续匹配效率：

library(rnoaa)
# 获取佛罗里达州的GHCN-D站点，筛选有2003-2023年数据的站点
stations <- ghcnd_stations() %>%
  filter(state == "FL") %>% # 只保留佛罗里达州站点，也可以扩大范围比如lat在24-32、lon在-86到-79
  select(id, name, lat, lon, elevation, start, end) %>%
  filter(start <= 2003, end >= 2023) # 确保站点覆盖近20年

用geosphere包计算球面距离，给每个公园找出距离最近的有效站点。我写了一个小函数批量处理：

library(geosphere)

# 定义函数：输入单个公园的经纬度，返回最近的站点信息
find_closest_station <- function(park_lat, park_lon, stations_df) {
  stations_df %>%
    # 计算每个站点到当前公园的球面距离（单位：米）
    mutate(distance = distHaversine(cbind(lon, lat), cbind(park_lon, park_lat))) %>%
    arrange(distance) %>%
    slice(1) %>% # 取距离最近的第一个站点
    select(id, station_name = name, distance_m = distance, elevation_m = elevation)
}

# 批量应用到所有公园
fl_parks_with_stations <- fl_parks %>%
  rowwise() %>%
  mutate(closest_station = list(find_closest_station(lat, lon, stations))) %>%
  unnest(closest_station)

# 查看匹配结果
print(fl_parks_with_stations)

用ghcnd_search()批量拉取每个匹配站点的气候数据，这里以气温（最高/最低）和降水为例，你可以按需调整变量：

# 设置时间范围：近20年（2003-2023）
start_date <- as.Date("2003-01-01")
end_date <- as.Date("2023-12-31")

# 批量获取数据，用group_map循环处理每个公园-站点组合
climate_data <- fl_parks_with_stations %>%
  group_by(park_name, id) %>%
  group_map(function(.x, .y) {
    ghcnd_search(
      id = .x$id,
      date_min = start_date,
      date_max = end_date,
      var = c("TMAX", "TMIN", "PRCP") # 按需选择气候变量
    ) %>%
      pluck("data") %>% # 提取数据部分
      mutate(park_name = .y$park_name) # 关联对应的公园名称
  }) %>%
  bind_rows()

# 清理数据：NOAA返回的数值是放大10倍的，转成实际单位
climate_data_clean <- climate_data %>%
  mutate(
    value = case_when(
      grepl("TMAX|TMIN", element) ~ value / 10, # 气温从0.1℃转成℃
      grepl("PRCP", element) ~ value / 10, # 降水从0.1mm转成mm
      TRUE ~ value
    ),
    # 给变量改个易懂的名字
    element = recode(element,
      TMAX = "max_temp_c",
      TMIN = "min_temp_c",
      PRCP = "precipitation_mm"
    )
  ) %>%
  pivot_wider(names_from = element, values_from = value) # 转成宽格式方便分析

# 查看清理后的结果
head(climate_data_clean)

• 如果请求数据时遇到限流，可以分批次处理，或者给ghcnd_search()加limit参数控制单次请求量；
• 部分站点可能有缺失值，你可以用tidyr::fill()做简单填充，或者用imputeTS包做专业插补；
• 如果需要更可靠的官方站点，可以在筛选stations时加上str_detect(id, "^USW|^USC")——USW是自动气象站，USC是合作观测站，数据质量更稳定。

内容的提问来源于stack exchange，提问作者antR