[R-sig-Geo] Cell declustering examples in Gstat

[Dylan Beaudette]

On Monday 22 February 2010, Edzer Pebesma wrote:
> Dylan Beaudette wrote:
> > On Monday 22 February 2010, Edzer Pebesma wrote:
> >> Carlos, please try:
> >>
> >> library(gstat)
> >> data(meuse)
> >> meuse$rn = 1:nrow(meuse)
> >> coordinates(meuse)=~x+y
> >> spplot(meuse["rn"])
> >> data(meuse.grid)
> >> gridded(meuse.grid)=~x+y
> >> out = krige(rn~1, meuse, meuse.grid, nmax=1)
> >> spplot(out[1])
> >> meuse$w = as.vector(table(out[[1]])/sum(table(out[[1]])))
> >> spplot(out[1],col.regions=bpy.colors(),sp.layout=list("sp.points",meuse)
> >>) spplot(meuse["w"], col.regions=bpy.colors())
> >>
> >> Having a somewhat finer interpolation grid might give slightly more
> >> satisfactory results.
> >>
> >> Best regards,
> >> --
> >> Edzer
> >
> > Thanks for the example Edzer. This is a really interesting demonstration
> > of cell declustering... however, pardon my lack of understanding, I can't
> > quite figure out how the resulting are obtained by interpolating the row
> > numbers. Can you please elaborate on how/why the use of row numbers and
> > setting nmax=1 produces these results?
>
> Hi Dylan,
>
> When "interpolating" with nmax=1, you basically assign the value of the
> nearest observation to each grid cell, so, honoustly, it's hard to call
> this interpolation, it is rather something of a discretized Thiessen
> polygon. At least, this is the case for inverse distance (used in the
> example above) and ordinary/universal kriging.
>
> Challenged to avoid the krige function, the following script assigns row
> indexes to all nearest grid cells (when executed after the one above),
> using spDists from sp to find the right index:
>
> d = spDists(meuse, meuse.grid)
> dim(d)
> out$x = apply(d, 2, function(x) { which(x == min(x)) })
> sum(abs(out[[1]]-out$x))
> # should give 0
>
> following that you could use the table stuff. Of course the spDists
> function computes all n x N distances, and is not recommended for large
> / massive data sets. The idw / krige functions in gstat use a PR bucket
> quadtree (not written by me!) to efficiently find the nearest nmax
> points; this does scale up.
>
> Best wishes,
> --
> Edzer

Excellent. Thank you for the extended discussion Edzer. I'll try and work up a 
small example comparing the weights derived from this method with those 
derived from Thiessen polygons.

Cheers,
Dylan




> > Thanks!
> >
> > Dylan
> >
> >> Carlos Alberto Gutierrez Carvajal wrote:
> >>> Hello everybody in the list.
> >>>
> >>> My question have been treated previously by Edzer Pebesma in his
> >>> communications with Stefano Pegoretti in july 2007, however I would
> >>> find very useful some examples.
> >>>
> >>> Basically, Edzer Pebesma pointed to two methods to perform cell
> >>> declustering.
> >>>
> >>> One of these methods is voronoi diagrams (avaliable in package deldir).
> >>>
> >>> The other method is, in words of Edzer Pebesma, by  'using the number
> >>> of nearest cells based on a regular discretization of the study area;
> >>> for the latter you could misuse package gstat, interpolate record
> >>> number with nmax=1 and compute a table of the resulting "prediction"
> >>> grid'.
> >>>
> >>> Could you provide me some examples for each one of this methods for an
> >>>  R user,please?.
> >>>
> >>> Regards,
> >>>
> >>> Carlos Gutierrez Carvajal.
> >>>
> >>> _______________________________________________
> >>> R-sig-Geo mailing list
> >>> R-sig-Geo at stat.math.ethz.ch
> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341