# stConstruct multivariable, time-wide
if (require(maps) && require(plm) && require(sf)) {
library(sp)
states.m <- map('state', plot=FALSE, fill=TRUE)
IDs <- sapply(strsplit(states.m$names, ":"), function(x) x[1])
sf = st_as_sf(states.m, IDs=IDs)
row.names(sf) = sf$ID # not needed if sf >= 1.0-13
states <- as(sf, "Spatial")
states=geometry(states)
yrs = 1970:1986
time = as.POSIXct(paste(yrs, "-01-01", sep=""), tz = "GMT")
data("Produc")
# deselect District of Columbia, polygon 8, which is not present in Produc:
Produc.st <- STFDF(states[-8], time, Produc[order(Produc[,2], Produc[,1]),])
# stplot(Produc.st[,,"unemp"], yrs, col.regions = brewer.pal(9, "YlOrRd"),cuts=9)
# example 1: st from long table, with states as Spatial object:
# use Date format for time:
Produc$time = as.Date(paste(yrs, "01", "01", sep = "-"))
# take centroids of states:
xy = coordinates(states[-8])
Produc$x = xy[,1]
Produc$y = xy[,2]
#using stConstruct, use polygon centroids for location:
x = stConstruct(Produc, c("x", "y"), "time", interval = TRUE)
class(x)
stplot(x[,,"unemp"])
# alternatively, pass states as SpatialObj:
Produc$state = gsub("TENNESSE", "TENNESSEE", Produc$state)
Produc$State = gsub("_", " ", tolower(Produc$state))
x = stConstruct(Produc, "State", "time", states)
class(x)
all.equal(x, Produc.st, check.attributes = FALSE)
}
if (require(sf)) {
fname = system.file("shape/nc.shp", package="sf")[1]
nc = as(st_read(fname), "Spatial")
timesList = list(
BIR=c("BIR74", "BIR79"), # sets of variables that belong together
NWBIR=c("NWBIR74", "NWBIR79"), # only separated by space
SID=c("SID74", "SID79")
)
t = as.Date(c("1974-01-01","1979-01-01"))
nc.st = stConstruct(as(nc, "data.frame"), geometry(nc), timesList,
TimeObj = t, interval = TRUE)
}
# stConstruct multivariable, space-wide
if (require(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"
# match station order to names in wide table:
stations = 4:15
wind.loc = wind.loc[match(names(wind[stations]), wind.loc$Code),]
row.names(wind.loc) = wind.loc$Station
# convert to utm zone 29, to be able to do interpolation in
# proper Euclidian (projected) space:
# create time variable
wind$time = ISOdate(wind$year+1900, wind$month, wind$day, 0)
w = STFDF(wind.loc, wind$time,
data.frame(values = as.vector(t(wind[stations]))))
space = list(values = names(wind)[stations])
wind.st = stConstruct(wind[stations], space, wind$time, SpatialObj = wind.loc, interval = TRUE)
all.equal(w, wind.st)
class(wind.st)
}
Run the code above in your browser using DataLab