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

Posted 2023-03-11

【中文标题】将 netcdf 时间变量转换为 R 日期对象

【英文标题】：convert a netcdf time variable to an R date object

【发布时间】：2018-02-10 14:20:30

【问题描述】：

我有一个带有时间序列的 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 日期对象？

【问题讨论】：

请注意，新提议且目前正在开发中的超赞stars 包将自动处理日期，请参阅第一篇博文以获取示例：r-spatial.org/r/2017/11/23/stars1.html

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

在我的答案中查看我的编辑以获取示例

【参考方案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 

numeric vector or netCDF connection: 要么是从 origin 开始的时间单位数要么是 netCDF 文件连接，在后一种情况下，时间矢量是从 netCDF 文件中提取的，这个文件，尤其是 time 变量，必须遵循CF netCDF 约定。

units   

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

origin  

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

因此，只需将 netcdf 连接作为第一个参数传递，然后函数处理其余部分就足够了。警告：这仅在 netCDF 文件遵循 CF 约定时才有效（例如，如果您的单位是“之后的年数”而不是“之后的秒数”或“之后的天数”，例如，它将失败）。

更多详情请点击此处： https://rdrr.io/cran/ncdf.tools/man/convertDateNcdf2R.html

【参考方案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

EDIT2：请注意，新提议且目前正在开发中的超赞stars 包将自动处理日期，请参阅the first blog post 示例。

EDIT3：另一种方法是直接使用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 的函数和 SnowFrog 的函数都不能正确处理 calendar=365_day 属性，而 ncdf4.helpers::nc.get.time.series 使用 365 天日历！

【参考方案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)

【讨论】：

非常感谢 - 我借用了一些 AF7 代码创意并将它们合并到我的 R 脚本中。我想知道这样的功能是否可以贡献给 ncdf4 包本身？将这样的东西作为标准内置会很棒。

请注意，这仅适用于规则间隔时间，并非所有 NetCDF 都如此。为什么我的功能对你不起作用？我会尽量让它更通用。

@AdrianTompkins。曾经有一个函数可以计算包中的日期，但是 netcdfs 的类型太多以至于它不适用于所有文件，因此开发人员将其删除（感谢 David Pierce 提供的信息）。因为它与我的功能相同，目前与 AF7 的功能相同，最好让这些功能非官方，至少可以帮助其他用户自定义他们自己的功能。

谢谢，知道这个很有用

我询问了tidync 开发者是否有兴趣。这是github问题，你可能想在那里表达你的意见：github.com/hypertidy/tidync/issues/54#issuecomment-331694920