Skip to content

block and output grid size in block kriging

8 messages · subash, Edzer Pebesma

Original
subash
# 
Dear Group,

I was using bock kriging  for groundwater level data. I was wondering, if
the output.grd size and the block size can be different.

For example, In the example for the case given below, I need block kriged
average of 5000x5000. The output grid size given is 1000m*1000m. Is it
correct?

Or should the output.grd size also be 5000mx5000m as in block size?


# IMP : There are the centre coordinates for block kriging and not corner
coordiantes.
Xmin = 581000 + 2500
/Xmax = 711000 - 2500
Ymin = 1266000 + 2500
Ymax = 1371000 - 2500

output.grd = expand.grid(x=seq(from=Xmin, to=Xmax,by=1000),y=seq(from=Ymin,
to=Ymax, by=1000))
coordinates(output.grd) = c("x", "y")
gridded(output.grd) = TRUE

krig =
krige(borewell.piezo[,3]~1,location~Easting+Northing,data=borewell.piezo,newdata=output.grd,model
= vario.fit,block=c(5000,5000))/


Thanks & Regards,
Subash



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/block-and-output-grid-size-in-block-kriging-tp7587770.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
Edzer Pebesma
#
On 02/09/2015 08:48 PM, subash wrote:
Yes, they can be different. Blocks are centered over grid cell centers,
and may overlap, be identical to grid cells, or smaller than grid cells
and have space between them.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      subash
    

    

      
      #
    

  


  

      
Dear Edzer,

Thanks for quick reply. A question still lingering in my mind is below

In a similar lines I was using block kriging of rainfall with block size of
25x25Km. The output grid  chosen is 1x1Km. This is resulting in kriging
estimate computed at 1x1Km . However, my interest is to get 25x25Km kriging
estimate.

The approach I tried was setting the output grid mesh as 25Km, but this is
wrong as the block average would be arrived from 1 point only.

output.grd = expand.grid(x=seq(from=tmp$min[1],
to=tmp$max[1],by=1),y=seq(from=tmp$min[2], to=tmp$max[2], by=1))
coordinates(output.grd) = c("x", "y")
gridded(output.grd) = TRUE

block_Krig = krige(stn.rain.data[,3]~1,location = ~easting+northing, data =
stn.rain.data,newdata=output.grd,model = clim.vrmod,block=c(25,25)) 

What is that I am missing.


Thanks in advance
Subash



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/block-and-output-grid-size-in-block-kriging-tp7587770p7587773.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

  
  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
On 02/10/2015 06:21 AM, subash wrote:

            

      
      
      
    

        
So, set up a prediction grid (newdata) with 1 km x 1 km, and specify
block = c(25, 25), assuming km is your distance unit - that seems to be
what you do below.

            

      
      
      
    

        
I don't understand this sentence.

            

      
      
      
    

        
I don't know. For more readable code, I'd suggest:

# assume col 3 in stn.rain.data is called "prec"
coordinates(stn.rain.data) = ~easting + northing
block_Krig = krige(prec~1, stn.rain.data, output.grd, clim.vrmod,
	block = c(25,25))

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    
2 days later

    

      
      

        


  

        

      subash
    

    

      
      #
    

  


  

      
Dear Edzer, 

To understand why block kriging results vary,with different grid sizes, to
keep gridsize identical to grid cells, and smaller than grid cells.

library(gstat)
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y

lzn.vgm<- variogram(log(zinc)~1, meuse)
lzn.fit1 <- fit.variogram(lzn.vgm, model=vgm(psill=1, model="Sph",
range=900, nugget=1))
plot(lzn.vgm, lzn.fit1)

# Ordinary kriging
lzn.ok <- krige(log(zinc)~1, meuse, meuse.grid, model = lzn.fit1)

#Block kriging 
lzn.bok <- krige(log(zinc)~1, meuse, meuse.grid, model = lzn.fit1,
block=c(40,40))

# Choosing a random location with the centre coordiantes(xx,yy) and find the
corresponding block kriged estimate
xx = 181140
yy = 333100
est.bk = lzn.bok$var1.pred[(lzn.bok$x==xx)&(lzn.bok$y==yy)]
[1] 6.499158
est.ok =lzn.ok$var1.pred[(lzn.ok$x==xx)&(lzn.ok$y==yy)]
[1] 6.499624

Why are the kriging estimates at a point different in block and ordinary
kriging, when the grid size is isentical to (40 x 40) or less than(20x20)
the meuse data.

Thanks







diff= m.bk - est.bk



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/block-and-output-grid-size-in-block-kriging-tp7587770p7587780.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

  
  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
On 02/12/2015 04:19 PM, subash wrote:

            

      
      
      
    

        
because lzn.ok contains ordinary point kriging estimates, for points
laid out on a grid, and lzn.bok contains ordinary 40 m x 40 m block
kriging estimates: a different quantity is predicted.

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      subash
    

    

      
      #
    

  


  

      
Dear Edzer,

Thanks for the clarification.Still it not clear to me.

I added an image below showing the block and the grids. If my visualization
is correct then for a block of (40,40), there is only 1 grid points. If Yes,
then ordinary kriging estimate and block kriging estimate for block size
(40,40) should be same.

<http://r-sig-geo.2731867.n2.nabble.com/file/n7587782/block_kriging_framework.png> 

Thanks in advance



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/block-and-output-grid-size-in-block-kriging-tp7587770p7587782.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

  
  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
The drawing (thanks!) suggests that you think that the point over which
a kriging block is centered is not the center point of a grid cell. This
is wrong: the coordinates of grid cells, obtained e.g. by

require(sp)
demo(meuse, ask = FALSE)
coordinates(meuse.grid)

are grid cell centers. You can verify this, e.g. for the first 20 grid
cells / pixels with

plot(as(meuse.grid[1:20,], "SpatialPolygons"))
points(meuse.grid[1:20,])

I also tried to point out in my previous answer that you are not looking
at "ordinary kriging" vs "block kriging", but at "ordinary point
kriging" vs. "ordinary block kriging". If your observations are not
constant and you don't use a pure nugget model, point kriging and block
kriging estimates will always be different.

        
On 02/12/2015 05:22 PM, subash wrote: