STplot package: How to slice y-coordinate when using mode 'xt' for spatio-temporal kriging?

Question

I'm trying to interpolate my dataset using spatio-temporal kriging approach from gstat and spacetime packages. For the plotting of the results, I use STplot function. I want to use the mode 'xt' which produce the space/time plots. But when I read about this from the STplot documentation, it states that:

" Beware: when the x-coordinate is plotted, and for each (x,t) [x-coordinate and time - author] element multiple y-coordinates are sent to the plot, it is not clear which (x,y,t) value becomes the plotted value, so slicing single y values is adviced -- no checking is done. "

So how do I do the slicing? It is very much appreciated if anyone can show some examples. Thanks.

my data:

structure(list(LAT = c(-27.6516, -28.125, -28.0402, -30.7834, 
-30.3186, -26.799, -29.1691, -27.8048, -32.0287, -30.7785, -26.8046, 
-28.4746, -33.7363, -27.8916, -27.1171, -28.5705, -31.5168, -26.2128, 
-29.4817, -30.6332, -27.502, -30.0559, -29.9606, -28.3372, -29.8045, 
-30.5696, -28.0601, -32.124, -27.901, -27.1366), LONG = c(122.3147, 
123.9281, 120.8122, 122.4, 119.0079, 117.9917, 117.6532, 120.7011, 
119.9546, 121.8267, 120.268, 117.8563, 117.2423, 121.2715, 119.2709, 
119.0009, 120.4255, 120.3706, 119.7581, 119.3093, 121.1075, 119.3646, 
119.4459, 121.1211, 119.8226, 121.4668, 119.5071, 119.024, 121.1539, 
118.2815), AGE_MA = c(2719L, 2717L, 2712L, 2711L, 2711L, 2711L, 
2710L, 2709L, 2708L, 2705L, 2704L, 2704L, 2702L, 2702L, 2701L, 
2701L, 2701L, 2700L, 2699L, 2699L, 2699L, 2698L, 2697L, 2695L, 
2693L, 2692L, 2691L, 2691L, 2690L, 2689L), EHFI_new = c(1.015, 
0.945, 1.02, -3.165, -2.565, -0.865, -0.65, 0.065, -4.52, -0.375, 
-2.76, 1.44, -2.03, 2.07, -7.915, -5.215, -2.115, 1.575, -8.785, 
-3.185, 2.015, -5.87, -7.855, 2.575, -8.095, 1.42, -6.665, -3.965, 
2.45, 0.665)), row.names = c(NA, 30L), class = "data.frame")

Preparing the session

library("spacetime")
library("gstat")
library("sp")
library("xts")
library("rgdal")

#setting the working directory and load the data
mydata <- read.csv("2627.csv")
df <- data.frame(mydata)
values <- df[, c('EHFI_new')]
time <- df[, c('AGE_MA')]

#Assigning the coordinates and projection
coordinates(df) <- c("LONG","LAT")
proj4string(df) <- CRS("+init=epsg:4326")

#Transform the projection 
sp.df <- spTransform(df,CRSobj = "+proj=utm +zone=51 +south 
+datum=WGS84 
+units=m +no_defs")

Creating object STIDF format

#STIDF consist of dataframe (spatialpoint, temporal, data) 
#spatialPoint 0bject
sp_sp <- SpatialPoints(sp.df@coords, CRS("+proj=utm +zone=51 +south 
+datum=WGS84 +units=m +no_defs"))
#temporal Object
t <- as.POSIXct(time*60*60*24, origin = "1970-01-01")
#data values
data <- data.frame(values=df$EHFI_new)

#'Merging' the objects into the STIDF format
#Creating the STIDF Object - spatial and temporal object, values
stidf <- STIDF(sp_sp, t, data = data)

Making the prediction grid - space and time grid

#space grid
x.range_1 <- as.numeric(c(-608380.858, 860489.815))
y.range_1 <- as.numeric(c(6080795, 7383164))
grd <- expand.grid(x = seq(from = x.range_1[1], to = x.range_1[2], by 
       = 10000), y = seq(from = y.range_1[1], to = y.range_1[2], by = 
       10000)) 
coordinates(grd) <- ~x + y
grd <- SpatialPixels(grd, proj4string = "+proj=utm +zone=51 +south 
+datum=WGS84 +units=m
+no_defs")
#time grid
n <- 5
tgrd <- xts(1:n, seq(min(t)-10, max(t)+100, length = n))
#space-time grid
pred.grd <- STF(grd, tgrd)

Variogram model

vs1 <- variogramST(values~1, stidf, tunit="days", 
               tlags=seq(0,100, length = 7),
               assumeRegular = F, 
               na.omit = T, 
               cores = 4)

#fitting Variogram modell
ssm <- vgmST("simpleSumMetric",
          space = vgm(9,"Gau", 136417, 1.6),
          time = vgm(3,"Gau", 18, 7.5), 
          joint = vgm(3,"Gau", 500, 2), 
          nugget=1.6, stAni= 18677.26)
ssm_vgm <- fit.StVariogram(vs1, ssm, fit.method = 7,
                       method = "L-BFGS-B")

STKriging interpolation

EHFI_krg <- krigeST(values~1, data=stidf, newdata=pred.grd,
                modelList=ssm_vgm, nmax = 30, stAni = 18677.26,
                bufferNmax=3, progress=TRUE) 

Space/time plot using STplot

STmap <- stplot(EHFI_krg, mode = "xt", as.table = T, scaleX = 1)