r-spatial / gstat

Spatial and spatio-temporal geostatistical modelling, prediction and simulation
http://r-spatial.github.io/gstat/
GNU General Public License v2.0
196 stars 50 forks source link

Question about "Sill" setting in vgmST #48

Closed bhurley6ehs closed 4 years ago

bhurley6ehs commented 5 years ago

Hi all:

Writing to ask what exactly the sill setting in "separable" vgmST is for. I admit I played around with it extensively but could not find a lot of rhyme or reason to it. Thank you!

bhurley6ehs commented 5 years ago

Hi all:

Another question: When fit.StVariogram is called, and then checked, You get a Space Component with a range, and a Time Component with a range. My Space component clearly is reporting in meters, which is great (I think), but is the time component reported back in a default of days? Or is it in t-units as given to STFDF? My T-units were months when the data were first loaded. I am just trying to figure out, if I have 340 months (340 temporal units given to gstat), what a final reported temporal range of 2813 would mean (I think it's days, which is fine. If it's months, that's an issue since I don't even have 2813 months worth of data). See attached for the readout I am looking at: Readout_Sent_to_Gstat

bhurley6ehs commented 5 years ago

Hi all:

Just checking in on this.

edzer commented 5 years ago

For the separable model, both components are essentially correlograms (psill values sum to 1); the sill: reported value is the value their product is multiplied with to get a covariance.

bhurley6ehs commented 5 years ago

Edzer:

Thank you. Under the time component, it reports a value in the range of 2813.319. What unit is this in? As I said I fed the STFDF a file broken into months (say 100 entries in one group for January, 100 for February, etc. ), for 340 something months. I get the wireframe back in days (which is perfectly fine and readable). I am just curious what the 2813.319 unit it. I apologize if this is obvious and something I have missed. My guess is days since I don't HAVE 2813.319 months of data.

bhurley6ehs commented 5 years ago

Hi all:

Just checking in on this.

edzer commented 5 years ago

What does the 2813 exactly refer to?

bhurley6ehs commented 5 years ago

Hi Edzer:

Please see attached. Refer to Spherical Range under

the "Time" component.

Sep_Exp_Sph

edzer commented 5 years ago

How is your time variable represented: POSIXct or Date?

bhurley6ehs commented 5 years ago

Hi Edzer:

(1) Thank you for getting back to me.

(2) I am currently converting my column to POSIXct.

(3) I have attached a sample of the excel table I am using for your reference. The column I am using is DTE_USED. It consists of a once-a-month sample, in the same spatial location, over a ~ 28 year period (341 months), with no gaps. As the table is laid out, and from your JSS (2012, page 2) paper, the table presents as a long-format ST table, of which I am currently feeding into STFDF with no problems.

ssamatab_platformed_for_R__to_Edzer.xlsx .

edzer commented 5 years ago

In case your time is POSIXct, I'd say the range parameter estimate is in (unit) seconds.

bhurley6ehs commented 5 years ago

Edzer:

Thank you much. I have checked and re-imported my data as Date.

I realize this may be too general a question, but when I am working to build a ProductSum, Metric, SumMetic, and SimpleSumMetric, model I am consistently getting range estimations that just reflect the settings I am inputting into fit.StVariogram.

As an example, the below code shows a given by me (under spaceVGMx2) spatial range of 300000m, and a temporal range of 2500 days (under timeVGMx2). When the fitted variogram is called, it consistently (and only) just reports back the ranges as given (the fitted model just shows the input ranges (300000m and 2500 days).

Is there a "dead give-away" why this would be consistently happening?

Thank you again.

spaceVGMx2 <- vgm(psill = .0004, model = "Exp", range = 300000, nugget = .0000175) ###Space timeVGMx2 <- vgm(psill = .0004, model = "Exp", range = 2500, nugget = .0000175)###Time ProdSum2 <- vgmST("productSum", space=spaceVGMx2, time =timeVGMx2, k=1) plot(var3, ProdSum2, wireframe=T, xlab = "distance (m)", zlab = "gamma", main = "Product Sum: Exponential (s) and Exponential (t)")

estiStAni(var3, c(1, 300000), method = "metric", spaceVGMx2, timeVGMx2, s.range = 2000000, t.range = 2500)

estiStAni = 150000

lower and upper bounds

pars.la <- c(sill.s = 0.000001, range.s = 10, nugget.s = 0, sill.t = 0.000001, range.t = 1, nugget.t = 0, sill.st = 0.000001, range.st = 2000, nugget.st = 0, anis = 0.000001) pars.ua <- c(sill.s = 600, range.s = 300000, nugget.s = 1000, sill.t = 7000, range.t = 5000, nugget.t = 1000, sill.st = 200000, range.st = 500000, nugget.st = 1000, anis = 200000)

ProdSum_Vgm2 <- fit.StVariogram(var3, ProdSum2, fit.method = 3, stAni = 150000, method = "L-BFGS-B") ProdSum_Vgm2

edzer commented 5 years ago

Well, it's a non-linear optimization: if the initial values are not good, or the problem is hard to optimize, then there's little hope you will get to a real optimum.

bhurley6ehs commented 5 years ago

Well the separable fit produced reasonable results after some tweeking. In your opinion would those results be somewhat trustworthy?

Thank you again!