r Terra创建一年的日长的狂人

Question

可以从library（geosphere）的日长函数找到特定纬度的日光长度。一天的长度在一年中有所不同。这是函数day -lengtth（lat，doy）。我想生成一个1/2度用366层生成spatraster，每个doy一层。下面的代码使我结束了一天。 COORDS DATAFRAME似乎与Rast（）示例集中的示例具有相同的结构，但在下面的错误和警告消息中失败。

library(data.table)
library(geosphere)
library(terra)
latVals <- seq(-60, 90, 0.5)
lonVals <- seq(-180, 180, 0.5)

dl_doy <- data.table(NULL)
temp <- vector()
for (j in 1:366) {
  for (i in 1:length(latVals)) {
    temp[i] <- daylength(latVals[i], j)
  }
  new_name <- paste0("day_", j)
  #  dl_doy[, (new_name) := .(temp)]
  dl_doy <- cbind(dl_doy, temp)
}
setnames(dl_doy,  new = paste0("day_", seq(1, 366, 1)))

coords <- data.frame(x = rep(-100, length(latVals)), y = latVals)
coords <- cbind(coords, dl_doy$day_1)  
names(coords) <- c("x", "y", "lyr.1")
test <- rast(coords, type = "xyz")
#Error: [round] invalid extent
#In addition: Warning message:
#In min(dx) : no non-missing arguments to min; returning Inf

Answer 1

我认为您可以得到这样的事情：

library(geosphere)
library(terra)
r <- rast(res=0.5, ymin=-60, nlyr=365)
lat <- yFromRow(r, 1:nrow(r))
d <- sapply(1:365, function(i) daylength(lat, i))
values(r) <- rep(d, each=ncol(r))

替代方法可以是：

r <- rast(res=0.5, ymin=-60)
y <- init(r, "y")
x <- app(y, \(i) apply(i, 1, \(j) daylength(j, 1:365)))

这比第一种方法效率不高，因为它计算每个经度的日程（而不是计算一个经度和使用对于所有纵向）。需要应用，因为day -lengthth仅部分矢量化。

发生的错误是因为您只有一个唯一的“ x”变量 - 使得无法猜测栅格的分辨率。无论如何，当您从常规的光栅开始时，不必使用rast（coords，type =“ xyz”）。看来您没有使用正确的坐标（单元格中心）。