Spatio-temporal Kriging and how extract specific values

Question

I have several years of data that needs to be converted to STDFDF. Shapiro.test and Adf.test talk about non-stationarity and non-normal distribution.

particolato_May2021 <- read.delim(file = "Data.txt")
dput(particolato_May2021) ##sample data

##set the stations of the region
stations <- read.csv("stations_CS.csv", header=TRUE)
names(stations)
coord <-cbind(stations$Longitude,stations$Latitude)
plot(coord)
save(stations,file='data/stations.Rdata')

coords <- SpatialPoints(particolato_May2021[, c("Longitude", "Latitude")]) 
summary(coords)
particolato_sp <- SpatialPointsDataFrame(coords, particolato_May2021) 
names(particolato_sp)
is(particolato_sp)

particolato_May2021$Data <- as.POSIXct(particolato_sp$Data, format = "%Y-%m-%d")
save(particolato_May2021,file='data/particolato.Rdata')

##The problem starts here
STFDF_day <- stConstruct(particolato,space=c('Longitude','Latitude'),time='Data',
                         SpatialObj=SpatialPoints(particolato[,c('Longitude','Latitude')]))
is(STFDF_day)
proj4string(STFDF_day) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
save(STFDF_day,file='data/STFDF_day.Rdata')
load('data/STFDF_day.Rdata') 

#Variogram
var <- variogramST(PM10~1,data=STFDF_day,tunit="days",
                   assumeRegular=F,na.omit=T)

This is an incomplete data set for seven years, and I replaced the missing data with averages for the same days in previous years. I want to conduct space-time kriging and extract data for other unknown points (long, lat) before kriging for the same period for each day. I have following problem:

Found more than one class "xts" in cache; using the first, from namespace 'spacetime'
Also defined by ‘quantmod’
Error in `vectbl_as_row_location()`:
! Must subset rows with a valid subscript vector.
x Subscript `time` must be a simple vector, not a matrix.
Run `rlang::last_error()` to see where the error occurred.

Sample date

 structure(list(NetC = c("Cosenza Provincia", "Cosenza Provincia", 
"Cosenza Provincia", "Reti Private", "Reti Private"), Stat = c("Citta dei Ragazzi", 
"Rende", "Acri", "Firmo", "Schiavonea"), ID = c("IT1938A", "IT2086A", 
"IT2110A", "IT1766A", "IT2090A"), Longitude = c(16.2452, 16.2433, 
16.3868, 16.194, 16.5468), Latitude = c(39.3134, 39.3389, 39.4896, 
39.7138, 39.6518), Data = structure(c(1619820000, 1619820000, 
1619820000, 1619820000, 1619820000), tzone = "", class = c("POSIXct", 
"POSIXt")), PM10 = c(28.48, 29.31, 32.3, 103.06, 19.83), day = c("01", 
"01", "01", "01", "01"), month = c("05", "05", "05", "05", "05"
), year = c("2021", "2021", "2021", "2021", "2021"), avePM10 = c(17.53875, 
20.6028571428571, 19.4285714285714, 29.46625, 18.7257142857143
), comment = c(NA_character_, NA_character_, NA_character_, NA_character_, 
NA_character_)), row.names = c(NA, -5L), class = c("tbl_df", 
"tbl", "data.frame"))

Stations:

structure(list(Location = c("IT1766A", "IT1938A", "IT2086A", 
"IT2090A", "IT2110A"), Longitude = c(16.19397, 16.24517, 16.24334, 
16.54677, 16.3868), Latitude = c(39.71376, 39.3134, 39.33893, 
39.65176, 39.48963), station_type = c("Industrial", "Background", 
"Traffic", "Industrial", "Background"), station_area = c("Rural", 
"Urban", "Urban", "Rural", "Urban")), class = "data.frame", row.names = c(NA, 
-5L))

Answer 1

If you really need a STFDF object then use the STFDF function. Your data must have n x m rows where n is the number of spatial locations and m is the number of time points. m=1 seems to be invalid so I can't illustrate this with your sample data but I'll extend it for an example.

The data must also be sorted into space-time order, so all the data for the first time point, then all the data for the second time point, and so on until you have n x m rows.

Starting with pts as your sample data, I'll first convert to a plain data frame because it seems to be a tbl_df and this is incompatible with real data frames - avoid as much as you can.

pts = structure(....
pts = as.data.frame(pts)

Since this is only one point I'll duplicate it a couple of times to 15 rows and set the times of points 6-10 and 11-15 to 1 day and 2 days past the first five data points which are all at the same time point:

pts = rbind(pts, pts, pts)
pts$Data[6:10] = pts$Data[1:5] + 24*60*60 # add one day
pts$Data[11:15] = pts$Data[1:5] + 2*24*60*60

I'll also jitter the measure variable a bit too so the data values aren't repeated so we can see what's happening. Don't do this on your real data of course...

pts$PM10 = pts$PM10 + rnorm(15, 0, 2)

Now we can work out the "grid" points in space and time:

n = 5
T = unique(pts$Data)
S = SpatialPoints(pts[1:n,c("Longitude","Latitude")])

And you can confirm there are 3 time points and 5 space points, and we know they are ordered correctly. The first 5 rows are day 1, the next are day 2, the last day 3, and the locations are always in the same order. Your data must be in this format.

Now all is set up to make a STFDF:

pts.st = STFDF(S, T, pts)

which we can check and plot:

> class(pts.st)
[1] "STFDF"
attr(,"package")
[1] "spacetime"
> stplot(pts.st[,,"PM10"])