Você realmente precisa fazer um pouco mais de pesquisa sobre sua metodologia e ler a documentação para entender a estrutura dos objetos da classe sp S4 e a interação dos objetos sp com as funções relevantes do gstat. Na vinheta sp, há uma explicação detalhada da diferença entre os objetos SpatialPolygons (apenas topologia de polígono) e SpatialPolygonDataFrame (polígonos com atributos).
O que você está explicando não é bloquear o Kriging e usar o tempo como um atributo não resulta em uma estimativa espaço-temporal. A metodologia conceitual que você descreve é bastante inválida. O uso de polígonos ou centróides de polígono viola as suposições de Kriging de um campo aleatório uniforme, anisotropia e não estacionariedade.
Aqui está uma bela vinheta gstat em modelos espaço-temporais usando a interface para o pacote espaço-tempo. Devo também observar que o pacote restritoKriging pode conduzir o Kriging de bloco em blocos de forma arbitrária usando uma função média não estacionária e um variograma isotrópico fracamente estacionário.
Dito isto, para responder sua pergunta, você pode passar um objeto SpatialPointsDataFrame sp diretamente para um modelo de variograma / Kriging no gstat. Nesse tipo de objeto sp, os atributos residem no slot "data" e já estão anexados às coordenadas por meio da estrutura interna da classe S4.
# COERCE meuse DATAFRAME TO sp SpatialPointsDataFrame OBJECT
require(gstat)
data(meuse)
coordinates(meuse) <- ~ x + y
head(meuse@data)
# CREATE SEMIVARIOGRAM USING THE zinc ATTRIBUTE
# NOTE: THERE IS NO ARGUMENT FOR A "4th DIM"
v <- variogram(log(zinc) ~ 1, meuse)
plot(v, type = "l")