Hi all,

I am new to geostatistics and need some help on projection. I am working on a bivariate model of a sphere and want to do cokriging with gstat package in R. My data look like

> head(mydata)
         lon         lat         v1         v2
1  2.1238698  0.87854240  1.8051750  1.2222044
2 -2.6524104 -0.14883920 -0.4271378 -0.1659977
3 -2.5863044 -0.18126753 -0.4142029 -0.2782104
4  1.1100337 -0.03636767 -1.1252408  0.4394822
5 -0.5421093 -0.60182344  0.8544659  0.6174832
6 -3.0686334  0.53252533  0.9143837  0.9566044

where lon and lat are longitude and latitude in radian. v1and v2 are two collocated variables. The total number of points is 100 (only 6 are shown here). I want to do cokriging at a couple of new locations saved as a data frame called newloc.

> newloc
  lon lat
1 0.2 0.3
2 0.5 1.0

My thought is to first get the sample variogram/cross variogram and use linear model of coregionalization to fit them and then do cokriging. Here is my code (basically following demo(cokriging)in gstat). I first simply neglected the sphere and just treated lon and lat as Cartesian coordinate.

g <- gstat(id = "val1", formula = v1~1, locations = ~lon+lat, data = mydata)
g <- gstat(g, id = "val2", formula = v2~1, locations = ~lon+lat, data = mydata)
g <- gstat(g, model = vgm(1, "Sph", 1, 1), fill.all = TRUE)
x <- variogram(g, cutoff = 2)
myfit <- fit.lmc(x, g)
plot(x, model = myfit)
newloc <- data.frame(lon = c(0.2, 0.5), lat = c(0.3, 1.0))
z <- predict(myfit, newdata = newloc)

I don't know if the above code is correct or not but it gives me some output. Now comes the question of how to taking into account the curvature of the sphere. I combed through some R-sig-Geo posts and found something using proj4string or proj4 but I can't make it work. I really need some help on how to do projection in computing variogram/cross variogram and doing cokriging.

Please help me on this. Thank you very much.

Zhaofeng
	[[alternative HTML version deleted]]

