为什么在R环境中操纵netCDF文件会产生与ArcGIS中创建的光栅不同的分辨率的光栅？

Question

我正在处理从哥白尼海洋环境监测处下载的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中打开它时，它与用"make netCDF光栅层“工具创建的光栅并不完全重叠。在R中创建的光栅分辨率稍小(0.049比0.050)。

你知道这里有什么问题吗？

谢谢你的帮助!

编辑：这里您可以下载数据。

Answer 1

有更简单的方法来处理具有空间数据的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。