r—不同压力水平下HDF数据的提取

Question

我试图从多个hdf文件中提取一个名为Data Fields/OzoneTropColumn的变量，在不同压力水平(825.40198,681.29102,464.16000,316.22699hPa)的点位置(lon=40,lat=34)

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)

我可以测试它，因为我没有数据，但我认为这将工作。