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

Edzer Pebesma edzer.pebesma at uni-muenster.de
Mon Feb 22 20:06:00 CET 2010



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
> 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
>>>       
>
>
>
>   

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster 
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251 
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de 
http://www.52north.org/geostatistics      e.pebesma at wwu.de



More information about the R-sig-Geo mailing list