Search code examples
rterrageosphere

R terra create spatraster of day lengths for a year


The length of daylight at a particular latitude can be found with the daylength function from library(geosphere). The day length varies by the day of the year. Here's the function daylength(lat, doy). I'd like to generate a 1/2 degree spatraster with 366 layers, one for each doy. The code below gets me close for one day. The coords dataframe appears to have the same structure as the example in the rast() set of examples but fails with the error and warning messages just below.

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

Solution

  • I think you can get what you are after like this:

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

    An alternative approach could be:

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

    This is less efficient than the first approach because it computes daylength for each longitude (instead of computing it for one longitude and using that for all longitudes). The apply is needed because daylength is only partially vectorized.

    The error you get occurs because you only have one unique "x" variable -- making it impossible to guess the resolution of the raster. It should not be necessary to use rast(coords, type = "xyz") anyway, as you start out with a regular raster. It seems you are not using the correct coordinates (the center of cells).