0

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.

enter image description here

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?

Community
  • 1
  • 1
Orpheus
  • 329
  • 6
  • 21
  • 2
    Because your question is rather about choosing and interpreting statistical methods than about programming, I believe http://stats.stackexchange.com/ is a more suitable site. I have voted to transfer it. – Henrik Jul 28 '15 at 17:49

1 Answers1

2

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))

Level Plot (directional semivariograms)

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))

Directional semivariograms

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))
Jared Smith
  • 280
  • 3
  • 12
  • How can I understand from above 7 directional variogram that there is anisotropy or not? I know about geometric and zonal anisotropy. But can I say from the above figures that sill is nearly same in all direction? and also I cant understand that range in different or not in these figures! Could you please comment on the anisotropy based on this figure? – Orpheus Jul 29 '15 at 07:15
  • What you're really looking for is which direction has the shortest `range` to the `sill`, and which direction has the longest. From the level plot and directional variograms it looks like 30 degrees could be the minor axis (has some increase in semivariance, potentially to a sill), which would make 120 degrees the major axis. I recommend plotting those angles separately and seeing what those variograms look like. – Jared Smith Jul 29 '15 at 11:52
  • Thanks. But I saw some example where they show the numeric value of sill and range in table for every direction. By which I can be clear about the major and minor axis. Link: http://www.techmat.vgtu.lt/~art/proc/file/BudrLi.pdf (page 365). Is there any built function in R by which I can see the numeric value of range and sill for every direction? – Orpheus Jul 29 '15 at 13:30
  • If the `sills` were reached in all directions, then yes, you would be able to make a table like the one you reference by using `fit.variogram` on a `vgm()` defined for the `variogram()` with the `alpha` of interest (with a reasonable `tol.hor` on `alpha`). Then use something like `sill = sum(model.3112$psill)` and `range = model.3112$range` to extract them. – Jared Smith Jul 29 '15 at 16:03
  • Did you mean this? `seoul3112.var1<-variogram(g,width=8000,cutoff=80000, alpha=seq(0,180,45),tol.hor=15)` // `model<- fit.variogram(seoul3112.var1,vgm(800,"Sph",40000,200),fit.method = 2)` // `sill = sum(model$psill)` // `range = model$range` // `sill` // `range` I tried this. but didn't got my desired out put. (here. object g which I created in my question) – Orpheus Jul 29 '15 at 17:16
  • I meant using something like `seoul3112.var1<-variogram(g,width=8000,cutoff=80000, alpha=MajorAngle,tol.hor=15)` first so that a variogram for only the major direction of anisotropy is created. Then repeat for the minor angle (or any other angle) and compare. – Jared Smith Jul 29 '15 at 17:36