time units

[edzer]

While experimenting with sf and stars support, I ran into this: the following script

library(sp)
library(spacetime)
library(gstat)

data(wind)
wind.loc$y = as.numeric(char2dms(as.character(wind.loc[["Latitude"]])))
wind.loc$x = as.numeric(char2dms(as.character(wind.loc[["Longitude"]])))
coordinates(wind.loc) = ~x+y
proj4string(wind.loc) = "+proj=longlat +datum=WGS84 +ellps=WGS84"

wind$time = ISOdate(wind$year+1900, wind$month, wind$day)
wind$jday = as.numeric(format(wind$time, '%j'))
stations = 4:15
windsqrt = sqrt(0.5148 * wind[stations]) # knots -> m/s
Jday = 1:366
daymeans = colMeans(
	sapply(split(windsqrt - colMeans(windsqrt), wind$jday), colMeans))
meanwind = lowess(daymeans ~ Jday, f = 0.1)$y[wind$jday]
velocities = apply(windsqrt, 2, function(x) { x - meanwind })
# match order of columns in wind to Code in wind.loc;
# convert to utm zone 29, to be able to do interpolation in
# proper Euclidian (projected) space:
pts = coordinates(wind.loc[match(names(wind[4:15]), wind.loc$Code),])
pts = SpatialPoints(pts)
#if (require(rgdal) && require(maps)) {
proj4string(pts) = "+proj=longlat +datum=WGS84 +ellps=WGS84"
utm29 = CRS("+proj=utm +zone=29 +datum=WGS84 +ellps=WGS84")
pts = spTransform(pts, utm29)
# note the t() in:
w = STFDF(pts, wind$time, data.frame(values = as.vector(t(velocities))))

library(mapdata)
library(maptools)
m = map2SpatialLines(
    map("worldHires", xlim = c(-11,-5.4), ylim = c(51,55.5), plot=F))
proj4string(m) = "+proj=longlat +datum=WGS84 +ellps=WGS84"
m = spTransform(m, utm29)

# setup grid
grd = SpatialPixels(SpatialPoints(makegrid(m, n = 300)),
	proj4string = proj4string(m))
# grd$t = rep(1, nrow(grd))
#coordinates(grd) = ~x1+x2
#gridded(grd)=TRUE

# select april 1961:
w = w[, "1961-04"]

covfn = function(x, y = x) { 
	du = spDists(coordinates(x), coordinates(y))
	t1 = as.numeric(index(x)) # time in seconds
	t2 = as.numeric(index(y)) # time in seconds
	dt = abs(outer(t1, t2, "-"))
	# separable, product covariance model:
	0.6 * exp(-du/750000) * exp(-dt / (1.5 * 3600 * 24))
}

n = 10
tgrd = seq(min(index(w)), max(index(w)), length=n)
pred = krige0(sqrt(values)~1, w, STF(grd, tgrd), covfn)
layout = list(list("sp.points", pts, first=F, cex=.5),
	list("sp.lines", m, col='grey'))
wind.pr0 = STFDF(grd, tgrd, data.frame(var1.pred = pred))

v = vgmST("separable",
          space = vgm(1, "Exp", 750000), 
          time = vgm(1, "Exp", 1.5 * 3600 * 24),
          sill = 0.6)
wind.ST = krigeST(sqrt(values)~1, w, STF(grd, tgrd), v)

all.equal(wind.pr0, wind.ST)
#}

# stars:
library(stars)
df = data.frame(a = rep(NA, 324*10))
tgrd = tgrd + 0.5 * diff(tgrd)[1]
s = STF(grd, tgrd)
newd = addAttrToGeom(s, df)

# 0
(wind.sta = krigeST(sqrt(values)~1, st_as_stars(w), st_as_stars(newd), v))

# attr(v, "temporal unit") = "secs"
(wind.sta = krigeST(sqrt(values)~1, st_as_stars(w), st_as_stars(newd), v))

# 1
plot(stars::st_as_stars(wind.ST), breaks = "equal", col = sf.colors(15), main = "", nbreaks = 16)
# 2
stplot(wind.ST)
# 3
plot(wind.sta, breaks = "equal", col = sf.colors(15), nbreaks = 16, main = "")
st_as_stars(wind.ST)[[1]][1:3,1:3,1]
(wind.sta)[[1]][1:3,1:3,1]
st_bbox(wind.sta)
bbox(wind.ST)
all.equal(wind.sta, stars::st_as_stars(wind.ST), check.attributes = FALSE)

# 4: roundtrip wind.sta->STFDF->stars
rt = stars::st_as_stars(as(wind.sta, "STFDF"))
plot(rt, breaks = "equal", col = sf.colors(15), nbreaks = 16, main = "")
# 5:
stplot(as(wind.sta, "STFDF"), main = "as(wind.sta, STFDF)")
st_bbox(rt)

# 6:
stplot(as(st_as_stars(wind.ST), "STFDF"), main = "ST->stars->ST", col = sf.colors(15), nbreaks = 16)

gives pretty strongly differing outputs when line 84 is outcommented or not. I suspect that the results are wrong when it is outcommented, and right when it is not, but this is not a requirement of the variogram specification. When outcommented, I see

[Using the following time unit: secs]
[Using the following time unit: hours]
[Using the following time unit: secs]

when not outcommented it is secs three times. Also, when tgrd is not shifted half a time step in line 77, the output is three times secs and output looks good (the six plots generated should all be identical).

Please use a stars version from github to test this; there are still some bugs in the stars <--> STFDF conversions in the now released 0.3-0. Any idea @BenGraeler ?

[BenGraeler]

This originates from stVariogramModels.R where in all models

  dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    message("[Using the following time unit: ", units(dt), 

temporal distances are calculated. R uses different units depending on the temporal differences. Subtracting e.g. hours from hours yields differences in hours and the differences can be re-scaled assigning a different unit. Without the attribute "temporal unit" krigeST has no chance to guess whether the distances in the vgm-model coincide with the distances between locations.

The line appears three times, as we calculate the co-variances for covfn.ST(data, model = modelList, separate=separate), v0 = covfn.ST(data, newdata, modelList) and covfn.ST(d0, d0, modelList, separate = FALSE). In order to be at least consistent (three times secs in the above example), krigeST now (#44) assigns the temporal difference unit of units(abs(outer(index(data@time[1]), index(newdata@time[1]), "-"))) to the vgm model. As an alternative, we could force users to provide an attribute "temporal unit" with the vgm model. A model obtained via fit.StVariogram will use the attribute "unit" of the timelag column of the empirical variogram.