Skip to content

GStat kriging using spherical geometry

2 messages · hgreatrex, Edzer Pebesma

Original
hgreatrex
# 
Dear all,

Does anyone know if the Gstat package in r offers kriging using 
spherical geometry rather than euclidean distances?

I couldn't find anything in ?krige or online,  but I've heard that the 
general gstat package includes this option.  From comparisons with a 
different kriging code outside r, the inclusion of spherical geometry 
does make a difference to my results, but I'd really like to use gstat!

For some context, I'm block kriging rainfall data using lat long 
coordinates.  An example of my code is below.  To make things simple, in 
this example I'm aiming to krige over a 1 degree block centred at (39.5, 
7.5)
(So the corners of the block would be [39,7],[40,7],[40,8],[39,8])
Note, In general I wouldn't be kriging over a block this size - it's 
just to make the code simple!

Best wishes and thank you for your time
Helen

---------------------------------------------------------------------------

 ># Read in the input files
 >##############################
 >  dayinput <- file("zsourcep.txt","r")
 >  mydata <- read.table(dayinput,header=TRUE)
 >  unlink("zsourcep.txt")
 >
 >
 ># Read in the output files
 >##############################
 >  dayoutput <- file("ztarget.txt","r")
 >  mydataout <- read.table(dayoutput,header=TRUE)
 >  unlink("ztarget.txt")
 >
 >
 ># Turn these into the correct format for krige
 >##############################
 >  cinput <- data.frame (Long = mydataout$Long, Lat= mydataout$Lat     
                              )
 >  input  <- data.frame  (Long = mydata$Long ,     Lat=mydata$Lat,      
Value=mydata$Value)
 >  coordinates(input)  <- ~ Long + Lat
 >  coordinates(cinput)<- ~ Long + Lat
 >
 >
 ># Now make the variogram model
 >##############################
 >  m <- vgm(.675, "Exp", 9.51, 0.237)
 > >
 ># And krige
 >##############################
 >  x <- krige(Value ~ 1, input,cinput , model = m, block=c(1,1))
  [using ordinary kriging]
 >  x
     coordinates var1.pred  var1.var
Edzer Pebesma
# 
I find it hard to believe that you did a thorough online search. 
Googling on "kriging spherical distance" gave me the following as first hit:

https://stat.ethz.ch/pipermail/r-sig-geo/2008-October/004457.html

and then there the search facility for the archives of this list, 
referenced at the end of every message.

In your example, don't forget to set the proj4string of youe objects. 
Further, the exponential model does, afaik, not give positive definite 
covariance matrices on the sphere. I have received no reactions to the 
question in the email mentioned above.

Attempts for block kriging should result in an error message, as I can't 
see what is meant by a constant block size on the sphere--something like 
one degree by one degree?
--
Edzer
Helen Greatrex wrote: