从不同压力水平的 HDF 中提取数据

Question

我正在尝试在不同压力水平（825.40198、681.29102、464.16000、316.22699 hPa）的点位置（lon=40、lat=34）处提取名为 Data Fields/OzoneTropColumn 的变量多个 hdf 文件

library(raster)
library(ncdf4)
library(RNetCDF)
# read file 
nc <- nc_open("E:/Ozone/test1.nc")
list_col1 <- as.list(list.files("E:/Ozone/", pattern = "*.hdf",
                                full.names = TRUE))
> attributes(nc$var) #using a single hdf file to check its variables
$names
 [1] "Data Fields/Latitude"                "Data Fields/Longitude"              
 [3] "Data Fields/O3"                      "Data Fields/O3DataCount"            
 [5] "Data Fields/O3Maximum"               "Data Fields/O3Minimum"              
 [7] "Data Fields/O3StdDeviation"          "Data Fields/OzoneTropColumn"        
 [9] "Data Fields/Pressure"                "Data Fields/TotColDensDataCount"    
[11] "Data Fields/TotColDensMaximum"       "Data Fields/TotColDensMinimum"      
[13] "Data Fields/TotColDensStdDeviation"  "Data Fields/TotalColumnDensity"     
[15] "HDFEOS INFORMATION/StructMetadata.0" "HDFEOS INFORMATION/coremetadata" 

> pres <- ncvar_get(nc, "Data Fields/Pressure") #looking at pressure level from single file of hdf
> pres
 [1] 825.40198 681.29102 464.16000 316.22699 215.44400 146.77901 100.00000  68.12950  46.41580  31.62290
[11]  21.54430  14.67800  10.00000   6.81291   4.64160

ncin <- raster::stack(list_col1,
                       varname = "Data Fields/OzoneTropColumn",
                       ncdf=TRUE)
#cannot extract using the following code
o3 <- ncvar_get(list_col1,attributes(list_col1$var)$names[9]) 
"Error in ncvar_get(list_col1, attributes(list_col1$var)$names[9]) : 
  first argument (nc) is not of class ncdf4!"
#tried to extract pressure levels
> prsr <- raster::stack(list_col1,varname = "Data Fields/Pressure",ncdf=TRUE)
"Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'stack': varname: Data Fields/Pressure does not exist in the file. Select one from:
Data Fields/O3, Data Fields/O3DataCount, Data Fields/O3Maximum, Data Fields/O3Minimum, Data Fields/O3StdDeviation, Data Fields/OzoneTropColumn, Data Fields/TotColDensDataCount, Data Fields/TotColDensMaximum, Data Fields/TotColDensMinimum, Data Fields/TotColDensStdDeviation, Data Fields/TotalColumnDensity"

#tried using index
#Point location can also be written as below 1 deg by 1 deg resolution
lonIdx <- which(lon >32 & lon <36)  
latIdx <- which(lat >38 & lat <42)  
presIdx <- which(pres >= 400 & pres <= 900)


#also tried
# Option 2 -- subset using array indexing
o3 <- ncvar_get(list_col1,'Data Fields/OzoneTropColumn')
"Error in ncvar_get(list_col1, "Data Fields/OzoneTropColumn") : 
  first argument (nc) is not of class ncdf4!"
extract2 <- o3[lonIdx, latIdx, presIdx, ]

如何在每个压力级别垂直提取这些值？ (SM=Some value)

我希望在位置（lon=40，lat=34）以以下方式输出：

Pressure             1           2         3          4          5 ....       10
825.40198          SM1           SM2      SM3        SM4         SM5...       SM10
681.29102          SM11          SM12     
464.16000 
316.22699          SM..          SM..      SM..       SM..        SM..        SM..

感谢任何帮助。 谢谢

Answer 1

这可能是 netcdf4 和栅格如何命名文件中每个图层的问题。尝试一次从多个 ncdf 创建多层对象时可能会有些混淆。

我会执行以下操作，仅使用光栅：加载单个 netCDF，使用 stack() 或 brick()。这会将文件加载为 R 中的多层对象。使用names() 来识别臭氧层的名称是什么根据光栅包。

firstraster <- stack("E:/Ozone/test1.nc")
names(firstraster)

一旦找到名称，您就可以将所有对象读取为stack()，提取兴趣点信息，甚至无需将所有层组装到一个堆栈中。

Ozonelayername <- "put name here"
files <- list.files("E:/Ozone/", pattern = "*.hdf", full.names = TRUE)
stacklist <- lapply(files, stack)
Ozonelayerlist <- lapply(stacklist, "[[", "Ozonelayername")
 

上面的行将输出一个栅格对象列表（不是堆栈或砖块，只是普通栅格），其中只有您想要的图层。

现在我们只需要在每个层上执行提取。 sapply() 将为我们将其整齐地格式化为矩阵。

pointsofinterest <-  expand.grid(32:36,38:42)
values <- sapply(Ozonelayerlist, extract, pointsofinterest)

我可以测试它，因为我没有数据，但我认为这会起作用。