[R-sig-Geo] Univeral kriging

[Edzer Pebesma]

Tankiso,

Here is some tips:
- make your code such that cut and paste works for us; for me, the ^2
below is vanished (and interpreted as superscript right away, I don't
know why)
- you didn't test your code, as the function inv does not exist
- your code doesn't compute the kriging variance
- the code is very hard to read, and uncommented, so we have to guess
what you want (first order linear trend in the coordinates?)
- use functions if you want to make something generic
- write and use functions such that matrices are allowed, e.g. for g()
you can pass a matrix vector; don't write loops over the individual elements
- If you want to solve Ax=b, write solve(A,b) instead of solve(A) %*% b.
- if you want to find out if your code works, don't post it here but
compare your output with that obtained from packages that do universal
kriging such as fields, geoR or gstat.
--
Edzer

TANKISO THEJANE wrote:
> R geo-helpers
>  
> i am new to r and would like you to help me as i am trying to go through basic steps. I need to compute a predicted value at point (3,3) and its variance given la set of locations and value. i want to make the code generic for any size such that i would be able to estimate universal krige at each point of grid (0,..,6) (0,..6).
>  
> the code
>  
> g = function(h) { return(0.1+3*(1-exp(-h/4)))}## variogram of exponential form
> x = rbind(c(1,1),c(2,5),c(4,1),c(5,4))##4 observed locations 
> x = rbind(x,c(3,3))  # Row five is the point at which we want to estimate.
> d = sqrt(x[,1]^2+x[,2]^2)##distance of the locations from the origin
> s = c(1,9,5,12)## values at the corresponding 4 locations
> ## pairsv = pair semivariogram
> pairsv = function(i,j) {
>   p1 = x[i,]
>   p2 = x[j,]
>   d = sqrt(sum((p1-p2)^2))
>   return(g(d))
> }
> G = matrix(0,4,4)
> for(i in 1:4) for(j in 1:4) G[i,j] = pairsv(i,j)
> X = matrix(0,4,3)
> for(i in 1:4) for(j in 1:3) X[i,j] = d[i]^(j-1)
> ##Asymmetric matrix 
> Gu = rbind(cbind(G,X), cbind(t(X),matrix(0,3,3)))
> x
> gu = (1:7)*0
> gu
> for(i in 1:4) gu[i] = pairsv(i,5)
> for(i in 1:3) gu[i+4] = d[5]^(i-1)
> lu = inv(Gu) %*% gu
>
>  
>  
>
>
>       
> 	[[alternative HTML version deleted]]
>
>   
> ------------------------------------------------------------------------
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>   

-- 
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