一、下载物种的经纬度、水深
# 1.从gbif下载
library(rgbif)
gbif_download <- occ_download(pred_default(),pred("taxonKey", 2333562),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)# 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)