# How to remove zonal anisotropy from directional variogram using gstat in R?

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: 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.

• possible duplicate of stackoverflow.com/questions/31778893/… – user32309 Aug 6 '15 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 Aug 6 '15 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 Aug 6 '15 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 Aug 10 '15 at 15:08

## 1 Answer

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)
``````
• It's working. However, could explain why we need to set very large ranges? How to decide the proportions, etc.? – Jot eN Feb 27 '16 at 12:57
• 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 Feb 28 '16 at 16:15