我希望能够提取每个网格点的时间序列,例如从特定纬度和经度提取10或30年的数据。然而,我发现这样做非常缓慢。例如,虽然我可以在0.002秒内读取全球10000个值的单个时间点的切片(维度顺序为纬度、经度、时间),但我需要0.01秒才能从一个点位置读取10个时间点上的值。
## a time series of 10 points from one location:
library(ncdf4)
met.nc <- nc_open("1901.nc")
system.time(a <- ncvar_get(met.nc, "lwdown", start = c(100,100,1),
count = c(1,1,10)))
user system elapsed
0.001 0.000 0.090
## close down session
## a global slice of 10k points from one time
library(ncdf4)
system.time(met.nc <- nc_open("1901.nc"))
system.time(a <- ncvar_get(met.nc, "lwdown", start = c(100,100,1),
count = c(100,100,1)))
user system elapsed
0.002 0.000 0.002
我猜测这些文件是为了优化空间层的阅读而编写的,因为a) 变量的顺序是纬度、经度、时间,b) 这是生成这些文件的气候模型的逻辑顺序,c) 因为全球范围是最常见的可视化方式。
我尝试重新排列变量,使时间排在第一位:
ncpdq -a time,lon,lat 1901.nc 1901_time.nc
(ncpdq来自NCO (netCDF operators)软件)
> library(ncdf4)
## first with the original data set:
> system.time(met.nc <- nc_open("test/1901.nc"))
user system elapsed
0.024 0.045 22.334
> system.time(a <- ncvar_get(met.nc, "lwdown", start = c(100,100,1), count = c(1, 1, 1000))
+ )
user system elapsed
0.005 0.027 14.958
## now with the rearranged dimensions:
> system.time(met_time.nc <- nc_open("test/1901_time.nc"))
user system elapsed
0.025 0.041 16.704
> system.time(a <- ncvar_get(met_time.nc, "lwdown", start = c(100,100,1), count = c(1, 1, 1000)))
user system elapsed
0.001 0.019 9.660
我该如何优化在某一时间点阅读时间序列,而不是一次性读取大面积的图层?例如,如果文件按时间、纬度、经度不同方式编写,那么读取速度是否更快?有没有“简单”的方法来转换netCDF-4文件中维度的顺序?
更新:
(由@mdsumner请求的基准测试)
library(rbenchmark)
library(ncdf4)
nc <- nc_open("1901.nc")
benchmark(timeseries = ncvar_get(nc, "lwdown",
start = c(1, 1, 50),
count = c(10, 10, 100)),
spacechunk = ncvar_get(nc, "lwdown",
start = c(1, 1, 50),
count = c(100, 100, 1)),
replications = 1000)
test replications elapsed relative user.self sys.self user.child
2 spacechunk 1000 0.909 1.000 0.843 0.066 0
1 timeseries 1000 2.211 2.432 1.103 1.105 0
sys.child
2 0
1 0
更新 2:
我已经开始在这里开发一个解决方案。这些代码碎片都在 github.com/ebimodeling/model-drivers/tree/master/met/cruncep 的一组脚本中。
这些脚本仍需要一些工作和组织 - 并非所有脚本都有用。但读取速度非常快。虽然与上述结果不完全相同,但到了一天结束时,我可以立即读取一个1.3TB文件(分辨率为0.5度,每6小时记录一次,时间跨度达100年):
system.time(ts <- ncvar_get(met.nc, "lwdown", start = c(50, 1, 1), count = c(160000, 1, 1)))
user system elapsed
0.004 0.000 0.004
(注:维度顺序已更改,如下所述:当使用ncdf4 :: ncvar_get时,我该如何指定维度顺序?)