1

I have a dataframe named seoul315 contains the 107 data of PM10 concentration at 1 march,2012,5.00 am. please, download. I have tried to plot semi-variogram for kriging bt the following code:

library(sp)
library(gstat)
library(rgdal)

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

hist(seoul315$PM10)   #Data are not in Gaussian form.So need to log transform
hist(log(seoul315$PM10))

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

#plot Omnidirectional Variogram
seoul315.var<-variogram(log(PM10)~1,data=seoul315,cutoff=80000, width=10000)
seoul315.var
plot(seoul315.var, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 315")

#Model fit
model.315<- fit.variogram(seoul315.var,vgm(psill=0.11,model="Gau",range=60000/sqrt(3),nugget=0.04),
                          fit.method = 6)
plot(seoul315.var,model=model.315, col="black", pch=16,cex=1.3,
     xlab="Distance",ylab="Semivariance",
     main="Omnidirectional Variogram for seoul 315")

#Directional Variogram
seoul315.var1<-variogram(log(PM10)~1,data=seoul315,width=10000,cutoff=80000,
                         alpha=seq(0,135,45),tol.hor=15)
plot(seoul315.var1,model=model.315, cex=1.1,pch=16,col=1,
     main="Anisotropic Variogram for PM10")

I got the directional variogram like this: enter image description here

But from this variogram I can see sill are not same in all direction that means its zonal anisotropy. I think in this case I can't use anis() function which I used to remove the geometric anisotropy. How can I remove the zonal anisotropy for using the variogram in kriging? Is it possible by using gstat in R?

Actually I can't find the proper way for doing this after rigorous searching. If you need further information please let me know.

Improve this question
4
  • possible duplicate of stackoverflow.com/questions/31778893/…
    – user32309
    Commented Aug 6, 2015 at 5:26
  • Sorry, I posted this question on SO but I didn't get response and someone suggested me this forum. Could you please let me know that this type of problem is appropriate for this page or not? Actually, I badly struck with this problem.
    – Orpheus
    Commented Aug 6, 2015 at 6:49
  • The problem is it looks more than a (geo)statistical issue than a programming issue. Let's see what happen, as it is correct that here is a better place than in SO.
    – user32309
    Commented Aug 6, 2015 at 6:58
  • This anisotropy appears to result from a trend, so the first thing to consider is either removing an underlying trend or using a method than incorporates trend fitting (such as Universal Kriging).
    – whuber
    Commented Aug 10, 2015 at 15:08

1 Answer 1

Reset to default
2

Zonal anisotropy is not something you can remove, but you can try to model it. R package gstat lets you model it by geometrically anisotropic models with very large ranges; try

library(gstat)
demo(zonal)
Improve this answer
2
  • It's working. However, could explain why we need to set very large ranges? How to decide the proportions, etc.?
    – Jot eN
    Commented Feb 27, 2016 at 12:57
  • 1
    We do that to emulate zonal anisotropic models with geometrically anisotropic ones. You set the range so large that effectively the component is absent in one direction, then the ratio such that it is present in the perpendicular direction.
    – Edzer Pebesma
    Commented Feb 28, 2016 at 16:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.