What KInd of variogram should I plot before doing Kriging in R?

Question

I have a csv file named seoul3112, contains PM10 concentrations. please, download..I have tried to plot sample variogram and fit a model on it.

library(sp)
    library(gstat)
    library(rgdal)
    seoul3112<-read.csv("seoul3112.csv",row.name=1)
    seoul3112<-na.omit(seoul3112)

#assign a CRS and reproject

coordinates(seoul3112)=~LON+LAT
proj4string(seoul3112) =  "+proj=longlat +datum=WGS84" 
seoul3112<-spTransform(seoul3112, CRS("+proj=utm +north +zone=52 +datum=WGS84"))

#plot semi-variogram

g<-gstat(id="PM10",formula=PM10~LON+LAT, data=seoul3112)
seoul3112.var<-variogram(g, cutoff=70000, width=6000)
seoul3112.var
plot(seoul3112.var, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 3112")

#Model fit

model.3112<- fit.variogram(seoul3112.var,vgm(700,"Gau",40000,400),fit.method = 2)
                       
plot(seoul3112.var,model=model.3112, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 3112")

After writing these code I got a semi-variogram like this.

As I am very new in geostatistics, so I am confused that my above variogram is okay or not for my dataset. Because, in typical variogram, semivariance value become horizontal at sill. But this variogram going upward! Should I make make some correction in my code?

Another thing is, actually my final goal is doing kriging interpolation on my dataset(seoul3112).I don't understand that, for doing kriging, this sample variogram is enough or should I plot directional variogram or any other thing? Could anyone please explain it in details?

Answer 1

If you look at your fitted model, it will have a sill parameter (i.e. nugget + psill), but it is beyond the extent of your samples.

 sill = sum(model.3112$psill)

It is possible that you did not have a point-pair distance that was far enough to arrive at the range to the sill. I do not think this is not a problem, so long as you are using this semivariogram to predict within the spatial domain of your data + 60000 m (the distance for which you have data support). I take this position because the most important portion of a semivariogram to fit is the beginning, and in the extent of the data the fit is good.

Something to check is how many point pairs (np) you have supporting the bins that you plotted (more point-pairs is better).

Another suggestion is looking for anisotropy using a level plot

seoul3112.var_map<-variogram(g, cutoff=70000, width=6000, map=TRUE)
plot(seoul3112.var_map, col.regions=terrain.colors(20))

or by looking at the fit in many directions by setting alpha and plotting the model fit. You only need to examine 180 degrees because it's symmetric (you can see this in the plot below, in which 0 and 180 are the same).

seoul3112.var_angles<-variogram(g, cutoff=70000, width=6000, alpha=seq(0,180,30))
plot(seoul3112.var_angles, model=model.3112, pch=16, ylim=c(0,3000))

If there are differences directionally, you can model the major and minor axes of anisotropy using fit.variogram with anis defined for the vgm() model. Example:

vgm(700,"Gau",40000,400, anis=c(someangle, someratio))