rspatialkriginggstatcovariogram

How to fit model with semivariogram using gstat in R?


I have a csv file file contains atmospheric PM10 concentration data of 1 march,12.00 pm.please, download. I want to draw a semivariogram using gstat package in R. I tried to write these code in R. but with these data, I cant fit the model.

    library(sp)
    library(gstat)

    seoul3112<-read.csv("seoul3112.csv")
    seoul3112<-na.omit(seoul3112)

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

model.3112<- fit.variogram(seoul3112.var,vgm(700,"Gau",0.5,200), 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")

Actually I am beginner in R and statistics. So, I am very ignorant even about variogram. I have some query :

a)When I plot the data as a semivariogram, it looks different not as typical semivariogram! why this is happening? should I do any other thing with my data like transforming?

b)How Can I fit the model with this data? I have tried different model like "Sph","Exp" but they look like linear! why?

c)How can I understand that what initial value of sill,range,nugget I should use in vgm() function?

d)How can I understand that model fits with data properly?

e)For using kriging, what kind of semivariogram I should plot? only Omnidirectional semivariogram? or I should plot directional semivariogram?

f)And how can I interpret the semivariogram? I mean what actually I can understand about the data from semivariogram?

Thanks in advance.


Solution

  • I'll provide answers to your code related questions. The rest of your questions (d, e, and f) are more theory related.

    First, in your comment, when you changed the proj4string, the distance units should have changed on the plot. Did they? Based on your comment, it sounds like that did not happen.

    a) In addition to playing around with the cutoff distance, also be careful of the np (point-pairs) that you have supporting each bin on the semivariogram. For example, using your updated proj4string information, I tried cutoff=80 and width=80/10 (10 bins instead of 15) to see how the semivariogram shape changed. Decreasing from 15 to 10 bins does not change where point-pairs exist, just increases the distance that each bin represents. Also, this approach is not necessarily what you should use, but it's an example of how to change the bins for smoother sample semivariograms (but smoother does not mean better).

    Comparison of 10 Bins to 15 Bins

    b) Using your code, the "Sph" and "Exp" models return Warning: singular model in variogram fit. This warning indicates that there is not enough data to fit some parameters of the spherical and exponential empirical models. See the gstat user manual for guidance on each of the empirical equations and their parameters.

    c) The vgm() function can be used, for example, to eye-fit sample semivariograms. If you're confused about how to plot the vgm() model with the sample data, try something like

    eye_vgm = vgm(psill=1200,model="Gau",range=60,nugget=350)
    plot(seoul3112.var,model=eye_vgm, col="black", pch=16,cex=1.3)
    

    You're using vgm() in a call to fit.variogram(), so as long as the parameters you give to vgm() are reasonable (e.g. based on the sample data) and the empirical model can have parameters fitted, fit.variogram() will find a fit according to the fit.method.