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

Convert a netcdf time variable to an R date object

提问人:ClimateUnboxed 提问时间:9/1/2017 最后编辑:zx8754ClimateUnboxed 更新时间:10/2/2023 访问量:8933

问:

我有一个带有时间序列的 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 日期对象。目前,我通过读取 units 属性并拆分字符串并使用第三个条目作为我的原点(因此假设间距为“天”,时间为 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 日期对象?

r 日期 netcdf netcdf4 unidata

评论

1赞 AF7 11/24/2017
请注意,新提出且目前正在开发的 awesome 包将自动处理日期,请参阅第一篇博客文章以获取示例:r-spatial.org/r/2017/11/23/stars1.htmlstars
1赞 AF7 12/29/2017
啊,我忘了补充一点,这个包似乎可以优雅地处理日期。值得一试。units
0赞 AF7 12/29/2017
有关示例,请参阅我在答案中的编辑

答:

3赞 AF7 9/2/2017 #1

没有,据我所知。我使用这个方便的函数,它与你的基本相同。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:请注意,新提出且目前正在开发的 awesome 包将自动处理日期,请参阅第一篇博客文章以获取示例。stars

EDIT3:另一种方式是直接使用包,这就是使用。可以做这样的事情:(仍然没有正确处理日历,我不确定是否可以)unitsstarsunits

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)
}

评论

2赞 tbc 11/19/2017
注意:AF7 的函数和 SnowFrog 的函数都无法正确处理该属性,而适用于 365 天日历!calendar=365_dayncdf4.helpers::nc.get.time.series
0赞 Patrick 9/28/2023
该包是 UDUNITS 的包装器,它不知道日历 - 这些日历在 CF 元数据约定中定义。使用套餐提供一站式解决方案。unitsCFtime
3赞 SnowFrog 9/20/2017 #2

我无法让@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)

评论

0赞 ClimateUnboxed 9/21/2017
非常感谢 - 我借用了一些 AF7 代码创意并将它们合并到我的 R 脚本中。我想知道这样的功能是否可以贡献给 ncdf4 包本身?如果将这样的东西内置为标准,那就太好了。
0赞 AF7 9/24/2017
请注意,这仅适用于规则间隔的时间,这不一定适用于所有 NetCDF。为什么我的函数不适合你?我会尽量让它更笼统。
1赞 SnowFrog 9/25/2017
@AdrianTompkins。曾经有一个计算包中日期的函数,但是 netcdfs 的类型太多了,它不适用于所有文件,因此开发人员将其删除(感谢 David Pierce 提供的信息)。由于我的功能与我的功能相同,并且目前与 AF7 相同,因此最好将这些功能设置为非官方的,并且至少可以帮助其他用户自定义自己的功能。
0赞 ClimateUnboxed 9/25/2017
谢谢,知道这一点非常有用
0赞 AF7 9/25/2017
我问开发人员他是否有兴趣。这是 github 问题,您可能想在那里表达您的意见: github.com/hypertidy/tidync/issues/54#issuecomment-331694920tidync
5赞 ClimateUnboxed 12/13/2019 #3

编辑 2023:似乎这个包/答案现在已经过时了,请参阅帕特里克的公认答案以获取执行此操作的新方法。


我刚刚发现(在发布问题两年后!)有一个名为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 连接作为第一个参数传递就足够了,其余参数由函数处理。 注意:仅当 netCDF 文件遵循 CF 约定时,这才有效(例如,如果您的单位是“此后几年”而不是“此后的秒数”或“此后的天数”,它将失败)。

有关该功能的更多详细信息,请访问:https://rdrr.io/cran/ncdf.tools/man/convertDateNcdf2R.html

评论

2赞 Patrick 9/28/2023
软件包已存档。取而代之的是,现在有一个包完全支持 CF 元数据约定的“时间”维度。ncdf.toolsCFtime
3赞 Patrick 9/28/2023 #4

您的希望已经通过CFtime软件包得到了满足。此软件包可以无缝处理 CF 元数据约定的“时间”维度,包括所有定义的日历。

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)

该函数为所有可能的日期提供正确的输出,包括“360_day”日历上的奇怪“2023-02-30”而不是“2023-03-31”。转换为 POSIXct 很棘手,但您真的需要使用还是字符表示会很好?CFtimestamp()Date

评论

0赞 ClimateUnboxed 9/28/2023
不知道为什么这被否决了。对我来说似乎是一个很好的答案!(当然,与此同时,我已经从 R 迁移到了 python......但无论如何;-))
0赞 Patrick 9/28/2023
@ClimateUnboxed 感谢您接受此内容和赞成票。祝你在蟒蛇的世界里好运!
0赞 AF7 10/20/2023
@ClimateUnboxed最终这可能是一个明智的选择;)
0赞 ClimateUnboxed 10/20/2023
我不太确定,我仍然想知道我是否应该直接去找朱莉娅,我老化的大脑在未来几年内无法面临另一次切换,而且自从我获得博士学位以来,我已经从 IDL->metview->NCL->R->python 开始进行处理/绘图:-D