[STATSGRASS] gstat probability maps

[Pierluigi De Rosa]

Dear Edzer & Wolf

I read the all thread and some papers... At the end I think Edzer is right, so 
you can use the kriging variance to evaluate a propability map of exceeding a 
value!
Thanks to all for your help!
Pierluigi

Alle 19:34, martedì 27 marzo 2007, Edzer J. Pebesma ha scritto:
> Dear Wolf & Pierluigi,
>
> I believe that it is another common misunderstanding to believe that
> indicator kriging is needed to obtain (estimates of) probabilities. You
> can of course do it, try something along the lines of
>
> v <- variogram(I(zinc>500)~1, meuse)
> m <- vgm(.15, "Sph", 874, .1) # just filling in some values here
> ik <- krige(I(zinc>500)~1, meuse, meuse.grid, m)
>
> but then you give up any information on how far an observation is below
> or above 500, as data are transformed to 0 (below 500) and 1 (above 500)
> values. IMO this is valuable information, and although you may criticize
> the Gaussian assumption (on the log-scale) underlying my previous
> suggestion, at least it leads to "true" probabilities, instead of
> interpolated 0 and 1 values that may (and will) as well lie outside the
> interval [0,1].
>
> I wouldn't go as far as advising never to use indicator kriging, but the
> automatic "for probabilities you need indicator kriging" clearly misses
> important essentials of geostatistics.
>
> Both issues, the need to simulate or the need for indicator transforms
> to obtain probabilities, have been pushed by certain authors that know
> better, and re-arise from time to time. A third and related issue is the
> suggestion that the (ordinary) kriging variance is not a measure for
> prediction error; see the following ai-geostats thread:
> http://www.mail-archive.com/ai-geostats%40jrc.it/msg02890.html
>
> Best wishes,
> --
> Edzer
>
> Wolf Bergenheim wrote:
> > Pierluigi,
> >
> > That would be indicator kriging:
> >
> > I use gstat, so I can't give you any specific R commands, but the
> > principle is the same: You have to create a new dataset which is 1 if
> > the threshold is exceeded, and 0 if not. Now your data is 0 or 1. Then
> > do normal ordinary kriging on that new data column. The resulting map
> > shows the probability of exceeding the threshold. In gstat you have to
> > add  extra checks to see that your result is strictly [0.0 - 1.0]
> >
> > --Wolf
> >
> > On 27.03.2007 18:18, Pierluigi De Rosa wrote:
> >> Dear All,
> >>
> >> I need to prepare a probability map of exceeding a value...
> >> I can make the simulation, for example I tried to make it with data
> >> meuse:
> >>
> >> v <- variogram(log(zinc)~1,loc= meuse)
> >> m <- vgm(.59, "Sph", 874, .04)
> >> sim <- krige(formula = log(zinc)~1, locations = meuse, model = m,
> >> newdata = meuse.grid, nmax = 15, nsim = 5)
> >> names(sim)
> >> [1] "sim1" "sim2" "sim3" "sim4" "sim5"
> >>
> >> Now How can I create a probabilility map of exceding a value (for
> >> example 5)??
> >>
> >> Thanks
> >> Pierluigi

-- 

Pierluigi De Rosa

Department of Enviromental & Civil engineering 
Faculty of Engineering 
Perugia (Italy)