[R-sig-Geo] GStat kriging using spherical geometry

[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:
> 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
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.springer.com/978-0-387-78170-9 e.pebesma at wwu.de