bk.grid.R 5.03 KB
Newer Older
1

agebhard's avatar
agebhard committed
2

3
bk.grid <- function(point.obj,
agebhard's avatar
agebhard committed
4
5
6
7
8
                    at,
                    prior,
                    var.mod.obj,
                    xsw=NULL,ysw=NULL,xne=NULL,yne=NULL,
                    dx=NULL,dy=NULL,nx=NULL,ny=NULL,
alge's avatar
alge committed
9
                    angle=0,
agebhard's avatar
agebhard committed
10
11
12
13
14
15
16
17
18
19
20
                    maxdist = NULL,
                    extrap = F,
                    border=NULL,
                    trend=0,
                    rsearch=0,
                    nsearch=0,
                    nsmin=-1,
                    nsmax=-1,
                    mode=3,
                    duplicate = "error",
                    dupfun = NULL,
21
22
                    method="gqr",
 		    get.lm=FALSE)
23
{
24
25
    tmp   <- check.gridparams(angle,xsw,xne,ysw,yne,
                              dx,dy,nx,ny)
26
27
28
29
30
31
32
33
34
    angle <- tmp$angle
    xsw   <- tmp$xsw
    xne   <- tmp$xne
    ysw   <- tmp$ysw
    ysw   <- tmp$ysw
    dx    <- tmp$dx 
    dy    <- tmp$dy 
    nx    <- tmp$nx 
    ny    <- tmp$ny 
35
    nz    <- tmp$nz 
agebhard's avatar
agebhard committed
36

alge's avatar
alge committed
37
    dog   <- check.border.grid(extrap,xsw,xne,ysw,yne,nx,ny,border,point.obj)
agebhard's avatar
agebhard committed
38

39
40
41
42
43
    
    point.obj <- remove.duplicates(point.obj,at,duplicate,dupfun)
    tmp       <- check.krigedata(point.obj,at,var.mod.obj,mode)
    n         <- tmp$n

44
45
46
47
48
    tmp     <- check.searchparams(maxdist,rsearch,nsearch,nsmin,nsmax,n)
    rsearch <- tmp$rsearch
    nsmin   <- tmp$nsmin
    nsmax   <- tmp$nsmax
    
agebhard's avatar
agebhard committed
49
50
51
52
53
54
55
56
57
58
59
60
61
62
    if(trend==0) ntrend<-1
    if(trend==1) ntrend<-3
    if(trend==2) ntrend<-6

    covtype<-switch(attr(var.mod.obj,"type"),
                    exponential=1,
                    gaussian=2,
                    spherical=3,
                    linear=4,
                    0)
    cov0<-0
    
#   determine optimum array sizes:
    if(!is.null(method)){
63
64
        if(method!="gqr" && method!="direct" && method!="ols")
            stop("method (used for glsfit) should be one of \"gqr\", \"ols\" or \"direct\"!")
agebhard's avatar
agebhard committed
65
    } else {
66
        method <-"gqr"
agebhard's avatar
agebhard committed
67
    }
68
    method<-switch(method,direct=2,gqr=1,ols=0)    
alge's avatar
alge committed
69
70
#    lwork <- glsfit.workquery(n,ntrend,method)
    lwork <- 3000
alge's avatar
alge committed
71

72
73
    if(prior$ntr!=ntrend)
        stop("model order of priors does not match!")
agebhard's avatar
agebhard committed
74
75
76
    npr<-prior$n
    typpr<-prior$info
    typpr[prior$type=="subjective"]<-typpr[prior$type=="subjective"]*(-1)
agebhard's avatar
??    
agebhard committed
77

alge's avatar
alge committed
78
79
80
81

        
    snbbit<-rep(0,1+n*nz)
    snbbit[1]<-1
82
83
84
85
86
87
88
89
90
    if(get.lm){
      mu <- double(ntrend*nz)
      lambda <- double(n*nz)
      lambd0 <- double(nz)
    } else {
      mu <- double(ntrend)
      lambda <- double(n)
      lambd0 <- double(1)
    }
alge's avatar
alge committed
91
92
    

agebhard's avatar
agebhard committed
93
    ans<-.C("bk_grid",
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
            xsw     = as.double(xsw),
            ysw     = as.double(ysw),
	    xne     = as.double(xne),
            yne     = as.double(yne),
            angle   = as.double(angle),
	    nx      = as.integer(nx),
	    ny      = as.integer(ny),
            dx      = as.double(dx),
            dy      = as.double(dy),
	    xg      = double(nx),
	    yg      = double(ny),
	    zg      = double(nz),
	    varg    = double(nz),
	    dog     = as.integer(dog),
            lon     = as.double(point.obj$x),
            lat     = as.double(point.obj$y),
            z       = as.double(point.obj[[match(at, names(point.obj))]]),
	    extrap  = as.integer(extrap),
	    n       = as.integer(n),
	    covtype = as.integer(covtype),
            covpar  = as.double(var.mod.obj$parameters),
agebhard's avatar
agebhard committed
115
116
117
	    covmat  = double(n*n), 	
	    ldcov   = as.integer(n),
            extcov  = as.integer(0), # no external cov matrix
118
119
120
	    trend   = as.integer(trend),
	    ntrend  = as.integer(ntrend),
            mupr    = as.double(matlist.cbind(prior$mu)),
agebhard's avatar
agebhard committed
121
            ldmpr   = as.integer(ntrend),
122
            phipr   = as.double(matlist.cbind(prior$phi)),
agebhard's avatar
agebhard committed
123
124
125
            ldphpr  = as.integer(ntrend),
            lonpr   = as.double(prior$lon),
            latpr   = as.double(prior$lat),
126
127
128
129
130
131
132
133
            npr     = as.integer(npr),
            typpr   = as.integer(typpr),
            rsearch = as.double(rsearch),
	    nsearch = as.integer(nsearch),
	    nsmin   = as.integer(nsmin),
	    nsmax   = as.integer(nsmax),
            lwork   = as.integer(lwork),
            mode    = as.integer(mode),
134
135
136
            mu      = as.double(mu),
            lambda  = as.double(lambda),
            lambd0  = double(lambd0),
alge's avatar
alge committed
137
            bits    = as.integer(c(integer(nz),snbbit)),
138
            ierr    = integer(1),
139
            get.lm  = as.integer(get.lm),
140
            glsmth  = as.integer(method),
141
142
#           ,.Package= "baykrig"
            )
agebhard's avatar
agebhard committed
143
#browser()
agebhard's avatar
agebhard committed
144
145
    retval<-list(x=ans$xg,
                 y=ans$yg,
agebhard's avatar
agebhard committed
146
147
148
                 z=matrix(ans$zg,nx,ny),
                 var=matrix(ans$varg,nx,ny),
                 done=matrix(ans$dog, nx, ny),
149
                 snb=matrix(ans$bits[(nz+2):(nz+n*nz+1)],nrow=n,ncol=nz,byrow=F)
alge's avatar
alge committed
150
                 )
151
152
153
154
155
    if(get.lm){
      retval$lambda <- matrix(ans$lambda,nrow=n,ncol=nz,byrow=FALSE)
      retval$lambda0 <- matrix(ans$lambd0,nx,ny)
      retval$mu <- matrix(ans$mu,nrow=n,ncol=ntrend,byrow=FALSE)
    }
agebhard's avatar
agebhard committed
156
157
    retval$z[retval$done<=0] <- NA
    retval$var[retval$done<=0] <- NA
alge's avatar
alge committed
158
    retval$lambda0[retval$done<=0] <- NA
agebhard's avatar
agebhard committed
159
    retval$data<-point.obj
160
    retval$at<-at
agebhard's avatar
agebhard committed
161
162
163
    class(retval)<-"krige.map"
    retval    
  }