Hi,

Straightforward answer: Anisotropy ratios for geostatistics are typically calculated, although depending on implementation, by range of main over range of minor. So:

r1 = main/minor1

r2 = main/minor2 (if an ellipsoid, instead of an ellipse)

minor1 direction would be the orthogonal of the main direction (ellipse). For 3D you must give two orthogonal to build an ellipsoid (this is given by three angles of orientation: azimuth, dip, rake).

Detailed answer: If you want to make a good study about anisotropy in your case study than you should make variograms for a lot of directions. Let's say one variogram for every 10º in which the tolerance typically does not exceed the 5º.

For each of those directions you will interpret a range (given by adjusting a model). This means you can plot in a radial axis one value for every 10º. If you do this for 180º (remember directions are symmetrical) you'll get a 360º study.

Look at your plot and try to put there an ellipse (check the illustration attached). The bigger axis of that ellipse will be  your main direction. The orthogonal will be your minor1.

About experimental variograms: Experimental variograms are interpreted which means that you need to be aware of the phenomena you're studying to use as criteria for that interpretation.

Let's take as an example the variograms you've put here. You say it's PM (particulate matter?!). I'm guessing this kind of phenomena is very continuous. Using that assumption I would also guess that the 20% nugget effect you've put there does not exist. I think that the "structures" you are studying have small focus points so the variability rises very fast following a more or less constant growth from that point (check the illustration of the variogram model).

The second thing I'm noticing is the fact that you are not plotting directly the sill. Another way to put this is: you are not making an assumption about the frontier between data you consider correlated or not. You did put a gaussian model there which leads me to think that you've chosen your sill somehow (typical default is the variance of the data).

Don't be afraid to use a composite of models for your interpretations. This will bring you a lot more flexibility in interpretation.

Hope it helps.

Dear  Pedro Correia, Thanks for your detail answer. But probably for being a newbie, I couldn't understand some part. Whatever, I have a specific question.

What I know about anisotropy modeling using gstat package is, after plotting the directional variogram, I need to find the major axis (the direction of maximum range) and the ratio (range of minor axis(orthogonal to major axis)/ range of major axis) by eye estimation. [because its complex to fit the separate model for every direction and determine the range]

(My concept might be wrong)

So, I have already attached my directional plot. now I have to tell gstat, the value of major axis and ratio. My question is, according to your opinion what could be the value of major range and ratio based on my attached plot?

[the model you can see on my plot is just a general model which I fitted for isotropic variogram.you can ignore this model]

[[I am attaching the 10 degree directional plot(as you said) and my R code using gstat package]]

Hi Uzzal,

Yes. If I understood your purpose that is incorrect. If you want to perform kriging than what you need is a variogram ellipsoid that represents the anisotropy that exists in the spatial dispersion of your variable.

Checking for the direction of maximum range will only tell you that (which is insufficient). You might even want to dis-consider that direction if you see that it does not fit the global behavior of the variable.

I've never worked with gstat but from your code it seems you are giving a tolerance of 15º (I'm guessing for each side) for every 10º. You should be aware that you are reusing point from the previous and posterior direction for every one of those plots. This is ok if you know what you are doing but be careful nevertheless.

As direct answer to your question I would say, from your plots, the direction of maximum range is about 60º or 70º (notice how most of the points are below the generalist model you've plotted there). I'm saying this because a bigger range means you have less variability (so smaller semi-variogram values) in a bigger distance (so most points are at the bottom*).

Consequently the orthogonal direction should be somewhere around (60/70+90 = 150/160). But again, I think building an ellipsoid purely based on this notions will give you a very poor model indeed. Notice that the directions of, what seems to me, smaller range (about 30/40º) are not 90º from the maximum. You have something else going on there.

Finding the range or ratio by eye estimation from those plots is very tricky. You didn't plot the sill (so I'll be assuming from the converge of the gaussian you put there) or assumed the function you want to use. In any case I'm assuming something more or less like 120000 for maximum range, and 60000 for minor1 which gives:

main/minor = 120000/60000 = 2

Again my advice is that you follow the guidelines I've given in my first post. Making these kind of studies like this can get you to to completely miss the behavior of your variable.

* by the way, I'm assuming these are all variogram values (so no covariances or correlograms). If you don't know what this is don't worry.

Coreria has given you good advice but several points/issues mentioned deserve additional emphasis. First of all gstat is not a single program, i.e. it has several components. In the kriging component you must specify the maximum and minimum range, i.e. the ratio. In contrast in the component to compute/plot sample variograms you must choose a direction (directions if in 3-D) and an angular tolerance. If you are computing multiple directional sample variograms, e.g. with 10 degree increments then you want the angular tolerance to be 5 degrees or less. Otherwise you will have some paired differences being used to compute more than one directional sample variogram, it will be very difficult to determine the effect of doing that so it is best to avoid it. Note also the total number of distinct paired differences is  n(n-1)/2 where n is the number of data location points, this means that the more different directional sample variograms you compute the fewer paired differences per direction and the fewer paired differences means it will be harder to interpret the plot of the directional sample variograms. I would suggest starting with directional sample variograms at 0, 45, 90,135 degrees and angular tolerances of 22.5 degrees.The objective here is to ask whether there is evidence of a (geometric) anisotropy, i.e. are the apparent ranges different for the different directions. Depending of the size of your data set it may be that none of the directional sample variograms is very good, i.e. it is hard to ascertain the shapes and the ranges. In that case you would assume there is no anisotropy (only one range and the ratio is 1)..

Note also that a geometric anisotropy only allows the range to change with the direction (not the sill).

Correia also suggests that you consider the particular phenomena generating the data and what you know about it, this is good advice. One additional way you can do this is to make a coded plot of the data locations. Each data location coded according to the value at that location, e.g. different colors. It is possible that a directional affect would show up on such a plot. Since you are using R. there is another package that works well with R called "ggobi"  also free and comes with an electronic textbook, it has very nice graphics. Just do a search on Google to find ggobi..

Finally estimating/modeling a variogram from data is not a "black box" process, you likely need to do some experimenting (trial and error) including the use of cross validation,e.g how do the cross validation statistics change if you use a model with an anisotropy and with different choices of the main direction and ratio, how sensitive are the statistics to the changes?

Bottom line, are the results useful, do they make sense given your understanding of the phenomena? Geostatistics is based on some assumptions but in general those assumptions are not fully (statistically) testable. The software will give you results based on what you put in (data and all the various choices) but it won't tell you whether the results are "correct/truet"

Hello

in my opinion you do not have anisotropy in your study.

 So you can use the SURPAC Packag

Hello

in my opinion you do not have anisotropy in your study.

 So you can use the SURPAC Packag

You don't mention doing any exploratory analysis of the data prior to trying to fit a variogram model. One of the basic underlying assumptions is intrinsic stationarity and in particular a constant mean. While it is not areal statistical test, you should look at a trend surface of the PM data, it is possible that the anisotropy you perceive in the directional empirical variograms is actually a directional effect in a non-constant mean. This could be a physical effect of wind velocity.

Your comments seem to suggest that you are using the documentation for gstat to understand geostatistics, that is not adequate. The documentation is only intended to tell you how to use the software, i.e. the mechanics of using the software It is not an adequate introduction to the theory. You need to consult one or more of the various textbooks on geostatistics. Two of the the best are "Mining Geostatistics" by A. Journel and Ch Huijbregts (Academic Press) and "Geostatistics:modeling spatial uncertainity", J.P. Chiles and P. Delfiner (J. Wiley). There is also a very short book called "Basic Linear Geostatistics", M. Armstrong.

There are case studies in the literature on modeling anisotropic variograms, one by A. Journel in Mathematical Geology, sometime in the early 1980's and a series of papers by R. Webster in the J. of Soil Science (now known as the European J. of Soil Science) that appeared in 1979-1980.

Unless you have a very large data set it is likely that the directional empirical variograms will not be very useful (too few paired differences for each lag distance) especially if you compute the empirical variograms for each 10 degrees using a tolerance of 5 degrees. Start out using every 45 degrees and a tolerance of 22.5 degrees.

In order to model the anisotropy, I personally use "anisotropy map". A very useful tool to identify the direction of maximum continuity (=long range).

"Anisotropy maps" are provided both in Surpac 6.2 (and newer versions) and GS+.

Here I put one of my sample works as the theory has explained by experts.

Like Hossein Molayemat said, you could try plotting variographic maps (e.g. with the varmap program from GSLib).

Uzzal

You didn't say how many data locations you have. This is important because for n data locations you will get n(n-1)/2 possible pairings. When you compute directional experimental variograms you will be splitting the pairings up between the different directions, you may end up with too few pairs for some or all of the directional experimental variograms and thus find it more difficult to select a model type. If there is a "trend" in the data that likely is directional and might be confused with an anisotropy. You are likely going to have to do some experimenting with the anisotropy parameters and then apply cross validation to see how sensitive the results are to different choices of the parameters and perhaps to the variogram model type. You are looking at air pollution data, what do you know about prevailing wind directions? The pattern of the data locations can result in seeming to need an anisotropy when it is not true.

I suggest you look at the book "Geostatistics" by Chiles and Delfiner.

A final comment, there is no assumption of normality for kriging, i.e. in the derivation of the kriging equations.

Note that for a geometric anisotropy the model type must be the same for all directions. Only the "range" changes with direction for a geometric anisotropy

Note also that " range has several different meanings

1. The distance at which the variogram value is a constant, e.g. a Spherical model

2. The distance at which the variogram value is 95% of the asymptotic value, e.g Exponential model and Gaussian model. 95% is an arbitrary choice but seems to be standard in the literature. Sometimes this range is called the "effective range" since neither the Exponential nor Gaussian models have a"true range" in the sense of 1.

Well written software will use the same name for both of these but the evidence on the plot will be slightly different.

Pedro suggests using the sample variance for the default value of the sill, this is not correct. The sill is the variance of the random function but the sample variance is not a good estimator since the data is not a random sample (in a geostatistics application). Random selection of the data locations is not random sampling.

Several respondents have suggested using Surpac. this is commercial software and the code is not accessible, as contrasted with gstat. GSLIB is a collection of FORTRAN codes and GS+ is essentially based on GSLIB.

Were the recordings all taken at the same time? If not you need to consider space time models.

If you want, you could try the software I've developed for my PhD experiments: https://github.com/PauloCarvalhoRJ/gammaray/releases . It is a GUI over GSLib (separate install) and it is free and open source. You can get vanilla GSLib executables for your OS here: http://www.statios.com/Quick/gslib.html

I forgot to mention that you will also need GhostScript to see GSLib plots (either in the program or by directly using GSLib executables): https://www.ghostscript.com/download/gsdnld.html .

If your data set has too many points, it may be interesting to grid it using the nearest neighbor estimation and compute the variographic map fast with FFT. Both steps you can do with the software I mentioned. The program's manual presents the theory in case you'd like to obtain a variographic map using FFT elsewhere (e.g. in Matlab, Octave or Python). Traditional variogram computing is O(n^2) while FFT reduces this figure to O(n log n). In addition to being fast, computing the variographic map with FFT is parameter-less, so no tiresome task of configuring lags, tolerances, bandwidths, etc. by trial-and-error.

In this notebook I demonstrate how you can compute variogram maps with GeoStats.jl: https://nbviewer.jupyter.org/github/juliohm/GeoStatsTutorials/blob/master/notebooks/DirectionalVariograms.ipynb