[R-sig-Geo] Spatio Temporal kriging in Gstat

Dan Turenne DanielTurenne at hotmail.com
Thu Jun 30 01:42:32 CEST 2016


Thank you very much for the tip Edzer, that was exactly what I was looking for.


Daniel Turenne

________________________________
From: Edzer Pebesma <edzer.pebesma at uni-muenster.de>
Sent: June 29, 2016 6:17 PM
To: Dan Turenne
Subject: Re: [R-sig-Geo] Spatio Temporal kriging in Gstat



On 29/06/16 11:52, Dan Turenne wrote:
> Hello Ben,
>
> Thank you for the tip about "tunit", it worked perfectly as you described.  As for my second question, the object returned by krigeST is indeed of the form STFDF.  After some digging I found the following information regarding the data from an STFDF object on inside-r.org:
>
> data<http://inside-r.org/r-doc/utils/data>:
data {utils} | inside-R | A Community Site for R<http://inside-r.org/r-doc/utils/data>
inside-r.org
Arguments... literal character strings or names. list a character vector. package a character vector giving the package(s) to look in for data sets, or NULL.



> Object of class data.frame<http://inside-r.org/r-doc/base/data.frame>, which holds the measured values; space index cycling first, time order preservedI would interpret this as follows: if I had x prediction locations, the first set of x values from the predictions vector will be for day 1, the second set of x values will be for day 2 etc.  The order of the stations will be the same as they were in the sp object used to create the prediction grid.
>
> Is this interpretation correct or am I missing something?  Once again your help is greatly appreciated.

Probably, but you might also want to use

as(object, "data.frame")

and try to understand what you get.

>
> Thanks,
> Daniel Turenne
> ________________________________
> From: R-sig-Geo <r-sig-geo-bounces at r-project.org> on behalf of Benedikt Graeler <benedikt.graeler at rub.de>
> Sent: June 29, 2016 3:04 AM
> To: r-sig-geo at r-project.org
> Subject: Re: [R-sig-Geo] Spatio Temporal kriging in Gstat
>
> Dear Daniel,
>
> On 28.06.2016 16:42, Dan Turenne wrote:
>> Thank you for the information, I was able to get around the problems
>> I was having by setting upper to Inf and lower to 0.0000001 (An
>> arbitrary small value to guarantee that the value will be greater
>> then 0).  After working some more on my code I had just a few other
>> general questions I was hoping someone could answer.
> Note, that you can specify a boundary for each parameter. The arguments
> upper and lower are recycled up to the length where they match the
> parameter vector length. "Inf" is fine as a bound, but narrowing down
> the parameter space also might improve/speed up the estimation.
>
>> 1.  Whenever I call variogramST to calculate the sample variogram I
>> get a bunch of repeating output that says $tunit [1] "days" over and
>> over again.  I have tried setting progress to FALSE but this output
>> still occurs.  Is there anyway to suppress this output?
> The argument "tunit" is only necessary if a STIDF is provided. In your
> script, you create a STSDF and a different routine is used in the
> background (assuming "tlags" refers to index differences). All arguments
> that are not recognized by "variogramST" are passed on to the space only
> function "variogram", that does not know how to handle it.
> Solution: just do not provide the argument "tunit" or provide a STIDF
> where a "tunit" is necessary.
>
>> 2.  When using krigeST to get predictions, the predictions are
>> returned in the form of a single column, however there is no index to
>> say which prediction corresponds to which location.  Is this column
>> organized with the spatial index moving faster or the temporal index
>> moving faster?
> The column is supposed to be in the same order as the target geometry.
> Which version of gstat are you using? Unless "fullCovariance=TRUE", the
> returned object should be of the form ST*DF.
>
> HTH,
>
>   Ben
>
>>
>> Thank you in advance for any advice, everyone on this list has been
>> very friendly and helpful and it is greatly appreciated.
>>
>> Thanks, Daniel Turenne ________________________________ From:
>> R-sig-Geo <r-sig-geo-bounces at r-project.org> on behalf of Benedikt
>> Graeler <benedikt.graeler at rub.de> Sent: June 28, 2016 5:23 AM To:
>> r-sig-geo at r-project.org Subject: Re: [R-sig-Geo] Spatio Temporal
>> kriging in Gstat
>>
>> Dear Daniel,
>>
>> the errors occur during the optimisation where optim tries different
>> (possibly non-positive) values to find the best fit. Providing
>> sensitive limits via "lower" and "upper" to optim should resolve the
>> issue.
>>
>> HTH,
>>
>> Ben
>>
>> On 26.06.2016 16:30, Dan Turenne wrote:
>>> Hello Listers,
>>>
>>>
>>> https://www.dropbox.com/sh/7ccy5wu68gxf6sf/AACmk0AFxQvetunq4YPc9DXla?dl=0
> [https://cf.dropboxstatic.com/static/images/icons128/folder_dropbox.png]<https://www.dropbox.com/sh/7ccy5wu68gxf6sf/AACmk0AFxQvetunq4YPc9DXla?dl=0>
>
> Spatio temporal kriging<https://www.dropbox.com/sh/7ccy5wu68gxf6sf/AACmk0AFxQvetunq4YPc9DXla?dl=0>
> www.dropbox.com<http://www.dropbox.com>
> Shared with Dropbox
>
>
>
>>>
>>>
>>>
>>>
> Thanks to the help of Dr. Graeler I have been having some success with
> my spatio temporal kriging script but I had some more questions I was
> hoping someone could help with.  When fitting the ST variogram with
> fit.stVariogram, I have noticed that the function is extremely sensitive
> to the choice of initial parameters.  I have included a link to my
> script at the beginning of this email.  If I use
>>>
>>>
>>> prior.stvgm=vgmST(stModel="sumMetric",
>>> space=vgm(0.1,"Exp",0.5,0.3), time=vgm(0.1,"Exp",0.5,0.3),
>>> joint=vgm(0.1,"Exp",0.5,0.3), stAni=0.5)
>>>
>>> fitted.stvgm=fit.StVariogram(sample.stvgm,prior.stvgm)
>>>
>>> fit.stVariogram works without producing any errors, however if I
>>> adjust any of these parameters I get some strange errors.  For
>>> example by using
>>>
>>> prior.stvgm=vgmST(stModel="sumMetric",
>>> space=vgm(0.1,"Exp",0.5,0.3), time=vgm(0.1,"Exp",0.5,0.3),
>>> joint=vgm(0.1,"Exp",0.5,0.3), stAni=1)
>>>
>>> Then I get this error, despite the fact that stAni is clearly
>>> positive:
>>>
>>> Error in vgmST("sumMetric", space = vgm(par[1],
>>> as.character(model$space$model[2]),  : "stAni" must be positive.
>>>
>>> As another example if I enter
>>>
>>> prior.stvgm=vgmST(stModel="sumMetric",
>>> space=vgm(0.1,"Exp",0.5,0.3), time=vgm(0.1,"Exp",0.5,0.3),
>>> joint=vgm(0.1,"Exp",0.5,0.3), stAni=0.99)
>>>
>>> It will produce the following error, despite the fact that all
>>> ranges are clearly positive:
>>>
>>> Error in vgm(par[1], as.character(model$space$model[2]), par[2],
>>> par[3],  : range should be positive
>>>
>>> The only thing that changes among these examples is the value of
>>> the stAni parameter.  When stAni=0.5 everything is fine, if stAni=1
>>> then apparently stAni isn't positive and if stAni=0.99 then
>>> apparently the range isn't positive.  When I look at the plot of
>>> the sample variogram from my data, it appears that the spatial sill
>>> =1, spatial range=0.5, temporal sill=0.85, temporal range=4, joint
>>> sill=1.1 ,joint range=1.0 however these values also produce an
>>> error saying that range should be positive.  I'm finding these
>>> errors to be very confusing and would appreciate any tips anyone
>>> might be able to give.
>>>
>>> Thank you,
>>>
>>> Daniel Turenne University of Manitoba
>>>
>>>
>>> ________________________________ From: R-sig-Geo
>>> <r-sig-geo-bounces at r-project.org> on behalf of Benedikt Graeler
>>> <benedikt.graeler at rub.de> Sent: June 23, 2016 2:54 AM To:
>>> r-sig-geo at r-project.org Subject: Re: [R-sig-Geo] Spatio Temporal
>>> kriging in Gstat
>>>
>>> Dear Dan,
>>>
>>> the error occurs because your prior.stvgm is far from the structure
>>> of your data. The following model provides better starting values:
>>>
>>> prior.stvgm=vgmST(stModel="sumMetric",
>>> space=vgm(0.1,"Exp",0.5,0.3), time=vgm(0.1,"Exp",0.5,0.3),
>>> joint=vgm(0.1,"Exp",0.5,0.3), stAni=1)
>>>
>>> To get a first "estimate", it is always useful to look at a plot of
>>> the sample variogram (as in the spatial case):
>>>
>>> plot(sample.stvgm, wireframe=T, scales=list(arrows=F))
>>>
>>> From the plot, one can read spatial and temporal ranges, the
>>> magnitude of variability and the ratios of nuggets and sills.
>>> Sometimes, one is better of by starting with a simpler model (e.g.
>>> "pure" metric), as the routines are more robust for less
>>> parameters, to get a first impression of the range of parameters.
>>> Furthermore, the plot suggests that a much smaller spatial upper
>>> boundary is sufficient, as the variogram surface is mainly flat in
>>> the spatial dimension.
>>>
>>> Spatio-temporal geostatistics remains tricky in terms of numerical
>>> routines and model interpretation, we try to give some guidance
>>> regarding gstat's functionalities in [1].
>>>
>>> HTH,
>>>
>>> Ben
>>>
>>> [1]
>>> https://journal.r-project.org/archive/accepted/na-pebesma-heuvelink.pdf
>>>
>>>
> Spatio-Temporal Interpolation using
> gstat<https://journal.r-project.org/archive/accepted/na-pebesma-heuvelink.pdf>
>>> journal.r-project.org CONTRIBUTED RESEARCH ARTICLE 1
>>> Spatio-Temporal Interpolation using gstat by Benedikt Gr�ler, Edzer
>>> Pebesma and Gerard Heuvelink Abstract We present new spatio ...
>>>
>>>
>>>
>>>
>>>
>>> On 21.06.2016 17:16, Dan Turenne wrote:
>>>> Hello again,
>>>>
>>>>
>>>> After some tinkering with my code last night I got the sample
>>>> variogram to work however now I have encountered another issue.
>>>> When attempting to use the fit.stVariogram function to fit the
>>>> sample variogram I receive this error:
>>>>
>>>>
>>>> Error in optim(extractPar(model), fitFun, ..., method = method,
>>>> lower = lower,  : non-finite value supplied by optim
>>>>
>>>>
>>>> even if I specify the upper and lower arguments to be finite.
>>>> Could this error be cause by the initial values that I set for
>>>> the space, time and joint variograms?  If so can anyone recommend
>>>> a good paper/article about how to fit these three variograms?  I
>>>> have attached a link to my script and data below, if anyone has
>>>> any insight as to why I am getting this error it would be
>>>> appreciated.  I tried to be as thorough as I could in my code
>>>> commenting, please let me know if anything is unclear.
>>>>
>>>>
>>>> https://www.dropbox.com/sh/7ccy5wu68gxf6sf/AACmk0AFxQvetunq4YPc9DXla?dl=0
>>>>
>>>>
>>>>
>>>>
> ________________________________
>>>> From: R-sig-Geo <r-sig-geo-bounces at r-project.org> on behalf of
>>>> Edzer Pebesma <edzer.pebesma at uni-muenster.de> Sent: June 21, 2016
>>>> 1:52 AM To: r-sig-geo at r-project.org Subject: Re: [R-sig-Geo]
>>>> Spatio Temporal kriging in Gstat
>>>>
>>>> This might be due to a bug (or feature) in the software caused by
>>>> the sparseness in your data, which might be different from that
>>>> used to test the software. Please make the data available
>>>> (off-list), along with an R script, so we can try to reproduce
>>>> the error message and look into it.
>>>>
>>>> On 21/06/16 04:11, Dan Turenne wrote:
>>>>> My apologies, I accidentally sent an unfinished email, here is
>>>>> the complete version of my question
>>>>>
>>>>>
>>>>> Hello R-Sig-Geo,
>>>>>
>>>>>
>>>>> As part of my masters thesis I am attempting to use
>>>>> spatio-temporal regression kriging to make predictions with
>>>>> temperature data, and I was hoping that someone might be able
>>>>> to give some insight as to how the algorithms work in gstat.
>>>>> My data consists of daily temperature observations from April 1
>>>>> to July 31, 2000.  There are observations from 164 stations
>>>>> across these 122 days, however not all stations have
>>>>> observations on all days, making for a total of 19282 records.
>>>>>
>>>>>
>>>>> I have tried to use an STSDF object but I have not had any
>>>>> success.  I created an sp object of length 164 with the station
>>>>> locations:
>>>>>
>>>>>
>>>>> sp = data.frame(long = stations$long, lat = stations$lat)
>>>>>
>>>>> coordinates(sp) = ~ long+lat
>>>>>
>>>>>
>>>>> Then I created a vector of length 122 with the times the
>>>>> observations were recorded and a data vector of length 19282:
>>>>>
>>>>>
>>>>> beginDate = as.Date(2000/04/01)
>>>>>
>>>>> endDate = as.Date(2000/07/31)
>>>>>
>>>>> times = as.POSIXct(seq(beginDate,endDate,by="days"))
>>>>>
>>>>>
>>>>> data=data.frame(temps$residual)
>>>>>
>>>>>
>>>>> And I also made an index detailing where observations are
>>>>> available, it looks like this with the first column
>>>>> representing spatial index and the second representing the time
>>>>> index
>>>>>
>>>>>
>>>>> 1   1
>>>>>
>>>>> 2   1
>>>>>
>>>>> 3   1
>>>>>
>>>>> 4   1
>>>>>
>>>>>
>>>>> st=STSDF(sp,time,data,index,endTime=delta(time))
>>>>>
>>>>>
>>>>> however when I try to calculate the sample variogram I get the
>>>>> following error:
>>>>>
>>>>>
>>>>> sample.stVariogram=variogramST(residual~1,data=st,
>>>>> tunit="days", tlags=1:7, progress=TRUE)
>>>>>
>>>>>
>>>>>
>>>>> Error in apply(do.call(cbind, lapply(ret, function(x) x$np)),
>>>>> 1, sum,  : dim(X) must have a positive length In addition:
>>>>> There were 50 or more warnings (use warnings() to see the first
>>>>> 50)
>>>>>
>>>>>
>>>>> All 50 of the errors are :
>>>>>
>>>>>
>>>>> In is.na(data[[as.character(as.list(formula)[[2]])]]) : is.na()
>>>>> applied to non-(list or vector) of type 'NULL'
>>>>>
>>>>>
>>>>> Can anyone see what I am doing wrong or give me any pointers?
>>>>> This error  is rather cryptic and I'm not quite sure what I'm
>>>>> doing wrong.  Any help would be appreciated.
>>>>>
>>>>>
>>>>> Many Thanks,
>>>>>
>>>>> Daniel Turenne
>>>>>
>>>>> University of Manitoba
>>>>>
>>>>>
>>>>> [[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
>>>> R-sig-Geo - R Special Interest Group on using Geographical
>>>> ...<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>>> stat.ethz.ch R-sig-Geo -- R Special Interest Group on using
>>>> Geographical data and Mapping About R-sig-Geo
>>>>
>>>>
>>>>
>>>>>
>>>>
>>>> -- Edzer Pebesma Institute for Geoinformatics  (ifgi),
>>>> University of M?nster Heisenbergstra?e 2, 48149 M?nster, Germany;
>>>> +49 251 83 33081 Journal of Statistical Software:
>>>> http://www.jstatsoft.org/ Computers & Geosciences:
>>>> http://elsevier.com/locate/cageo/ Spatial Statistics Society
>>>> http://www.spatialstatistics.info
>>>>
>>>>
>>>> [[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
>>>>
>>>
>>> -- Dr. Benedikt Gr�ler Institute of Hydrology Ruhr University
>>> Bochum Universit�tsstra�e 150 44801 Bochum
>>>
>>> ben.graeler.org
>>>
>>> +49 (0234) 32 - 27619
>>>
>>> _______________________________________________ R-sig-Geo mailing
>>> list R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>> [[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
>>>
>>
>> -- Dr. Benedikt Gräler Institute of Hydrology Ruhr University Bochum
>> Universitätsstraße 150 44801 Bochum
>>
>> ben.graeler.org
>>
>> +49 (0234) 32 - 27619
>>
>> _______________________________________________ R-sig-Geo mailing
>> list R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>> [[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
>>
>
> --
> Dr. Benedikt Gräler
> Institute of Hydrology
> Ruhr University Bochum
> Universitätsstraße 150
> 44801 Bochum
>
> ben.graeler.org
>
> +49 (0234) 32 - 27619
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>        [[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
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/
Spatial Statistics Society http://www.spatialstatistics.info


	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list