Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

wrong prediction with krigeST() #112

Open
zheangzhao98 opened this issue Aug 25, 2022 · 0 comments
Open

wrong prediction with krigeST() #112

zheangzhao98 opened this issue Aug 25, 2022 · 0 comments

Comments

@zheangzhao98
Copy link

I have no idea what is going on. I use data in Jan to fill the missing data in 27th Jan
My data:
result.csv

My code:

library(tidyverse) 
library(sf) 
library(sp) 
library(raster) 
library(gstat)  
library(automap) 
library(lattice)
library(patchwork)
library(viridis)
library(spacetime) 
library(ggplot2) 
library(dplyr) 
library(gstat) 
library(tidyr)

data_tidy.subset <- read.csv("result.csv")
coordinates(data_tidy.subset) <- ~lat.zone + lon.zone

sp = unique(data_tidy.subset@coords)
sp = SpatialPoints(sp)
time = as.Date("2003-12-31")+c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,20,21,22,23,25,27,29,30,31)
data = as.data.frame(cbind(data_tidy.subset@coords,rep(2004,21600),rep(1,21600),data_tidy.subset@data$time.period,data_tidy.subset@data[["SST.zone.period"]]))
colnames(data) <- c("lon","lat","year","month","time.period","SST.zone.period")
data.STFDF <- STFDF(sp,time,data)
data_tidy.subset = na.omit(data_tidy.subset)
data_tidy.subset.vgm <- variogram(object = SST.zone.period ~ 1 + lon + lat,
                                  data = data.STFDF,
                                  width = 0.25, cutoff = 2)
plot(data_tidy.subset.vgm,wireframe=T)

sepVgm <- vgmST(stModel = "separable",
                space = vgm(10, "Exp", 5, nugget = 0),
                time = vgm(10, "Gau", 50, nugget = 0),
                sill = 20)
sepVgm <- fit.StVariogram(data_tidy.subset.vgm, sepVgm)

metricVgm <- vgmST(stModel = "metric",
                   joint = vgm(100, "Exp", 100, nugget = 0),
                   sill = 10,
                   stAni = 0.4)
metricVgm <- fit.StVariogram(data_tidy.subset.vgm, metricVgm)

prodSumVgm <- vgmST("productSum",
                    space = vgm(10, "Exp", 5, nugget = 0),
                    time = vgm(10, "Gau", 50, nugget = 0),k = 11)
prodSumVgm <- fit.StVariogram(data_tidy.subset.vgm, prodSumVgm)

sumMetricVgm <- vgmST("sumMetric", 
                      space = vgm(10, "Exp", 5, nugget = 0),
                      time = vgm(10, "Gau", 50, nugget = 0), 
                      joint = vgm(100, "Exp", 100, nugget = 0), stAni=0.4) 
sumMetricVgm <- fit.StVariogram(data_tidy.subset.vgm, sumMetricVgm)

attr(sepVgm, "optim")$value
attr(metricVgm, "optim")$value
attr(prodSumVgm, "optim")$value
attr(sumMetricVgm, "optim")$value


plot(data_tidy.subset.vgm, 
     list(sepVgm,metricVgm,prodSumVgm,sumMetricVgm), 
     main = "Semi-variance",all=T,wireframe=T)

###prediction st kriging
spat_pred_grid <- expand.grid(lon=seq(20.125, 29.875, 0.25), 
                              lat=seq(-44.875,-40.125, 0.25))%>%
  SpatialPoints(proj4string = CRS(proj4string(data.STFDF)))
gridded(spat_pred_grid) <- TRUE
temp_pred_grid <- as.Date("2003-12-31") + c(26,27)
DE_pred <- STF(sp = spat_pred_grid, # spatial part
               time = temp_pred_grid) # temporal part
data.STIDF <- as(data.STFDF, "STIDF") # convert to STIDF
data.STIDF <- subset(data.STIDF, !is.na(data.STIDF$SST.zone.period))
pred_kriged1 <- krigeST(formula = SST.zone.period ~ 1  + lon + lat, data = data.STIDF, 
                        newdata = DE_pred, modelList = sumMetricVgm, computeVar = TRUE)
View(pred_kriged1)
stplot(pred_kriged1)
# for free to join this conversation on GitHub. Already have an account? # to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant