时间戳

首页 > 解决方案 > 为什么在 R 环境中操作的 netCDF 文件生成的栅格与在 ArcGIS 中创建的栅格具有不同的分辨率？

r - 为什么在 R 环境中操作的 netCDF 文件生成的栅格与在 ArcGIS 中创建的栅格具有不同的分辨率？

问题描述

我正在处理从哥白尼的海洋环境监测服务下载的 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）。

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

感谢您的帮助！

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

标签： rarraysgisnetcdf4

解决方案

处理带有空间数据的 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进行基于单元格的回归的示例

推荐阅读