【发布时间】：2021-02-28 12:38:07
【问题描述】：
我正在处理从哥白尼的海洋环境监测服务下载的 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）。
你知道这里有什么问题吗？
感谢您的帮助！
编辑：here你可以下载数据。
【问题讨论】：
-
能否请您重写您的问题以专注于手头的问题。这似乎是正确的分辨率。如果是这样，您展示的大部分内容都是无关紧要的，但您没有展示太多可以帮助回答您的问题的内容，那么我们怎么知道。可以提供文件吗？你能显示
raster::brick("file.nc")吗？ -
嗨罗伯特。我只是想展示我是如何创建我的栅格的。如果过程是正确的，那么我想知道为什么 R 和 ArcGIS 会产生不同分辨率的栅格，就是这样。我添加了问题中数据的链接。
-
您能否展示您使用
brick("file.nc")获得的结果以及您的方法获得的结果----以便我们进行比较（并删除不相关的位）？