将一个netcdf时间变量转换为R日期对象

Question

我有一个带有时间序列的netcdf文件，时间变量具有以下典型的元数据：

    double time(time) ;
            time:standard_name = "time" ;
            time:bounds = "time_bnds" ;
            time:units = "days since 1979-1-1 00:00:00" ;
            time:calendar = "standard" ;
            time:axis = "T" ;

在R中，我想将时间转换为R日期对象。目前我是通过读取单位属性并拆分字符串，然后使用第三个条目作为起始点（假设间隔是“天”，时间是00:00等）来实现这一点。

require("ncdf4")
f1<-nc_open("file.nc")
time<-ncvar_get(f1,"time")
tunits<-ncatt_get(f1,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
dates<-as.Date(time,origin=unlist(tustr)[3])

这个硬编码的解决方案适用于我的特定示例，但我希望在R中有一个包可以很好地处理UNIDATA netcdf日期约定的时间单位，并安全地将它们转换为R日期对象。

- Adrian Tompkins

---

请注意，新提出的、目前正在开发中的神奇 stars 包将自动处理日期，有关示例，请参见第一篇博客文章：http://r-spatial.org/r/2017/11/23/stars1.html。 - AF7

啊，我忘记了补充一点，units 包似乎可以优雅地处理日期。值得一试。 - AF7

请查看我的回答中的修改示例。 - AF7

Answer 1

我刚发现（发布问题两年后！）有一个名为ncdf.tools的软件包，其中包含以下函数：

convertDateNcdf2R

该函数可以将来自netCDF文件或从指定起始日期以来的朱利安日向量（或秒、分钟、小时）转换为POSIXct R向量。

用法：

convertDateNcdf2R(time.source, units = "days", origin = as.POSIXct("1800-01-01", 
    tz = "UTC"), time.format = c("%Y-%m-%d", "%Y-%m-%d %H:%M:%S", 
    "%Y-%m-%d %H:%M", "%Y-%m-%d %Z %H:%M", "%Y-%m-%d %Z %H:%M:%S"))

参数：

time.source 

数字向量或netCDF连接：可以是自起始时间以来的时间单位数，也可以是netCDF文件连接。在后一种情况下，时间向量从netCDF文件中提取。此文件，特别是时间变量，必须遵循CF netCDF约定。

units   

字符串：时间源的单位。如果源是netCDF文件，则此值将被忽略并从该文件中读取。

origin  

POSIXct对象：时间来源的起点或零时刻。 如果源是netCDF文件，则忽略此值并从该文件中读取。

因此，只需将netcdf连接作为第一个参数传递即可，函数会处理其余部分。注意：这仅适用于遵循CF规范的netCDF文件（例如，如果您的单位是“years since”而不是“seconds since”或“days since”，则会失败）。

有关该功能的更多详细信息，请参见此处： https://rdrr.io/cran/ncdf.tools/man/convertDateNcdf2R.html

- Adrian Tompkins

---

ncdf.tools包已被归档。现在有一个名为CFtime的包，它完全支持CF元数据约定中的“时间”维度。 - undefined

Answer 2

据我所知，没有这种功能。我有一个方便的函数，使用 lubridate 实现，基本上与你的一样。

getNcTime <- function(nc) {
    require(lubridate)
    ncdims <- names(nc$dim) #get netcdf dimensions
    timevar <- ncdims[which(ncdims %in% c("time", "Time", "datetime", "Datetime", "date", "Date"))[1]] #find time variable
    times <- ncvar_get(nc, timevar)
    if (length(timevar)==0) stop("ERROR! Could not identify the correct time variable")
    timeatt <- ncatt_get(nc, timevar) #get attributes
    timedef <- strsplit(timeatt$units, " ")[[1]]
    timeunit <- timedef[1]
    tz <- timedef[5]
    timestart <- strsplit(timedef[4], ":")[[1]]
    if (length(timestart) != 3 || timestart[1] > 24 || timestart[2] > 60 || timestart[3] > 60 || any(timestart < 0)) {
        cat("Warning:", timestart, "not a valid start time. Assuming 00:00:00\n")
        warning(paste("Warning:", timestart, "not a valid start time. Assuming 00:00:00\n"))
        timedef[4] <- "00:00:00"
    }
    if (! tz %in% OlsonNames()) {
        cat("Warning:", tz, "not a valid timezone. Assuming UTC\n")
        warning(paste("Warning:", timestart, "not a valid start time. Assuming 00:00:00\n"))
        tz <- "UTC"
    }
    timestart <- ymd_hms(paste(timedef[3], timedef[4]), tz=tz)
    f <- switch(tolower(timeunit), #Find the correct lubridate time function based on the unit
        seconds=seconds, second=seconds, sec=seconds,
        minutes=minutes, minute=minutes, min=minutes,
        hours=hours,     hour=hours,     h=hours,
        days=days,       day=days,       d=days,
        months=months,   month=months,   m=months,
        years=years,     year=years,     yr=years,
        NA
    )
    suppressWarnings(if (is.na(f)) stop("Could not understand the time unit format"))
    timestart + f(times)
}

编辑：大家还可以看看 ncdf4.helpers::nc.get.time.series

编辑2: 注意，新提出并目前正在开发中的神奇 stars 包将自动处理日期，请参见第一篇博客文章以获取示例。

编辑3: 另一种方法是直接使用 units 包，这也是 stars 所使用的。你可以像这样做：（仍无法正确处理日历，我不确定 units 是否可以）

getNcTime <- function(nc) { ##NEW VERSION, with the units package
    require(units)
    require(ncdf4)
    options(warn=1) #show warnings by default
    if (is.character(nc)) nc <- nc_open(nc)
    ncdims <- names(nc$dim) #get netcdf dimensions
    timevar <- ncdims[which(ncdims %in% c("time", "Time", "datetime", "Datetime", "date", "Date"))] #find (first) time variable
    if (length(timevar) > 1) {
        warning(paste("Found more than one time var. Using the first:", timevar[1]))
        timevar <- timevar[1]
    }
    if (length(timevar)!=1) stop("ERROR! Could not identify the correct time variable")
    times <- ncvar_get(nc, timevar) #get time data
    timeatt <- ncatt_get(nc, timevar) #get attributes
    timeunit <- timeatt$units
    units(times) <- make_unit(timeunit)
    as.POSIXct(time)
}

- AF7

---

注意：AF7的函数和SnowFrog的函数都无法正确处理calendar=365_day属性，而ncdf4.helpers::nc.get.time.series可以使用365天日历！ - tbc

units 包是对 UDUNITS 的封装，它不了解日历 - 这些是在 CF 元数据约定中定义的。使用 CFtime 包来实现一站式解决方案。 - undefined

Answer 3

我无法使用@AF7的函数处理我的文件，所以我编写了自己的函数。下面的函数创建了一个POSIXct日期向量，其中开始日期，时间间隔，单位和长度从nc文件中读取。它适用于许多（但可能不是所有...）形状或形式的nc文件。

 ncdate <- function(nc) {
    ncdims <- names(nc$dim) #Extract dimension names
    timevar <- ncdims[which(ncdims %in% c("time", "Time", "datetime", "Datetime",
                                          "date", "Date"))[1]] # Pick the time dimension
    ntstep <-nc$dim[[timevar]]$len
    tm <- ncvar_get(nc, timevar) # Extract the timestep count
    tunits <- ncatt_get(nc, timevar, "units") # Extract the long name of units
    tspace <- tm[2] - tm[1] # Calculate time period between two timesteps, for the "by" argument 
    tstr <- strsplit(tunits$value, " ") # Extract string components of the time unit
    a<-unlist(tstr[1]) # Isolate the unit .i.e. seconds, hours, days etc.
    uname <- a[which(a %in% c("seconds","hours","days"))[1]] # Check unit
    startd <- as.POSIXct(gsub(paste(uname,'since '),'',tunits$value),format="%Y-%m-%d %H:%M:%S") ## Extract the start / origin date
    tmulti <- 3600 # Declare hourly multiplier for date
    if (uname == "days") tmulti =86400 # Declare daily multiplier for date
    ## Rename "seconds" to "secs" for "by" argument and change the multiplier.
    if (uname == "seconds") {
        uname <- "secs"
        tmulti <- 1 }
    byt <- paste(tspace,uname) # Define the "by" argument
    if (byt == "0.0416666679084301 days") { ## If the unit is "days" but the "by" interval is in hours
    byt= "1 hour"                       ## R won't understand "by < 1" so change by and unit to hour.
    uname = "hours"}
    datev <- seq(from=as.POSIXct(startd+tm[1]*tmulti),by= byt, units=uname,length=ntstep)
}

编辑

针对 @AF7 评论中指出的上述代码仅适用于间隔规则文件的缺陷，datev 可以计算如下：

 datev <- as.POSIXct(tm*tmulti,origin=startd)

- SnowFrog

---

非常感谢 - 我借鉴了一些AF7代码的想法，并将它们合并到我的R脚本中。我想知道是否可以将这样的功能贡献给ncdf4包本身？如果能够内置这样的功能，那将是非常棒的。 - ClimateUnboxed

请注意，这仅适用于时间间隔规则的NetCDF文件，而并非所有NetCDF文件都是如此。为什么我的函数对您无效？我会尝试使它更通用。 - AF7

@AdrianTompkins。以前包中有一个计算日期的功能，但是由于有很多种类型的netcdf文件，它无法与所有文件一起使用，因此开发人员将其删除（感谢David Pierce提供信息）。由于我的功能也会遇到相同的问题，目前AF7的情况也是如此，最好是使这些功能非官方化，并且至少可以帮助其他用户定制自己的功能。 - SnowFrog

谢谢，知道这个非常有用。 - ClimateUnboxed

我问了 tidync 的开发者是否对此感兴趣。这是 Github 上的问题，你可能想在那里表达你的意见：https://github.com/hypertidy/tidync/issues/54#issuecomment-331694920 - AF7

Answer 4

你的期望已经得到满足，package CFtime 可以无缝处理 CF 元数据约定中的 "time" 维度，包括所有定义的日历。

f1 <- nc_open("file.nc")
cf <- CFtime(f1$dim$time$units, f1$dim$time$calendar, f1$dim$time$vals)
dates <- CFtimestamp(cf)

# This works reliably only for 3 of the 9 defined calendars
dates <- as.Date(dates)

CFtimestamp()函数可以正确输出所有可能的日期，包括奇怪的日期"2023-02-30"，但在"360_day"日历上不能输出"2023-03-31"。将其转换为POSIXct类型有些棘手，但你真的需要一个Date类型来处理吗？字符表示形式是否足够？

- Patrick

---

不知道为什么这个回答被踩了。在我看来，这是一个很好的回答！（当然，与此同时我已经从R迁移到了Python...但无论如何;-)） - undefined

@ClimateUnboxed 感谢您接受这个并点赞。祝您在Python的世界中好运！ - undefined

@ClimateUnboxed 最后这可能是一个明智的选择 ;) - undefined

我不太确定，我还在犹豫是否应该直接找朱莉娅，我这颗老化的大脑已经无法再承受接下来几年的转换了。自从我博士期间，我已经从IDL->metview->NCL->R->python用于数据处理/绘图，哈哈。 - undefined