Skip to content

Cell declustering examples in Gstat

5 messages · Carlos Alberto Gutierrez Carvajal, Dylan Beaudette, Edzer Pebesma

Original
Carlos Alberto Gutierrez Carvajal
# 
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.
Edzer Pebesma
# 
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
Carlos Alberto Gutierrez Carvajal wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Dylan Beaudette
    

    

      
      #
    

  


  

      
On Monday 22 February 2010, Edzer Pebesma wrote:

            

      
      
      
    

        
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?

Thanks!

Dylan

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
Dylan Beaudette wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Dylan Beaudette
    

    

      
      #
    

  


  

      
On Monday 22 February 2010, Edzer Pebesma wrote:

            

      
      
      
    

        
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