[R-sig-Geo] How to develop variogram models in R: setting minumum lag distance and modelling 3D anisotropy

[Edzer Pebesma]

On 12/06/2012 01:02 PM, Michele Di Marcantonio wrote:
> Dear All,
> 
> 
> I have some questions about the development of variogram models in R, in
> paerticular for setting a minimum lag distance and for modelling 3D
> anisotropy.
> 
> 
> The aim of our research is to analyze a set of spatial data Z(x; y)
> distributed on a regular grid (x; y). We want to develop a Kriging
> forecasting model trying different variogram models (both isotropic and
> anisotropic) fitted using sample observations (empirical variograms).
> 
> If we use gstat (NOT STRICTLY NECESSARY) the initial codes to be run in R
> are the following:
> 
> 
> library(gstat)
> DATA<- read.table("DATA.txt",header=TRUE)
> 
> 
> In the following questions we specify some R (gstat) functions and
> parameters that we consider appropriate for reaching our objectives of
> analysis - although obviously still to be integrated/corrected considering
> the answers to our questions. However, in case you consider more
> appropriate to use different libraries, functions or parameters, we would
> be grateful if you could specify the codes we should use to run the
> analyses in R.
> 
> 
> 
> 
> 
> *1) *In the <variogram> function, for calculating the empirical variogram
> (we call it “V.EMP”), we can set the [cutoff], that is the maximum lag
> distance between the observations to be considered. For example, for a
> maximum distance equal to 0.07:
> 
> 
> 
> V.EMP<- variogram(z~1, locations=~x+y, DATA, cutoff=0.07)
> 
> 
> *How should we do to set also a minimum distance below which observations
> must not be considered to determine the variogram? Which parameter should
> we add in the <variogram> function (if that function is the correct one)?*
> 
> *Or should we add a minimum distance parameter later in the <vgm> function
> for estimating the variogram model (and in case which one)?*

V.EMP is a data.frame, so if you want you can throw out records e.g. by

V.EMP <- V.EMP[V.EMP$dist > 0.01,]

> 
> The empirical variogram is used for fitting a variogram model and finally
> for Kriging. For example:
> 
> *
> *model.variog <- vgm(psill=0.00725878, model="Exp", nugget=0,
> range=0.2326746)
> 
> fit.variog <- fit.variogram(V.EMP, model.variog, fit.sills=TRUE,
> fit.ranges=TRUE)
> 
> OK <- krige(z~1, DATA, FORECASTS.COORDINATES, model=fit.variog)
> 
> 
> *2)* *How should we do to calculate the anisotropic empirical variogram
> (that is the empirical variogram in the three dimensions x, y and z)?*

look into arguments alpha and beta in ?variogram, as well as tol.hor and
tol.ver
> 
> *How should we do to fit a 3D anisotropic variogram model (estimating also
> the anisotropy parameters) to be used for Kriging?*

There is no high-level interface for this; you can use optim to do this,
but it's a bit of work; I have no example lying around.

> 
> *What are the R codes for calculating i) the Ordinary Kriging and ii) the
> Universal Kriging using that anisotropic variogram model (if gstat is
> appropriate, can we use the same <krige> function?)?
> 
> *

Yes.

> 
> We have tried to add one by one the following parameters (assuming an
> anisotropy angle of 90°) in the <variogram> function: [dir=pi/2],
> [dir.hor=90], [alpha=90]; nevertheless, those insertions produced no
> effects, as all the resulting variograms we obtained coincide with the one
> calculated without adding any of those parameters.
> 
>  Similarly, also the insertion of the [anis=c(k1, k2, k3, k4, k5)]
> parameter in the <vgm> function (with realistic values of the ks) produced
> no effects.
> 

If you find unexpected results, it might help to present a reproducible
example to the list.

> 
> model.variog <- vgm(psill=0.00725878, model="Exp", nugget=0,
> range=0.2326746, *anis=c(k1, k2, k3, k4, k5)*)
> 
> 
> 
> Furthermore, even if it produced effects, considering that – if we have
> understood well – the <fit.variog> function does not estimate the
> anisotropy parameters (the ks), we would have still the problem of
> estimating them.
> 
> 
> Once more, we have given suggestions for possible functions in gstat that
> we consider appropriate, but in case you suggest other libraries or
> functions, consider that we just need to know a code in R to calculate: i)
> empirical variogram, ii) variogram model (with minimum lag, both isotropic
> and anisotropic), iii) fit the variogram model and iv) develop Kriging on
> the base of the variogram model.
> 

Thank you for your suggestions.

I suggest using optim, then writing a well documented high-level wrapper
function, and contributing this back to the developers of gstat or some
other geostats package where this would be useful, so they can integrate
it for future users with similar needs.

> 
> 
> Thanks in advance for your help
> 
> 
> 
> Best regards,
> 
>  Michele Di Marcantonio
> 
> 	[[alternative HTML version deleted]]
> 
> 
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de