数据提取

Question

我试图从netcdf文件中提取2000-2003年一个区域(经度-45W到10E，纬度30-40°)和一个点(lat=37 deg，lon=-40)以及整个地理区域的垂直剖面变质量混合比(MMR)的月平均值。任务：

1. 在2000-2003年期间提取该地区和点的月MMR值分别为400-450,450-500,500-550,550-600 hpa。例如，在400 hPa get (纬度、经度、MMR作为光栅层)
2. 作为光栅格式在两个压力级别之间提取月臭氧MMR。

library(ncdf4)
nc = nc_open("E:/ERA5/ERA_2000_2003.nc")
lat=ncvar_get (nc,"latitude") ; lon=ncvar_get(nc,"longitude")
t = ncvar_get(nc, "time")
pres = ncvar_get(nc, "level")
length(lat); length(lon);length(time);length(pres)
t <- as.POSIXct("1900-01-01 00:00")+as.difftime(t,units="hours")
x <- c(400, 450, 500, 550, 600) 

#I want ozone for longitude -45W to 10E and latitude (30-40E) at pressure level 500 
#for all years (12months*4 yrs= 48 months) 

oz = ncvar_get(nc,'o3',start=c(1,1,1,1),count=c(100,-1,x[3],-1))
lonIdx <- which( lon >= -40 & lon <= 10.00);
latIdx <- which( lat >= 30.00 & lat <= 45.00)

我得到以下错误：

1. Error: cannot allocate vector of size 558.5 Mb

oz = ncvar_get(nc,'o3',start=c(lonIdx,latIdx,1,1),count=c(length(lonIdx),length(latIdx),x[1),-1))

# and sometimes error messages like
2.Error in ncvar_get_inner(ncid2use, varid2use, nc$var[[li]]$missval, addOffset,  : 
 Error: 
variable has 4 dims, but start has 264 entries.  They must match!

3. Error in ncvar_get_inner(ncid2use, varid2use, nc$var[[li]]$missval, addOffset, : Error: variable has 4 dims, but start has 324 entries.  They must match!
Traceback:

1. ncvar_get(nc, "o3", start = c(lonIdx, latIdx, 1, 1), count = c(length(lonIdx), 
 .     length(latIdx), -1, -1))
2. ncvar_get_inner(ncid2use, varid2use, nc$var[[li]]$missval, addOffset, 
 .     scaleFact, start = start, count = count, verbose = verbose, 
 .     signedbyte = signedbyte, collapse_degen = collapse_degen, 
 .     raw_datavals = raw_datavals)
3. stop(paste("Error: variable has", ndims, "dims, but start has", 
 .     length(start), "entries.  They must match!"))

任何帮助都将不胜感激。

谢谢

​

​

Answer 1

从netcdf对象中提取感兴趣子集的两个选项是：( 1)使用具有start和count定义的count函数；或( 2)使用ncvar_get提取整个数组，然后使用标准base R数组对子集进行索引。我不太确定如何计算MMR，所以我将展示如何提取o3值并进行一些基本的总结。如果这没有具体地完成这两个任务，希望它有足够的提示让您开始。

读取数据

library(ncdf4)

# read file 
nc <- nc_open("ERA5_2000_03.nc")

# extract variable levels
lat <- ncvar_get (nc, "latitude")
lon <- ncvar_get(nc, "longitude")
t <- ncvar_get(nc, "time")
t <- as.POSIXct("1900-01-01 00:00") + as.difftime(t, units = "hours")
pres <- ncvar_get(nc, "level")

定义感兴趣区域的指数

lonIdx <- which(lon >= -40 & lon <= 10)
latIdx <- which(lat >= 30 & lat <= 45)
presIdx <- which(pres >= 475 & pres <= 525) # subset to only 500
tIdx <- which(t == t) # i.e., no subset of t

提取所有月份和压力水平的o3

# Option 1 -- use ncvar_get()
extract1 <- ncvar_get(nc,
                      'o3',
                      start = c(lonIdx[1],latIdx[1], presIdx[1], tIdx[1]),
                      count = c(length(lonIdx),length(latIdx),length(presIdx), length(tIdx)))


# Option 2 -- subset using array indexing
o3 <- ncvar_get(nc,'o3')
extract2 <- o3[lonIdx, latIdx, presIdx, ]


# Check the outputs are identical
identical(extract1, extract2)
#>TRUE

还可以选择命名数组维度，这有助于使结果直观。

# Name the array dimensions
dimnames(o3) <- list('lon' = lon, 
                     'lat' = lat,
                     'pres' = pres,
                     'time' = as.character(as.Date(t)))

extract3 <- o3[lonIdx,latIdx, , ]

# Monthly mean for each level of pres across the entire region of interest
apply(extract3, c('pres', 'time'), mean)
#       time
# pres    2000-01-01   2000-02-01   2000-03-01
#   400 8.627874e-08 8.303606e-08 9.857230e-08
#   450 7.911534e-08 7.790138e-08 9.043681e-08
#   500 7.421088e-08 7.398450e-08 8.510393e-08
#   550 7.074213e-08 7.102885e-08 8.128490e-08
#   600 6.817185e-08 6.840323e-08 7.853805e-08
#   ...

# o3 levels at a specific long/lat
o3['-40', '37', ,]
#       time
# pres    2000-01-01   2000-02-01   2000-03-01
#   400 8.160814e-08 8.172732e-08 1.006738e-07
#   450 7.688993e-08 7.681632e-08 9.276743e-08
#   500 7.274836e-08 7.514602e-08 8.791778e-08
#   550 6.989851e-08 7.359140e-08 8.330298e-08
#   600 6.867163e-08 7.110260e-08 8.087377e-08