[R-sig-Geo] Spatio Temporal kriging in Gstat
Benedikt Graeler
benedikt.graeler at rub.de
Tue Jun 28 12:23:54 CEST 2016
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
>
>
> 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
More information about the R-sig-Geo
mailing list