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

[Edzer Pebesma]

On 12/09/2012 11:03 PM, Michele Di Marcantonio wrote:
> Dear Edezer,
> 
> thank you for your help and for all your suggestions.
> 
> I have tried to develop the analyses adding the "$dist > D" parameter in
> order to throw out records but I have got an error message, surely for my
> fault. Anyway, if possible, I would like to calculate the empirical
> variogram just for lag distances higher than a "D" value without throwing
> away records.
> 
> I enclose you a txt file "DATA" with [x, y, Z] sample values, where (x, y)
> are the coordinates of the Z(x, y) observations.
> Below is the base for a possible code in gstat to get what I need following
> those sample data.
> 
> I would be grateful if somebody could help me showing me directly how
> should I integrate and/or correct that code in order to:
> 
> 1) get a variogram model just for lag distances higher than a "D" level
> (for the sample: e.g. for lag distances between 0.24 and 0.42, so
> cutoff=0.42 and minimum distance (function?)=0.24);

variog1=variog1[variog1$dist > 0.24 & variog1$dist < 0.42,]
fit.variog1 <- fit.variogram(variog1, model.variog1)
plot(variog1,fit.variog1)

selects the sample variogram values in the range you mention, and fits a
model to those. Of course, this model again is valid for every distance,
as a variogram model is supposed to.

> 2) calculate a 3D anisotropic variogram (both empirical and a model) which
> can be used for kriging (krige(z~1, DATA, LOCATIONS, model=fit.variog1)).

I gave suggestions in my previous answer; did you try them?

> 
> Thank you again,
> 
> Michele
> 
> ---
> 
> ### Libraries
> library(gstat)
> ### Data
> DATA<- read.table("DATA.txt",header=TRUE)
> names(DATA)
> attach(DATA)
> ### Empiricial variogram
> variog1 <- variogram(z~1, locations=~x+y, DATA, cutoff=0.42)
> variog1
> plot(variog1)
> ### Variogram fitting
> model.variog1 <- vgm(psill=0.7, model="Gau", nugget=0, range=1)
> fit.variog1 <- fit.variogram(variog1, model.variog1, fit.sills=TRUE,
> fit.ranges=TRUE)
> plot(variog1, model=fit.variog1)
> fit.variog1
> 
> 
> 2012/12/6 Edzer Pebesma <edzer.pebesma at uni-muenster.de>
> 
>>
>>
>> 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
>>
>> _______________________________________________
>> 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