将多维 NetCDF 读取为 R 中的数据框

Question

我使用存储一个变量并具有以下维度的 netCDF 文件：经度、纬度、时间。 一般来说，我希望将它与我已经在 R 中存储为数据框的不同数据进行比较 - 前两列是 WGS84 中的坐标，接下来是特定时间的值。

所以我写了以下代码。

# since # ncFile$dim$time$units  say:  [1] "days since 1900-1-1"
daysFromDate <- function(data1, data2="1900-01-01")
{
  round(as.numeric(difftime(data1,data2,units = "days")))
}

#study area:
lon <- c(40.25, 48)
lat <- c(16, 24.25)
myTime <- c(daysFromDate("2008-01-16"), daysFromDate("2011-12-31"))
varName <- "spei"



require(ncdf4)
require(RCurl)
x <- getBinaryURL("http://digital.csic.es/bitstream/10261/104742/3/SPEI_01.nc")
ncFile <- nc_open(x)

LonIdx <- which( ncFile$dim$lon$vals >= lon[1] | ncFile$dim$lon$vals <= lon[2])
LatIdx <- which( ncFile$dim$lat$vals >= lat[1] & ncFile$dim$lat$vals <= lat[2])
TimeIdx <- which( ncFile$dim$time$vals >= myTime[1] & ncFile$dim$time$vals <= myTime[2])
MyVariable <- ncvar_get( ncFile, varName)[ LonIdx, LatIdx, TimeIdx]

我认为将返回数据框，以便我能够轻松地操作数据（例如 - 检查相关性或创建绘图）。 不幸的是，已返回 3 维列表。 如何将其重新格式化为具有以下列 X-Y-Time1-Time2-...的数据框...

所以，示例数据如下所示

X Y 2014-01-01 2014-01-02 2014-01-02

50 17 0.5 0.4 0.3

其中 0.5、0.4 和 0.3 是示例变量值

或者可能有不同的解决方案？

Answer 1

好的，请尝试以下代码，但它假定范围是密集填充的。我将lon 测试从or 更改为and

require(ncdf4)

nc <- nc_open("SPEI_01.nc")
print(nc)

lon <- ncvar_get(nc, "lon")
lat <- ncvar_get(nc, "lat")
time <- ncvar_get(nc, "time")


lonIdx <- which( lon >= 40.25 & lon <= 48.00)
latIdx <- which( lat >= 16.00 & lat <= 24.25)

myTime <- c(daysFromDate("2008-01-16"), daysFromDate("2011-12-31"))
timeIdx <- which(time >= myTime[1] & time <= myTime[2])

data <- ncvar_get(nc, "spei")[lonIdx, latIdx, timeIdx]

indices <- expand.grid(lon[lonIdx], lat[latIdx], time[timeIdx])
print(length(indices))
class(indices)
summary(indices)
str(indices)

df <- data.frame(cbind(indices, as.vector(data)))
summary(df)
str(df)

更新

好的，看起来我知道你想要什么，但有直接的解决方案。到目前为止，我得到的是：使用 split() 函数或 data.table 包拆分数据帧。通过 X&Y 拆分后，您将获得小数据帧的列表，其中 X&Y 是给定帧的常数。可能有可能将它们转置并重新组合回来，但我不知道如何。继续将数据作为列处理可能是个好主意，列表是嵌套的，可以展平，这里是 R 中的拆分链接：http://www.uni-kiel.de/psychologie/rexrepos/posts/dfSplitMerge.html

代码，如上例继续

require(data.table)

colnames(df) <- c("X","Y","Time","spei")
df$Time <- as.Date(df$Time, origin="1900-01-01")

dt <- as.data.table(df)
summary(dt)

# Taken from https://github.com/Rdatatable/data.table/issues/1389
# x data.table
# f use `by` argument instead - unlike data.frame
# drop logical default FALSE will include `by` columns in resulting data.tables - unlike data.frame
# by character column names on which split into lists
# flatten logical default FALSE will result in recursive nested list having data.table as leafs
# ... ignored
split.data.table <- function(x, f, drop = FALSE, by, flatten = FALSE, ...){
    if(missing(by) && !missing(f)) by = f
    stopifnot(!missing(by), is.character(by), is.logical(drop), is.logical(flatten), !".ll" %in% names(x), by %in% names(x), !"nm" %in% by)
    if(!flatten){
        .by = by[1L]
        tmp = x[, list(.ll=list(.SD)), by = .by, .SDcols = if(drop) setdiff(names(x), .by) else names(x)]
        setattr(ll <- tmp$.ll, "names", tmp[[.by]])
        if(length(by) > 1L) return(lapply(ll, split.data.table, drop = drop, by = by[-1L])) else return(ll)
    } else {
        tmp = x[, list(.ll=list(.SD)), by=by, .SDcols = if(drop) setdiff(names(x), by) else names(x)]
        setattr(ll <- tmp$.ll, 'names', tmp[, .(nm = paste(.SD, collapse = ".")), by = by, .SDcols = by]$nm)
        return(ll)
    }
}

# here is data.table split
q <- split.data.table(dt, by = c("X","Y"), drop=FALSE)
str(q)

# here is data frame split
qq <- split(df, list(df$X, df$Y))
str(qq)