我正在处理一个从哥白尼海洋环境监测服务下载的netcdf4文件，我想我错过了一些东西。。。
考虑一下我在以下代码中从nc文件中提取值的方法。我的范围是用线性回归重新分析存储在时间序列中的数据，以获得每个像素每天的平均SST变化(不太确定我是否得到了我想要的统计数据！我愿意接受建议！无论如何，这不是我的重点(。
library(ncdf4)
library(raster)
nc_data <- nc_open('L4_analyzed_sst_005res_25081981_31122018.nc')
print(nc_data) #it shows metadata
lon <- ncvar_get(nc_data, "lon")
lat <- ncvar_get(nc_data, "lat", verbose = F)
t <- ncvar_get(nc_data, "time")
t1 <- as.POSIXct(t, origin="1981-01-01", format="%Y-%m-%d", tz = "UTC")
sst.array <- ncvar_get(nc_data, "analysed_sst") # store the data in a 3-dimensional array (long, lat, time)
dim(sst.array)
fillvalue <- ncatt_get(nc_data, "analysed_sst", "_FillValue")
fillvalue
nc_close(nc_data)
sst.array[sst.array == fillvalue$value] <- NA #fill the no value with NA
str(sst.array)
##################################################################################################
#1. sst change for the entire period
#######################################################################################
slope_matrix <- matrix(nrow = length(lon), ncol = length(lat))
for (i in 1:dim(sst.array)[1]){
for (j in 1:dim(sst.array)[2]){
value <- sst.array[i,j,]
if (anyNA(value)==FALSE){ #this "if" speed up the process because in my case the only NA are noValue
extracted_data <- data.frame(value=value, t=t1)
extracted_data$quarter <- "null"
extracted_data$quarter <- ifelse(substr(extracted_data$t, 6,7)=="01"|substr(extracted_data$t, 6,7)=="02"|substr(extracted_data$t, 6,7)=="03", "Q1", extracted_data$quarter)
extracted_data$quarter <- ifelse(substr(extracted_data$t, 6,7)=="04"|substr(extracted_data$t, 6,7)=="05"|substr(extracted_data$t, 6,7)=="06", "Q2", extracted_data$quarter)
extracted_data$quarter <- ifelse(substr(extracted_data$t, 6,7)=="07"|substr(extracted_data$t, 6,7)=="08"|substr(extracted_data$t, 6,7)=="09", "Q3", extracted_data$quarter)
extracted_data$quarter <- ifelse(substr(extracted_data$t, 6,7)=="10"|substr(extracted_data$t, 6,7)=="11"|substr(extracted_data$t, 6,7)=="12", "Q4", extracted_data$quarter)
model <- lm(value ~ time(t) + quarter, data = extracted_data)
slope_matrix[i,j] <- mean(predict(model) - extracted_data$value)
}
}
}
slope_raster <- raster(t(slope_matrix), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
plot(slope_raster)
slope_raster <- flip(slope_raster, direction='y')
writeRaster(slope_raster, "sst_change_1981-2018.tif", "GTiff", overwrite=TRUE)
我使用nc文件作为数组，并构建了一个嵌套的for循环，以便对每个像素中的时间序列进行建模。然后，我将值矩阵转换为光栅，将其翻转以重新定向，并将其写成GeoTIFF以供进一步使用。现在，当我在ArcGIS中打开它时，它与用"；使netCDF光栅层"；工具在R中创建的光栅具有稍小的分辨率(0.049比0.050(
你知道这里出了什么问题吗？
谢谢你的帮助！
编辑：您可以在此处下载数据。
有许多更简单的方法来处理带有空间数据的ncdf文件。下面，我打开了带有raster包(具有基于R的代码(和terra包(使用GDAL库(的文件。他们给出了相同的决议；因此可以肯定地假设这是正确的。
f <- "L4_analyzed_sst_005res_25081981_31122018.nc"
library(raster)
b <- brick(f)
b
#class : RasterBrick
#dimensions : 34, 48, 1632, 13643 (nrow, ncol, ncell, nlayers)
#resolution : 0.05004599, 0.05015772 (x, y)
#extent : 13.22879, 15.631, 39.52957, 41.23494 (xmin, xmax, ymin, ymax)
#crs : +proj=longlat +datum=WGS84 +no_defs
#source : L4_analyzed_sst_005res_25081981_31122018.nc
#names : X1981.08.25, X1981.08.26, X1981.08.27, X1981.08.28, X1981.08.29, X1981.08.30, X1981.08.31, X1981.09.01, X1981.09.02, X1981.09.03, X1981.09.04, X1981.09.05, X1981.09.06, X1981.09.07, X1981.09.08, ...
#Date/time : 1981-08-25, 2018-12-31 (min, max)
#varname : analysed_sst
或使用terra
x <- rast(f)
x
#class : SpatRaster
#dimensions : 34, 48, 13643 (nrow, ncol, nlyr)
#resolution : 0.05004599, 0.05015772 (x, y)
#extent : 13.22879, 15.631, 39.52957, 41.23494 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs
#data source : L4_analyzed_sst_005res_25081981_31122018.nc
#names : L4__1, L4__2, L4__3, L4__4, L4__5, L4__6, ...
你使用的分辨率是错误的，因为你这样做了：
raster(nrow=34, ncol=48, xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat))
#class : RasterLayer
#dimensions : 34, 48, 1632 (nrow, ncol, ncell)
#resolution : 0.04900336, 0.04868249 (x, y)
#extent : 13.25381, 15.60598, 39.55465, 41.20986 (xmin, xmax, ymin, ymax)
#crs : +proj=longlat +datum=WGS84 +no_defs
您使用的是角单元中心的lon/lat值。但您需要提供的是角单元外缘的坐标。
与你的问题无关，但你现在可以继续计算四分之一(在任何循环之外——不需要为每个单元格反复计算？(
z <- getZ(b)
month <- as.numeric(substr(z, 6 ,7 ))
quarter <- ((month-1) %% 4 ) + 1
# or quarter <- trunc((month-1) / 3) + 1
table(quarter)
#quarter
# 1 2 3 4
#3434 3333 3435 3441
关于基于细胞的回归的示例，请参见?calc