当变量 lat/lon 在 R 中存储为矩阵时如何提取一组 NetCDF 值答案

【问题标题】：how to extract a set of NetCDF values when variable lat/lon stored as matrix in R当变量 lat/lon 在 R 中存储为矩阵时如何提取一组 NetCDF 值
【发布时间】：2018-04-16 08:39:37
【问题描述】：

我正在处理 3 维 (x,y,time) NetCDF 文件，其中包含一年中每小时 PM10 浓度估计值。我的目标是提取几个坐标的每小时估计值——这样将是 365days*24hrs=8760 估计值/年/坐标 ——然后平均到每日（365）估计值。

我的脚本（见下文）在 2013 年运行良好，但在 2012 年输出有很多 NA。我注意到的区别是 2012 年文件中的 lon/lat 以矩阵形式存储...

File E:/ENSa.2012.PM10.yearlyrea_.nc (NC_FORMAT_CLASSIC):

     3 variables (excluding dimension variables):
        float lon[x,y]   
            long_name: Longitude
            units: degrees_east
        float lat[x,y]   
            long_name: Latitude
            units: degrees_north
        float PM10[x,y,time]   
            units: ug/m3

     3 dimensions:
        x  Size:701
        y  Size:401
        time  Size:8784   *** is unlimited ***
            units: day as %Y%m%d.%f
            calendar: proleptic_gregorian

head(lon) 
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]
[1,] -25.0 -25.0 -25.0 -25.0 -25.0 -25.0 -25.0 -25.0 -25.0
[2,] -24.9 -24.9 -24.9 -24.9 -24.9 -24.9 -24.9 -24.9 -24.9

对于 2013 文件，lon 像这样“正常”

File E:/ENSa.2013.PM25.yearlyrea.nc (NC_FORMAT_NETCDF4):

     1 variables (excluding dimension variables):
        float PM25[lon,lat,time]   (Chunking: [701,401,1])  
            long_name: PM25
            units: ug
            _FillValue: -999

     3 dimensions:
        lon  Size:701
            standard_name: longitude
            long_name: longitude
            units: degrees_east
            axis: X
        lat  Size:401
            standard_name: latitude
            long_name: latitude
            units: degrees_north
            axis: Y
        time  Size:8760   *** is unlimited ***
            standard_name: time
            long_name: time at end of period
            units: day as %Y%m%d.%f
            calendar: proleptic_gregorian

head(lon) 
[1] -25.0 -24.9 -24.8 -24.7 -24.6 -24.5

我正在使用以下脚本：

# Command brick reads all layers (time slices) in the file
  pm102013 <- brick("ENSa.2013.PM10.yearlyrea.nc", varname = "PM10")

# Get date index from the file
  idx <- getZ(pm102013)

# Put coordinates and extract values for all time steps
  coords <- matrix(c( -2.094278,    -1.830583,  -2.584482,  -0.175269,  -3.17625,   0.54797,    -2.678731,  -1.433611,  -1.456944,  -3.182186,  
 57.15736,  52.511722,  51.462839,  51.54421,   51.48178,   51.374264,  51.638094,  53.230583,  53.231722,  55.945589),
 ncol = 2) # longitude and latitude
 vals <- extract(pm102013, coords, df=T)

# Merge dates and values and fix data frame names
 df.pm102013 <- data.frame(idx, t(vals)[-1,])
 rownames(df.pm102013) <- NULL
 names(df.pm102013) <- c('date','UKA00399', 'UKA00479', 'UKA00494', 'UKA00259', 'UKA00217', 'UKA00553', 'UKA00515', 'UKA00530', 'UKA00529', 'UKA00454')

#output
 options(max.print=100000000)
 sink("PM10_2013.txt")
 print(df.pm102013)
 sink()

有人知道有办法“解决”经度/纬度问题吗？或者还有另一种有效的方法来提取和平均数据？

【问题讨论】：

标签： r raster netcdf

【解决方案1】：

您可以从 bash 的命令行中提取到位置 lon/lat 的最近点并使用 CDO 计算每日平均值：

lon=34.4
lat=22.1
cdo daymean -remapnn,lon=${lon}/lat=${lat} input.nc output_${lon}_${lat}.nc

remapnn 上的减号表示将结果通过管道传送到 daymean 命令。您可以将其放入 bash 中的每个所需点的循环中。

【讨论】：