在 R 中多年从 netcdf 层获取两个变量的月平均值和总体平均值

Question

提前道歉，我对使用 netcdf 文件很陌生，想象这很简单，但我似乎无法破解它。我对变量、维度和实际数据感到困惑。我想要尽可能简单的代码。

我有一个来自 Copernicus 的 netcdf 文件，当通过以下方式进入 R 时，该文件的结构如下： chl <- nc_open(file) 使用 ncdf4 包。

所以我有 2 个变量 - 叶绿素和净初级生产力。 4个维度：时间、深度、纬度、经度。

 2 variables (excluding dimension variables):
    float chl[longitude,latitude,depth,time]   
        long_name: Total Chlorophyll
        standard_name: mass_concentration_of_chlorophyll_a_in_sea_water
        units: mg m-3
        unit_long: milligram of Chlorophyll per cubic meter
        _FillValue: 9.96920996838687e+36
        _ChunkSizes: 1
         _ChunkSizes: 15
         _ChunkSizes: 137
         _ChunkSizes: 288
    float nppv[longitude,latitude,depth,time]   
        long_name: Total Primary Production of Phyto
        standard_name: net_primary_production_of_biomass_expressed_as_carbon_per_unit_volume_in_sea_water
        units: mg m-3 day-1
        unit_long: milligrams of Carbon per cubic meter per day
        _FillValue: 9.96920996838687e+36
        _ChunkSizes: 1
         _ChunkSizes: 15
         _ChunkSizes: 137
         _ChunkSizes: 288

 4 dimensions:
    time  Size:323
        long_name: Time (hours since 1950-01-01)
        standard_name: time
        calendar: gregorian
        units: hours since 1950-01-01 00:00:00
        axis: T
        _ChunkSizes: 1024
        _CoordinateAxisType: Time
        valid_min: 377316
        valid_max: 612503.9375
    depth  Size:1
        valid_min: 0.505760014057159
        valid_max: 0.505760014057159
        units: m
        positive: down
        unit_long: Meters
        long_name: Depth
        standard_name: depth
        axis: Z
        _ChunkSizes: 75
        _CoordinateAxisType: Height
        _CoordinateZisPositive: down
    latitude  Size:153
        valid_min: 34
        valid_max: 72
        step: 0.25
        units: degrees_north
        unit_long: Degrees North
        long_name: Latitude
        standard_name: latitude
        axis: Y
        _ChunkSizes: 681
        _CoordinateAxisType: Lat
    longitude  Size:185
        valid_min: -27
        valid_max: 19
        step: 0.25
        units: degrees_east
        unit_long: Degrees East
        long_name: Longitude
        standard_name: longitude
        axis: X
        _ChunkSizes: 1440
        _CoordinateAxisType: Lon

17 global attributes:
    product: GLOBAL_REANALYSIS_BIO_001_029
    producer: CMEMS - Global Monitoring and Forecasting Centre
    title: Monthly mean fields for product GLOBAL_REANALYSIS_BIO_001_029
    area: GLOBAL
    quality_information_document: http://marine.copernicus.eu/documents/QUID/CMEMS-GLO-QUID-001-029.pdf
    Conventions: CF-1.6
    credit: E.U. Copernicus Marine Service Information (CMEMS)
    contact: servicedesk.cmems@mercator-ocean.eu
    references: http://marine.copernicus.eu
    source: MERCATOR FREEBIORYS2V4
    licence: http://marine.copernicus.eu/services-portfolio/service-commitments-and-licence/
    dataset: global-reanalysis-bio-001-029-monthly
    institution: Mercator Ocean
    product_user_manual: http://marine.copernicus.eu/documents/PUM/CMEMS-GLO-PUM-001-029.pdf
    _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
    comment: 
    history: Data extracted from dataset http://localhost:8080/thredds/dodsC/global-reanalysis-bio-001-029-monthly

我想获得一个简单的代码：

1. 将此 netcdf 文件放入光栅堆栈，每个月一个用于 CHL，一个用于 PP。
2. 将每个“月”（时间属性）拆分为年平均值。
3. 给出整个数据集所有月份的 Chl 的总体平均值和 PP 的平均值。
4. 包含一行代码，该代码对 CHL 和 PP 采用 深度平均平均值。目前有 1 个深度级别，深度大小：1，但我最多可以有 72 个，并且希望在我的代码中尽可能包含深度平均值。**

到目前为止...我已经在代码中做到了这一点：

    # Opened the netcdf file
    nc <- nc_open(file) # Open file
    # Use attributes to get variable attributes
    attributes(nc)$names ```
    # Got the nc variable names
    attributes(nc$var)$names
    [1] "chl"  "nppv"
    
    #Get attributes of variables
    ncatt_get(nc, attributes(nc$var)$names[1])
    ncatt_get(nc, attributes(nc$var)$names[2])
    
    # Get  latitude and longitude values - dimensions
    nc_lat <- ncvar_get( nc, attributes(nc$dim)$names[3])
    nc_lon <- ncvar_get( nc, attributes(nc$dim)$names[4])
    nc_time <- ncvar_get( nc, attributes(nc$dim)$names[1])

但是从这里开始，我很迷茫如何将它变成可用的格式，例如一个光栅。任何帮助将不胜感激。

Answer 1

应该没有必要像你一样直接使用 ncdf4 包（在这个网站上可以看到几个类似的问题）。正常情况下可以b <- brick(file)；但我在这里使用较新的terra 包

 library(terra)
 chl <- rast("filename", "chl")
 pp <- rast("filename", "nppv")

从那里您可以使用例如tapp。但是，如果这可行，您也许可以编辑您的问题和show(chl)，以便更容易看到您正在处理的内容，并且可能更容易描述您想要完成的内容。分享一个示例文件也会很有帮助。