如何在R中更改分辨率(或重新网格化)数据

5

我有一个数据集,包含经度、纬度和月均变量(如温度或降水量),涵盖了1961年至1970年。该数据集的分辨率为0.5度经/纬度,覆盖全球,并以.NC文件的形式下载,我使用以下命令在R中提取了数据:

library(ncdf)
f <- open.ncdf("D:/CRU/cru_ts3.21.1961.1970.tmp.dat.nc")
A <- get.var.ncdf(nc=f,varid="tmp")
B <- get.var.ncdf(nc=f,varid="lon")
C <- get.var.ncdf(nc=f,varid="lat")
D <- cbind(expand.grid(B, C))
E <- expand.grid(A)

扩展网格(E)是一个数据表,由31,104,000行的变量组成,扩展网格(D)是一个由259,200个经度/纬度行组成的数据表。如果你将259,200乘以10年乘以12个月,就会得到31,104,000。因此,可以使用以下方法将表E分解为月份值:

Month <- 1
Start <- (Month-1)*(259200)+1
Finish <- (Month*259200)
G <- E[Start:Finish,]
H <- expand.grid(G)
I <- cbind(D,H) 

因此,现在我有一个数据表,包括经度、纬度和变量,是第一个月(即1961年1月)的数据。以下是数据示例:
        lon    lat tmp
49184 -68.25 -55.75 7.5
49185 -67.75 -55.75 7.6
49186 -67.25 -55.75 7.6
49899 -70.75 -55.25 6.8
49900 -70.25 -55.25 7.0
49901 -69.75 -55.25 6.9
49902 -69.25 -55.25 7.1
49903 -68.75 -55.25 6.8
49904 -68.25 -55.25 7.6
49905 -67.75 -55.25 8.2

现在,我有一个问题。当前格网的分辨率为0.5*0.5度,我希望将数据"重新调整",以便分辨率为0.25*0.25度。我不想对数据进行任何特别复杂的处理,所以我只希望0.25的网格取代其所在的0.5网格的值,也就是说每个0.5*0.5网格包含4个0.25*0.25网格,我只想让这4个0.25*0.25网格具有与0.5*0.5网格相同的值。
我已经研究过栅格,但似乎无法进行任何操作。
3个回答

2

在R包raster中存在一种解决方案。具体步骤如下:

library("ncdf4")
library("raster")
nc <- nc_open("my_file.nc")
lon <- ncvar_get(nc, "lon")
lat <- ncvar_get(nc, "lat")
time <- ncvar_get(nc, "time")
dname <- "pre"        ## pre for the short name of precpitation 
nlon <- dim(lon)
nlat <- dim(lat)
nt <- dim(time)
lonlat <- expand.grid(lon, lat)    # make grid of given longitude and latitude 
pr.array <- ncvar_get(nc, dname)
dlname <- ncatt_get(nc, dname, "long_name")
dunits <- ncatt_get(nc, dname, "units")
fillvalue <- ncatt_get(nc, dname, "_FillValue")

pr.vec.long <- as.vector(pr.array)
pr.mat <- matrix(pr.vec.long, nrow = nlon * nlat, ncol = nt)
pr.df <- data.frame(cbind(lonlat, pr.mat))

pr_c <- pr.df[ ,-c(1:2)]
 ### Specific region have been clipped out from global datafile by 
## selecting lon and lat range and extract regridded data at 1lon 1lat
 ## resolution.  

x0 <- seq(67.5, 98.5, by = 1) ## choose different resolution, eg. by = 0.5 
y0 <- seq(6.5, 37.5, by = 1)


m <- cbind(x0, y0)
m <- as.data.frame(m)
s <- rasterFromXYZ(m)
pts <- expand.grid(x0, y0)
pos <- pr.df[ ,c(1:2)]
l_pr <- apply(pr_c, 2, function(x) cbind(pos, x))
colnm = c("x","y","z")
for (j in seq_along(l_pr)){
  colnames(l_pr[[j]]) <- colnm
}

pr_rstr <- lapply(l_pr, function(x) rasterFromXYZ(x))
## Use resample command to regrid the data, here nearest neighbor method can also be chosen by setting method = "ngb"
pr_bn <- lapply(pr_rstr, function(x) resample(x, s, method = "bilinear"))
pr_extr <- lapply(pr_bn, function(x) extract(x, pts))
df_pr <- do.call("cbind", lapply(pr_extr, data.frame))
## write dataframe in csv format
write.csv(df_pr, "my_data_regridded_1.csv")

我希望这可以达到目的。


你如何在netcdf中保存文件? - user2543
@user2543,您可以访问此链接将.csv文件转换为NetCDF格式。 - Pankaj

1

这里有一种使用plyr::ddply()的方法-取决于您想要重新分组的频率以及您的表格大小,可能会比较慢。我将考虑一种使用data.table的方法,应该会更快:

require(plyr)
# make your data frame
I<-data.frame(lat=seq(0.5,1000,0.5),lon=1,tmp=sample(1:100,2000,replace=T))

# make an adjustment grid
k<-expand.grid(c(0,0.25),c(0,0.25),0)

# use plyr:ddply() to expand out each entry into the correponding 4 entries
new_I<-ddply(I,.(lat,lon),function(x)as.list(x)+k)
colnames(new_I)<-c("lat","lon","newlat","newlon","tmp")

head(new_I)

  lat lon newlat newlon tmp
1 0.5   1   0.50   1.00  64
2 0.5   1   0.75   1.00  64
3 0.5   1   0.50   1.25  64
4 0.5   1   0.75   1.25  64
5 1.0   1   1.00   1.00  31
6 1.0   1   1.25   1.00  31

实际上,从时间角度来看,这里有一种更好的方式(虽然它有点像黑客,而且在未来你可能希望进行其他数据处理方面的控制时会给你带来更少的控制),但它需要6.5秒来处理2m >> 8M行。

# make your data frame
I<-data.frame(lat=seq(0.5,1000000,0.5),lon=1,tmp=sample(1:100,2000000,replace=T))

# make an adjustment vector
v<-rep(0.25,times=2000000)

# make 3 new tables, apply the vector appropriately, and rbind
I_latshift<-I
I_lonshift<-I
I_bothshift<-I

I_latshift$lat<-I_latshift$lat+v
I_lonshift$lon<-I_lonshift$lon+v
I_bothshift$lat<-I_bothshift$lat+v
I_bothshift$lon<-I_bothshift$lon+v

I<-rbind(I,I_bothshift,I_latshift,I_lonshift)

# sort it for neatness
I<-I[with(I, order(lat, lon)), ]


head(I)

         lat  lon tmp
1       0.50 1.00   3
6000001 0.50 1.25   3
4000001 0.75 1.00   3
2000001 0.75 1.25   3
2       1.00 1.00  88
6000002 1.00 1.25  88

1
这不是R语言的解决方案,只是想指出您可以在Linux/MAC OS环境下使用CDO轻松地从命令行重新对NetCDF文件进行网格化。根据您的描述,您似乎想要使用最近邻插值,对于一个0.25度的规则网格来说,最近邻插值就是最好的选择。
cdo remapnn,r1440x720 in.nc out.nc

然而,你也可以使用一阶或二阶保守重映射。例如,对于一阶:

cdo remapcon,r1440x720 in.nc out.nc

您可以按照当前的方式读取重新网格化后的数据,并将其导入R中。

如何获取不同度数(0.25、0.5、1、1.5、2、2.5度)的分辨率(r)等效值。 - UseR10085
很简单,只需要计算相应的点数即可。例如,分辨率为1度的常规经纬网具有360个经度点和180个纬度点,因此您只需要使用cdo remapcon,r360x180 in.nc out.nc。 - ClimateUnboxed

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接