当前位置: 首页 > 知识库问答 >
问题:

当变量lat/lon以矩阵形式存储在R中时,如何提取一组NetCDF值

程枫
2023-03-14

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

我的脚本(见下文)在2013年运行良好,但2012年的输出有很多NAs。我注意到的区别是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()

有人知道有办法“解决”纵向/横向问题吗?或者有另一种有效的方法来提取和平均数据?

共有1个答案

徐鸿文
2023-03-14

您可以从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中的一个循环中,用于每个期望的点。

 类似资料:
  • 我对R相对较新。我正在尝试从 netCDF 文件中获取温度数据的不同点(纬度、纬度)的时间序列。我的示例数据文件在这里,这里是小文件。我已经尝试了netCDF包和到目前为止我使用的代码 有人能帮我得到一个时间序列的数据帧(第一列)和另一列中某个特定点(lat,lon)的数据值吗。在这种情况下,我正在寻找一个特定纬度点(并对许多兴趣点重复)和给定变量(在本例中为tasmin)的时间序列(1950-0

  • 这是我第一次使用堆栈溢出,我的编码技术非常糟糕,我正在使用一个历史tos的NetCDF文件。我想提取特定lat和lon的tos数据。我有一个三维数组中的tos数据,lon和lat分别在一个二维矩阵中。问题是,我选择的lon和lat的行列组合与tos数组的行-列组合不一致。下面是我目前掌握的代码 我被困在这里,因为我的纬度和纬度矩阵的行和列数与 tos day1 数组的行和列号不对应。 如果你不明白

  • 我想知道在R中存储(和处理)多元(特别是矩阵值)时间序列的最佳选择是什么。 我有一个大数据框,它存储了所有数据和时间变量(在本例中,作为一列名为年) 以下是我可以想到的,但两种选择都有各自的缺点: > 数据帧列表,例如通过

  • 问题: 我在R中有一个代码,可以从单个Aqua Modis网络CDF文件中提取每月海面温度(SST)值(见下文)。但是,我在一个文件夹中有一批 59 个 Aqua Modis netCDF 文件。 目的: 我的目标是从所有59个netCDF文件的每个netCDF中提取变量的经度、纬度和SST,使用函数stack::raster()将它们转换为光栅文件,然后处理这些文件。 我的数据框有 650 行,

  • 在具有sigmoid激活函数的神经网络的反向传播中, 权重更新规则由以下公式给出: 其中alpha是学习速率,A是前一层的激活, 其中,Y=给定值,Y'由输出层在神经网络中计算 在我的例子中,Y是4x1=[0.3,0.2,0.4,0.1],Y的一个实例是4x1=[0.2,0.1,0.1,0.2] 如何计算D=(Y-Y’’Y’(1-Y’) (Y-Y’)=4x1和Y’=4x1,和1 还是元素乘法?

  • 我已经搜索了很长时间,但仍然无法弄清楚这一点。似乎光栅包是要提取的,但只能从二维数据中提取。 这个四维数据的例子,一个netCDF文件包含连续三天(72小时)的每小时压力水平(4级)气温。https://drive.google.com/file/d/1UIiX9-xHrtH2FT1torg53iPxyzLxSYQu/view?usp=sharing。 我只想提取一些点位置(xy)的温度,以及相应