物种的栖息温度信息下载（GBIF—OBIS—WOA2018）

一、下载物种的经纬度、水深

# 查物种的taxonKey 
library(rgbif)
name_backbone(name="Akarotaxis nudiceps")

# 1.从gbif下载
gbif_download <- occ_download(pred_default(),pred("taxonKey", 2393155),format = "SIMPLE_CSV")
occ_download_wait(gbif_download)
d <- occ_download_get(gbif_download) %>%
  occ_download_import()
library(dplyr)
# 筛选列：物种名，十进制经度decimalLongitude，十进制纬度decimalLatitude，水深
d1 <- select(d,species,decimalLongitude,decimalLatitude,depth)

# 2.从obis下载
install.packages("remotes")
remotes::install_github("iobis/robis")
library(robis)
records <- occurrence(scientificname = "Aeoliscus strigatus")
r1 <- select(records,species,decimalLongitude,decimalLatitude,depth)

# 3.合并结果。去空值、去重复
library(tidyr)
allocc <-bind_rows(d1,r1) %>%
    drop_na() %>%
    distinct()

二、下载温度数据

下载WOA2018十年平均温度数据（论文使用版本）
论文说明：WOA2018，0.25×0.25弧分分辨率，1955-2017年平均，年平均温度场

三、栖息地温度提取

install.packages(c("ncdf4", "terra", "sf"))
library(ncdf4)   # 读取WOA的NetCDF格式数据
library(terra)   # 空间栅格处理
library(dplyr)
library(tidyr)
library(sf)

# 4. 读取并预处理WOA温度数据
nc <- nc_open("woa18_decav_t00_04.nc")

# 提取变量：温度、经度、纬度、深度
temp_array <- ncvar_get(nc, "t_an")  # 年平均温度数组 [lon, lat, depth]
lon <- ncvar_get(nc, "lon")          # 经度：-179.875 ~ 179.875
lat <- ncvar_get(nc, "lat")          # 纬度：-89.875 ~ 89.875
depth <- ncvar_get(nc, "depth")      # 标准深度层（共102层，单位：米）

nc_close(nc)  # 关闭文件释放内存

# 5. 为每个出现记录匹配对应温度
# 5.1 预处理出现记录：转换为空间点（经度、纬度列转换为经纬列）
allocc_sf <- st_as_sf(allocc, 
                      coords = c("decimalLongitude", "decimalLatitude"), 
                      crs = 4326)  # WGS84坐标系

# 5.2 定义温度匹配函数（按经纬度+最近深度层）
match_woa_temp <- function(point_sf, depth_m, lon_grid, lat_grid, depth_grid, temp_array) {
  # 提取点的经纬度
  pt_lon <- st_coordinates(point_sf)[1]
  pt_lat <- st_coordinates(point_sf)[2]
  
  # 找到最近的经纬度网格索引
  lon_idx <- which.min(abs(lon_grid - pt_lon))    # 找lon_grid里离pt_lon最近的点
  lat_idx <- which.min(abs(lat_grid - pt_lat))
  
  # 找到最接近记录深度的WOA标准深度层
  depth_idx <- which.min(abs(depth_grid - depth_m))
  
  # 提取对应温度值（处理NA值）
  temp_val <- temp_array[lon_idx, lat_idx, depth_idx]
  return(ifelse(is.na(temp_val), NA, temp_val))
}

# 5.3 批量匹配所有记录（耗时取决于记录数，单物种约几秒）
# mapply：逐行循环 + 并行传参
allocc$temperature <- mapply(match_woa_temp,
                             point_sf = split(allocc_sf, 1:nrow(allocc_sf)), # 按行切分
                             depth_m = allocc$depth,
                             MoreArgs = list(lon_grid = lon,    # MoreArgs：所有行共用的大背景数据
                                             lat_grid = lat,
                                             depth_grid = depth,
                                             temp_array = temp_array))

# 5.4 过滤匹配失败的记录（温度为NA）
allocc_clean <- allocc %>% drop_na(temperature)

# 6. 计算物种栖息地温度统计量（完全匹配论文方法）
# 论文说明：
# - 平均栖息地温度：所有记录温度的平均值
# - 下热限(T₀₁)：栖息地温度的第1百分位
# - 上热限(T₉₉)：栖息地温度的第99百分位
species_temp_stats <- allocc_clean %>%
  group_by(species) %>%
  summarise(
    avg_habitat_temp = mean(temperature, na.rm = TRUE),
    T01_lower_thermal_limit = quantile(temperature, 0.01, na.rm = TRUE),
    T99_upper_thermal_limit = quantile(temperature, 0.99, na.rm = TRUE),
    n_records = n()  # 记录数，用于数据质量评估
  ) %>%
  ungroup()

# 6. 结果输出与保存
print(species_temp_stats)
write.csv(species_temp_stats, "species_thermal_stats.csv", row.names = FALSE)
write.csv(allocc_clean, "all_occurrences_with_temp.csv", row.names = FALSE)

总

library(rgbif)
library(dplyr)
library(robis)
library(tidyr)
library(ncdf4)
library(terra)
library(sf)
# 读取并预处理WOA温度数据
nc <- nc_open("woa18_decav_t00_04.nc")
temp_array <- ncvar_get(nc, "t_an")
lon <- ncvar_get(nc, "lon")
lat <- ncvar_get(nc, "lat")
depth <- ncvar_get(nc, "depth")
nc_close(nc)
rm(nc)
match_woa_temp <- function(point_sf, depth_m, lon_grid, lat_grid, depth_grid, temp_array) {
  pt_lon <- st_coordinates(point_sf)[1]
  pt_lat <- st_coordinates(point_sf)[2]
  lon_idx <- which.min(abs(lon_grid - pt_lon))
  lat_idx <- which.min(abs(lat_grid - pt_lat))
  depth_idx <- which.min(abs(depth_grid - depth_m))
  temp_val <- temp_array[lon_idx, lat_idx, depth_idx]
  return(ifelse(is.na(temp_val), NA, temp_val))
}

# 查taxonKey【修改物种】
name_backbone(name="Anarrhichthys ocellatus")

# 1.从gbif下载【修改taxonKey】
gbif_download <- occ_download(pred_default(),pred("taxonKey", 2392219),format = "SIMPLE_CSV")
occ_download_wait(gbif_download)
d <- occ_download_get(gbif_download) %>%
  occ_download_import()
d1 <- select(d,species,decimalLongitude,decimalLatitude,depth)

# 2.从obis下载【修改物种】
records <- occurrence(scientificname = "Anarrhichthys ocellatus")
r1 <- select(records,species,decimalLongitude,decimalLatitude,depth)

# 3.合并结果
allocc <-bind_rows(d1,r1) %>%
    drop_na() %>%
    distinct()

# 4. 匹配温度
allocc_sf <- st_as_sf(allocc, 
                      coords = c("decimalLongitude", "decimalLatitude"), 
                      crs = 4326)
allocc$temperature <- mapply(match_woa_temp,
                             point_sf = split(allocc_sf, 1:nrow(allocc_sf)),
                             depth_m = allocc$depth,
                             MoreArgs = list(lon_grid = lon,
                                             lat_grid = lat,
                                             depth_grid = depth,
                                             temp_array = temp_array))
allocc_clean <- allocc %>% drop_na(temperature)

# 5. 计算统计量
species_temp_stats <- allocc_clean %>%
  group_by(species) %>%
  summarise(
    avg_habitat_temp = mean(temperature, na.rm = TRUE),
    T01_lower_thermal_limit = quantile(temperature, 0.01, na.rm = TRUE),
    T99_upper_thermal_limit = quantile(temperature, 0.99, na.rm = TRUE),
    n_records = n()
  ) %>%
  ungroup()

# 6. 结果输出与保存
print(species_temp_stats)
write.csv(species_temp_stats, "species_thermal_stats.csv", row.names = FALSE)
write.csv(allocc_clean, "all_occurrences_with_temp.csv", row.names = FALSE)