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
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 useThank 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-boun...@r-project.org> on behalf of Benedikt Graeler <benedikt.grae...@rub.de> Sent: June 28, 2016 5:23 AM To: r-sig-geo@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
Spatio-Temporal Interpolation using gstat<https://journal.r-project.org/archive/accepted/na-pebesma-heuvelink.pdf>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-boun...@r-project.org> on behalf of Benedikt Graeler <benedikt.grae...@rub.de> Sent: June 23, 2016 2:54 AM To: r-sig-geo@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
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-boun...@r-project.org> on behalf of Edzer Pebesma <edzer.pebe...@uni-muenster.de> Sent: June 21, 2016 1:52 AM To: r-sig-geo@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@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geoR-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@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@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@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@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@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@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
