不同层和几个netcdf文件的平均值与R

Question

从 2000 年到 2014 年，我每年有 15 个 netCDF 文件 (.nc)。在一个 nc 文件中，我有 8760 层中一个变量的每小时数据。 3个维度是： 时间（尺寸 8760）， 纬度（尺寸 90）和 经度（大小 180）（2° 分辨率）。

我想计算 2000-2014 年 4 月至 9 月上午 8 点到晚上 7 点之间变量的平均值。

对于一个 .nc 文件，这对应于之间的平均值

• 层时间从 2169（即 01/04/2000 8am）到 2180（即 01/04/2000 7pm）（到 i=2169 到 i+11），
• 然后从 2193（即 2000 年 2 月 4 日上午 8 点）到 2204（即 2000 年 2 月 4 日晚上 7 点）（i+22，i+33）
• 等等……
• ...从 6537（即 2000 年 9 月 30 日上午 8 点）到 6548（即 2000 年 9 月 30 日晚上 7 点）
• 然后是所有 nc 的平均值。文件。

结果应以 3 个维度的 .nc 文件呈现： - 时间（只有一个值作为平均值）， - 纬度（尺寸 90）和 - 经度（大小 180）（2° 分辨率）

然后我可以绘制 2000 年至 2014 年（4 月至 9 月，上午 8 点至晚上 7 点）平均变量的地图。 我能够读取每个 nc 文件，为每个 nc 文件的每个小时做一张地图，但我知道如何制作所需的平均值。如果有人可以帮助我，那就太好了。

我的变量名：dname

这是我的代码作为第一读：

ncin <- nc_open("sfvmro3_hourly_2000.nc")
print(ncin)

lon <- ncvar_get(ncin, "lon")
lon[lon > 180] <- lon[lon > 180] - 360
nlon <- dim(lon)
head(lon)

lat <- ncvar_get(ncin, "lat", verbose = F)
nlat <- dim(lat)
head(lat)

print(c(nlon, nlat))

t <- ncvar_get(ncin, "time")
tunits <- ncatt_get(ncin, "time", "units")
nt <- dim(t)

dname <- "sfvmro3"
var.array <- ncvar_get(ncin, dname)*10^9  # from mol.mol-1 to ppb
dlname <- ncatt_get(ncin, dname, "long_name")
dunits <- ncatt_get(ncin, dname, "units")
fillvalue <- ncatt_get(ncin, dname, "_FillValue")
var.array[var.array == fillvalue$value] <- NA
dim(var.array)

tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tyear = as.integer(unlist(tdstr)[1])
tmonth = as.integer(unlist(tdstr)[2])
tday = as.integer(unlist(tdstr)[3])
chron = chron(t, origin = c(tmonth, tday, tyear))

以下是年度 file.nc 之一的详细信息：

 4 variables (excluding dimension variables):
    double time_bnds[bnds,time]   
    double lat_bnds[bnds,lat]   
    double lon_bnds[bnds,lon]   
    float sfvmro3[lon,lat,time]   
        standard_name: mole_fraction_of_ozone_in_air
        long_name: Ozone Volume Mixing Ratio in the Lowest Model Layer
        units: mole mole-1
        original_name: O_x
        original_units: 1
        history: 2016-04-22T05:20:31Z altered by CMOR: Converted units from '1' to 'mole mole-1'.
        cell_methods: time: point (interval: 30 minutes)
        cell_measures: area: areacella
        missing_value: 1.00000002004088e+20
        _FillValue: 1.00000002004088e+20
        associated_files: ...

 4 dimensions:
    time  Size:8760   *** is unlimited ***
        bounds: time_bnds
        units: days since 1850-01-01
        calendar: noleap
        axis: T
        long_name: time
        standard_name: time
    lat  Size:90
        bounds: lat_bnds
        units: degrees_north
        axis: Y
        long_name: latitude
        standard_name: latitude
    lon  Size:180
        bounds: lon_bnds
        units: degrees_east
        axis: X
        long_name: longitude
        standard_name: longitude
    bnds  Size:2

26 global attributes:
    institution: aaaa
    institute_id: aaaa
    experiment_id: aaaa
    source: aaaa
    model_id: aaaa
    forcing: HG, SA, S
    parent_experiment_id: N/A
    parent_experiment_rip: N/A
    branch_time: 0
    contact: aaa
    history: aaa
    initialization_method: 1
    physics_version: 1
    tracking_id: aaa
    product: output
    experiment: aaa
    frequency: hr
    creation_date: 2016-04-22T05:20:31Z
    Conventions: aaa
    project_id: aaa
    table_id:aaa
    title: aaaa
    parent_experiment: N/A
    modeling_realm: aaa
    realization: 1
    cmor_version: 2.7.1

Answer 1

我知道您的问题有两种不同的可能解决方案。一种是基于对每个 .nc 文件取平均值，然后对其进行加权平均，另一种是获得一个非常大的数组并使用该数组进行平均。

• 第一个可能的解决方案

您阅读的每个 .nc 都会为您提供数组、数组 1、数组 2 等等。同样对于每个数组，您将有一个与数组的一维相关联的时间序列。这意味着 time_serie1 具有 array1 的 POSIXct 格式的所有不同时间。所以首先你必须在那个向量中构建。你有一个，你可以得到一个你想要用于平均值的时间的矢量索引。为此，我会使用 lubridate 包，但这不是必需的。

index1 <- month(time_serie1) < 10 & month(time_serie1) > 3 # this make an index from april to septembre
index1 <- index1 & hour(time_serie1) <= 19 & hour(time_serie1) >= 8 # then you add the hour restriction
mean1 <- apply(array1[,,index1],1:2,mean)

此代码将为您提供第一年平均值的二维数组，您可以将数组和 time_series 放入列表并循环它。然后，您每年都会有一个该年平均值的二维数组，您可以对这些数组进行平均。我所说的“体重”平均值部分是因为如果你这样做，并且在你的平均值中包括二月份，你的平均值将在不同的天数内完成，对于你的例子来说，这是没有必要的，但是如果你使用二月份，那么你必须对用于每个平均值的数据量进行加权。

• 第二种可能的解决方案

因为这个解决方案与另一个解决方案几乎相同，但我更喜欢它。您可以按顺序将所有数组合并为一个大数组，以便时间索引按递增顺序排列，我将此数组称为 BigArray。然后合并每个数组关联的时间序列，我称之为BigTime。并寻找你想要平均的指数，它就完成了。最大的优势是您不必对列表中的数据进行循环，并且您不必关心 2 月份的大小变化。

Index <- month(BigTime) < 10 & month(BigTime) > 3 # this make an index from april to septembre
Index <- Index & hour(BigTime) <= 19 & hour(BigTime) >= 8 # then you add the hour restriction
Mean <- apply(BigArray[,,Index],1:2,mean) 

然后你的价值观就完成了。

在这两种可能的情况下，都会构建一个 2d 数组，如果您想要一个 3d 数组，其一维（时间）只有一个值，请添加该维。如果您想查找更多信息，则取特定时间值的平均值通常称为气象科学中的复合技术。

我希望这能解决你的问题。