我已从 NASA 海洋颜色项目下载了多个 Lm3 4 km netcdf 文件 (n = 4640)，用于了解海面温度。
class : SpatRaster
dimensions : 766, 709, 1 (nrow, ncol, nlyr)
resolution : 0.04165021, 0.04165796 (x, y)
extent : 24.61, 54.14, 5.8, 37.71 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
source : AQUA_MODIS.20120101.L3m.DAY.SST.x_sst.nc:sst
varname : sst (Sea Surface Temperature)
name : sst
unit : degree_C
我的数据涵盖 2012 年至 2024 年红海 82 公里区域，我想提取特定日期的海面温度点数据。
我已使用下面的代码成功从 netcdf 文件中提取了海面温度值，但我不确定是否只是从第一个 netcdf 层提取值而不是循环遍历所有文件（1:4640），我不确定如何将时间元素（“日期”列）与多年来特定日期的点数据合并。
由于所有权问题，我无法共享数据。
library(terra)
library(lubricate)
#Download Sea Surface Temperature
folder_SST<-"~/Documents/GIS_Data/SST/requested_files"
files_SST <-list.files(folder_SST, pattern='*.nc', full.names ="TRUE")
#Loop through all the file paths and the AQUA MODIS netcdf4 files to gain access to them:
for (file in files_SST) {
nc_SST=open.nc(files_SST)
print.nc(nc_SST)
}
#Select the sea surface temperature values from AQUA MODIS files
SSTs <- terra::rast(files_SST[1], "sst")
#Make the dataframe a spatial object of class = "sf" with a CRS of 4326
Ds_Points <- st_as_sf(x=MyDf,
coords = c("Longitude_E_DD", "Latitude_N_DD"),
crs = 4326)
#Format the dates in the dataframe
#Get dates
MyDf <- MyDf %>% mutate(Dates = as.Date(as.character(MyDf$Date), format = "%d/%m/%Y"))
#Extract the SST values from the AQUA MODIS files
SSTs_Data <- terra::extract(SSTs, Ds_Points)
所需输出：
Survey_Number Dates Longitude_E_DD Latitude_N_DD SST
1 2012-08-01 33.89083 27.26778 23.635
2 2012-06-02 33.86782 27.40854 23.640
3 2012-02-07 33.86230 27.44623 23.690
4 2012-02-12 33.88653 27.26957 23.635
5 2012-02-13 33.88766 27.26848 23.635
6 2012-02-14 33.85000 27.36111 23.780
7 2012-02-15 33.86177 27.41302 23.640
在您的代码中，您确实只从第一个 netcdf 文件中提取值。脚本中的这一行：
SSTs <- terra::rast(files_SST[1], "sst")
意味着您正在加载 files_SST 列表中的第一个文件，并且仅在脚本的其余部分中使用该文件。您将需要使用“for 循环”来迭代每个文件。我在您的 Ds_Points 定义上方没有看到对 MyDf 的引用，因此我使用您所需结果中显示的前 4 个点的纬度/经度值手动创建了 MyDf。我还从给定来源下载了 3 个最新的 SST 文件以用于此目的。
这是我的解决方案：
library(tidyverse)
library(terra)
library(tidyterra)
library(sf)
library(ncdf4)
library(lubridate)
library(stringr)
# create vectors of longs and lats for sample points
# add these to MyDf, then convert MyDf to sf and save sf to Ds_Points
longs <- c(33.89083, 33.86782, 33.86230, 33.88653)
lats <- c(27.26778, 27.40854, 27.44623, 27.26957)
MyDf <- data.frame(longs, lats) %>%
rename(
Longitude_E_DD = longs,
Latitude_N_DD = lats
)
Ds_Points <- st_as_sf(x = MyDf,
coords = c("Longitude_E_DD", "Latitude_N_DD"),
crs = 4326)
######################################################################
# set path to data folder
# get sea surface temperature data files
folder_SST <- paste(here::here(), '/data', sep = '')
files_SST <- list.files(folder_SST, full.names ="TRUE")
# set column names in expected output df,
# then add them to new sst_df
cols <- c('Survey Number', 'Date', 'Longitude_E_DD', 'Latitude_N_DD', 'SST')
sst_df <- data.frame()
for (col in cols) {
sst_df[1, col] <- NA
}
######################################################################
# initialize index and survey_num
index <- 1
survey_num <- 1
# loop through all the AQUA MODIS netcdf4 files
for (file in files_SST) {
# open file as raster and get sst layer
nc_SST = rast(file, "sst")
# get formatted date from file name
date <- as.Date(as.character(str_split(tail(str_split(file, '/')[[1]], 1), '\\.')[[1]][2]), format = "%Y%m%d")
# iterate over each point
for (i in 1:nrow(Ds_Points)) {
# get point
d_point <- Ds_Points[i,]
# get point's lat and long
lon <- st_coordinates(d_point)[,'X']
lat <- st_coordinates(d_point)[,'Y']
# get SST at that point from SST file
# must convert point geometry to spatvector to use in terra::extract function
SST_data <- terra::extract(nc_SST, as_spatvector(d_point, geom = 'geometry', crs = st_crs(Ds_Points)))
# insert data into sst_df at row index
sst_df[index, 'Survey Number'] <- survey_num
sst_df[index, 'Date'] <- as.character(date)
sst_df[index, 'Longitude_E_DD'] <- lon
sst_df[index, 'Latitude_N_DD'] <- lat
sst_df[index, 'SST'] <- SST_data$sst
# next row index
index <- index + 1
}
# next survey number/file
survey_num <- survey_num + 1
}
# show completed sst_df
sst_df
这是输出：
Survey Number Date Longitude_E_DD Latitude_N_DD SST
1 2024-10-28 33.89083 27.26778 26.980
1 2024-10-28 33.86782 27.40854 27.000
1 2024-10-28 33.86230 27.44623 27.075
1 2024-10-28 33.88653 27.26957 26.980
2 2024-10-30 33.89083 27.26778 27.055
2 2024-10-30 33.86782 27.40854 27.005
2 2024-10-30 33.86230 27.44623 27.060
2 2024-10-30 33.88653 27.26957 27.055
3 2024-11-01 33.89083 27.26778 27.265
3 2024-11-01 33.86782 27.40854 27.220
3 2024-11-01 33.86230 27.44623 27.165
3 2024-11-01 33.88653 27.26957 27.265
此脚本查找每个 netcdf 文件的 Ds_Points 中每个点的 SST 值，然后将数据添加到 sst_df。有 4 个点，因此 3 个日期/SST 文件中的每一个都有 4 个 SST 值。该日期是从 netcdf 文件的文件名中提取的。我认为每个日期的调查编号都是唯一的。
请记住，考虑到您有多少个文件，这个特定的解决方案会有点慢，甚至更慢，具体取决于您有多少个点（因为它将迭代每个文件的每个点）。