Skip to content

GitLab

Commit b022db9c authored by agebhard's avatar agebhard
Browse files

Imported sources

parents
Hide whitespace changes
Inline Side-by-side
R/bk.R 0 → 100644
View file @ b022db9c
bk.tiles <- function(xsw,ysw,xne,yne,
dx,dy,itx,ity,
angle=NULL,
point.obj,
at,
var.mod.obj,
maxdist = NULL,
extrap = FALSE,
trend=0,
rsearch=0,
nsearch=0,
nsmin=-1,
nsmax=-1,
mode=1)
{
if(is.null(angle))angle<-0
dgx<-xne-xsw
dgy<-yne-ysw
nx<-ceiling(dgx/dx)
ny<-ceiling(dgy/dy)
ntx<-ceiling(nx/itx)
nty<-ceiling(ny/ity)
at <- point.obj[[match(at, names(point.obj))]]
n<-length(point.obj$x)
nz<-nx*ny
nt<-ntx*nty
ipt<-itx*ity
dog<-matrix(1,nx,ny)
if (!extrap) {
dog <- in.convex.hull(tri.mesh(point.obj$x, point.obj$y,
duplicate = "remove"),
point.obj$x, point.obj$y)*dog
dog<-matrix(dog,nx,ny)
}
extrap<-1*extrap
if (!inherits(point.obj, "point"))
stop("point.obj must be of class, \"point\".\n")
if (!inherits(var.mod.obj, "variogram.model"))
stop("var.mod.obj must be of class, \"variogram.model\".\n")
if(mode==1 && length(at)!=n) stop("length of x and z differ\n")
if(rsearch>0 & nsearch>0)
stop("specify only one of rsearch and nsearch\n")
if(nsmin>nsmax)
stop("nsmin>nsmax\n")
# if(rsearch>0){
# if(nsmin==0) nsmin<-ceiling(n*0.1)
# if(nsmax==0) nsmax<-ceiling(n*0.9)
# }
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)
c0<-0
covmat<-matrix(0,n,n)
# print(c(nx,ny,nz))
ans<-.Fortran("krgtil",
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),
nz=as.integer(nz),
ntx=as.integer(ntx),
nty=as.integer(nty),
nt=as.integer(nt),
dx=as.double(dx),
dy=as.double(dy),
itx=as.integer(itx),
ity=as.integer(ity),
ipt=as.integer(ipt),
xg=double(nx),
yg=double(ny),
zg=double(nz),
varg=double(nz),
dog=as.integer(dog),
xgwork=double(nz),
ygwork=double(nz),
lon=as.double(point.obj$x),
lat=as.double(point.obj$y),
z=as.double(at),
extrap=as.integer(extrap),
n=as.integer(n),
covtype=as.integer(covtype),
covpar=as.double(var.mod.obj$parameters),
cov=as.double(covmat),
c0vec=double(n*ipt),
c0=as.double(c0),
trend=as.integer(trend),
ntrend=as.integer(ntrend),
rsearch=as.double(rsearch),
nsearch=as.integer(nsearch),
nsmin=as.integer(nsmin),
nsmax=as.integer(nsmax),
fwork=double(n*ntrend),
f0work=double(ipt*ntrend),
dist=double(n),
indsnb=integer(n),
indsna=integer(n),
indsrt=integer(n),
kwork=double((n+ntrend)*(n+ntrend)),
nkwork=as.integer(n+ntrend),
rhswork=double(ipt*(n+ntrend)),
ipiv=integer(ipt*(n+ntrend)),
mode=as.integer(mode),
mu=double(ipt*ntrend),
lambda=double(n*ipt),
x0=double(ipt),
y0=double(ipt),
z0=double(ipt),
do0=integer(ipt),
inddo=integer(ipt),
var0=double(ipt),
ierr=integer(1))
# ans<-krige.solve(s$x,s$y,point.obj$x,point.obj$y,
# at,covmat,c0vec,c0,trend,rsearch,nsmin,nsmax,mode)
retval<-list(x=ans$xg,
y=ans$yg,
z=matrix(ans$zg,nx,ny),
var=matrix(ans$varg,nx,ny))
retval
}
bk.grid <- function(xsw,ysw,xne,yne,
dx,dy,
angle=NULL,
point.obj,
at,
var.mod.obj,
maxdist = NULL,
extrap = F,
trend=0,
rsearch=0,
nsearch=0,
nsmin=-1,
nsmax=-1,
mode=1)
{
if(is.null(angle))angle<-0
dgx<-xne-xsw
dgy<-yne-ysw
nx<-ceiling(dgx/dx)
ny<-ceiling(dgy/dy)
at <- point.obj[[match(at, names(point.obj))]]
n<-length(point.obj$x)
nz<-nx*ny
dog<-matrix(1,nx,ny)
if (!extrap) {
dog <- in.convex.hull(tri.mesh(point.obj$x, point.obj$y,
duplicate = "remove"),
point.obj$x, point.obj$y)*dog
dog<-matrix(dog,nx,ny)
}
extrap<-1*extrap
if (!inherits(point.obj, "point"))
stop("point.obj must be of class, \"point\".\n")
if (!inherits(var.mod.obj, "variogram.model"))
stop("var.mod.obj must be of class, \"variogram.model\".\n")
if(mode==1 && length(at)!=n) stop("length of x and z differ\n")
if(rsearch>0 & nsearch>0)
stop("specify only one of rsearch and nsearch\n")
if(nsmin>nsmax)
stop("nsmin>nsmax\n")
# if(rsearch>0){
# if(nsmin==0) nsmin<-ceiling(n*0.1)
# if(nsmax==0) nsmax<-ceiling(n*0.9)
# }
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)
c0<-0
covmat<-matrix(0,n,n)
ans<-.Fortran("krggrd",
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),
nz=as.integer(nz),
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(at),
extrap=as.integer(extrap),
n=as.integer(n),
covtype=as.integer(covtype),
covpar=as.double(var.mod.obj$parameters),
cov=as.double(covmat),
c0vec=double(n),
c0=as.double(c0),
trend=as.integer(trend),
ntrend=as.integer(ntrend),
rsearch=as.double(rsearch),
nsearch=as.integer(nsearch),
nsmin=as.integer(nsmin),
nsmax=as.integer(nsmax),
fwork=double(n*ntrend),
f0work=double(ntrend),
dist=double(n),
indsnb=integer(n),
indsna=integer(n),
indsrt=integer(n),
kwork=double((n+ntrend)*(n+ntrend)),
nkwork=as.integer(n+ntrend),
rhswork=double(n+ntrend),
ipiv=integer(n+ntrend),
mode=as.integer(mode),
mu=double(ntrend),
lambda=double(n),
ierr=integer(1),
.Package="rgeostat")
# ans<-krige.solve(s$x,s$y,point.obj$x,point.obj$y,
# at,covmat,c0vec,c0,trend,rsearch,nsmin,nsmax,mode)
retval<-list(x=ans$xg,
y=ans$yg,
z=matrix(ans$zg,nx,ny),
var=matrix(ans$varg,nx,ny))
retval
}
bk.cell.pts <- function (s,
point.obj,
at,
var.mod.obj,
maxdist = NULL,
extrap = F,
trend=0,
rsearch=0,
nsearch=0,
nsmin=-1,
nsmax=-1,
mode=1)
{
if (!inherits(s, "point"))
stop("s must be of class, \"point\".\n")
if (!inherits(point.obj, "point"))
stop("point.obj must be of class, \"point\".\n")
if (!inherits(var.mod.obj, "variogram.model"))
stop("var.mod.obj must be of class, \"variogram.model\".\n")
s$do <- c(rep(1, length(s$x)))
if (!extrap) {
s$do <- in.convex.hull(tri.mesh(point.obj$x, point.obj$y,
duplicate = "remove"), s$x, s$y)*s$do
}
at <- point.obj[[match(at, names(point.obj))]]
# krige.maxdist(s, point.obj, at, var.mod.obj, maxdist)
n0 <- length(s$x)
c0 <- var.mod.obj$parameters["X1"]+var.mod.obj$parameters["X2"]
xy0.dist <- as.matrix(dist(rbind(cbind(s$x,s$y),
cbind(point.obj$x,
point.obj$y))
,diag=T,upper=T))[-c(1:n0),c(1:n0)]
c0vec <- c0-var.mod.obj$model(xy0.dist,var.mod.obj$parameters)
xx.dist<-as.matrix(dist(cbind(point.obj$x,point.obj$y),
diag=T,upper=T))
covmat <- c0-var.mod.obj$model(xx.dist,var.mod.obj$parameters)
ans<-krige.solve(s$x,s$y,point.obj$x,point.obj$y,
at,covmat,c0vec,c0,trend,s$do,rsearch,nsearch,
nsmin,nsmax,mode)
ans
}
bk.cell <- function (x,
y,
point.obj,
at, var.mod.obj,
maxdist = NULL,
extrap = F,
trend=0,
rsearch=0,
nsearch=0,
nsmin=-1,
nsmax=-1,
mode=1)
{
s<-data.frame(x=x,y=y)
if (!inherits(point.obj, "point"))
stop("point.obj must be of class, \"point\".\n")
if (!inherits(var.mod.obj, "variogram.model"))
stop("var.mod.obj must be of class, \"variogram.model\".\n")
s$do <- c(rep(1, length(s$x)))
if (!extrap) {
s$do <- in.convex.hull(tri.mesh(point.obj$x, point.obj$y,
duplicate = "remove"), s$x, s$y)*s$do
}
at <- point.obj[[match(at, names(point.obj))]]
# krige.maxdist(s, point.obj, at, var.mod.obj, maxdist)
n0 <- length(s$x)
c0 <- var.mod.obj$parameters["X1"]+var.mod.obj$parameters["X2"]
xy0.dist <- as.matrix(dist(rbind(cbind(s$x,s$y),
cbind(point.obj$x,
point.obj$y))
,diag=T,upper=T))[-c(1:n0),c(1:n0)]
c0vec <- c0-var.mod.obj$model(xy0.dist,var.mod.obj$parameters)
xx.dist<-as.matrix(dist(cbind(point.obj$x,point.obj$y),
diag=T,upper=T))
covmat <- c0-var.mod.obj$model(xx.dist,var.mod.obj$parameters)
ans<-krige.solve(s$x,s$y,point.obj$x,point.obj$y,
at,covmat,c0vec,c0,trend,s$do,rsearch,nsearch,
nsmin,nsmax,mode)
ans
}
bk.solve <- function(x0,
y0,
x,
y,
z,
prior,
var.mod.obj,
covmat,
# c0vec=NULL,
# c0=NULL,
trend=0,
do0=NULL,
rsearch=0,
nsearch=0,
nsmin=-1,
nsmax=-1,
mode=1)
{
n<-length(x)
if(length(y)!=n) stop("length of x and y differ\n")
n0<-length(x0)
if(length(y0)!=n0) stop("length of x0 and y0 differ\n")
if(is.null(do0)) do0<-rep(1,n0)
if(length(do0)!=n0) stop("length of x0 and do0 differ\n")
if(mode==1 && length(z)!=n) stop("length of x and z differ\n")
if(dim(covmat)[1]!=n | dim(covmat)[2]!=n)
stop("wrong dimension in covmat\n")
#if(is.vector(c0vec))
# if(length(c0vec)!=n) stop("wrong dimension in c0vec\n")
#if(is.matrix(c0vec))
# if(dim(c0vec)[1]!=n) stop("wrong dimension in c0vec\n")
if(rsearch>0 & nsearch>0)
stop("specify only one of rsearch and nsearch\n")
if(nsmin>nsmax)
stop("nsmin>nsmax\n")
# if(rsearch>0){
# if(nsmin==0) nsmin<-ceiling(n*0.1)
# if(nsmax==0) nsmax<-ceiling(n*0.9)
# }
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)
nadd<-dim(prior$mean)[2]
c0<-0
c0vec<-rep(0,n)
ans<-.Fortran("bk",
lon0=as.double(x0),
lat0=as.double(y0),
do0=as.integer(do0),
inddo=integer(n0),
n0=as.integer(n0),
lon=as.double(x),
lat=as.double(y),
z=as.double(z),
n=as.integer(n),
covtype=as.integer(covtype),
covpar=as.double(var.mod.obj$parameters),
cov=as.double(covmat),
c0vec=as.double(c0vec),
c0=as.double(c0),
trend=as.integer(trend),
mupr=as.double(prior$mean),
phipr=as.double(prior$var),
nadd=as.integer(nadd),
edim=as.integer(prior$dim),
ntsqr=as.integer(ntrend*ntrend),
ntrend=as.integer(ntrend),
rsearch=as.double(rsearch),
nsearch=as.integer(nsearch),
nsmin=as.integer(nsmin),
nsmax=as.integer(nsmax),
fwork=double(n*ntrend),
f0work=double(n0*ntrend),
dist=double(n),
indsnb=integer(n),
indsna=integer(n),
indsrt=integer(n),
kwork=double((n+ntrend)*(n+ntrend)),
nkwork=as.integer(n+ntrend),
rhswork=double(n0*(n+ntrend)),
fpwork=double(n*ntrend),
fpfwork=double(n*n),
fpf0wrk=double(n),
ipiv=integer(n0*(n+ntrend)),
mode=as.integer(mode),
mu=double(n0*ntrend),
z0=double(n0),
lambda=double(n*n0),
lambda0=double(1),
var=double(n0),
ierr=integer(1)
)
if(mode==1)
ret<-list(z0=ans$z0,
lambda=matrix(ans$lambda,n,n0),
var=ans$var,
mu=matrix(ans$mu,ntrend,n0))
if(mode==2)
ret<-list(lambda=matrix(ans$lambda,n,n0),
var=ans$var,
mu=matrix(ans$mu,ntrend,n0))
ret
}
glsfit.solve<-function(formula,x,weight=diag(rep(1,length(x)))){
if(!is.data.frame(x))stop("x is not a data frame!")
n <- dim(x)[1]
p <- dim(x)[2]
if(is.null(formula))stop("formula not given!")
ft<-terms(formula)
vars<-attr(ft,"variables")
fmat<-model.matrix.default(formula, data=x)
if(length(vars)!=4 & length(vars)!=2) stop("need 1 or 3 variables (\"z ~ f(x, y)\" or \"z ~ 1\") in formula!")
if(length(vars)==4){
namx<-as.character(vars[3])
namy<-as.character(vars[4])
}
namz<-as.character(vars[2])
ldc<-dim(weight)[1]
ntrend<-dim(fmat)[2]
.Fortran("glsfit",
fmat=as.double(fmat),
n=as.integer(n),
ntrend=as.integer(ntrend),
ldf=as.integer(n),
y=x[,namz],
weight=as.double(weight),
ldc=as.integer(ldc),
beta=double(ntrend),
u=double(n),
work=double(n+ntrend+n),
lwork=as.integer(n+ntrend+n),
info=integer(1)
)
}
"empirical.prior" <-
function (x, ..., formula=NULL,var.mod=NULL)
{
args <- list(x, ...)
n <- length(args)
# nx <- length(x)
nx<-NULL
for(i in 1:n)
nx<-c(nx,length(args[[i]][,1]))
if(is.null(formula))stop("formula not given!")
ft<-terms(formula)
vars<-attr(ft,"variables")
ntr<-dim(model.matrix.default(formula=formula,data=args[[1]]))[2]
pmean <- matrix(0,ntr,n)
pvar <- matrix(0,ntr*ntr,n)
if(length(vars)!=4 & length(vars)!=2) stop("need 1 or 3 variables (\"z ~ f(x, y)\" or \"z ~ 1\") in formula!")
if(length(vars)==4){
namx<-as.character(vars[3])
namy<-as.character(vars[4])
}
namz<-as.character(vars[2])
for (i in 1:n) {
if(!is.null(var.mod)){
c0 <- var.mod$parameters["X1"]+var.mod$parameters["X2"]
xx.dist<-as.matrix(dist(cbind(args[[i]][,namx], args[[i]][,namy])),
diag=T,upper=T)
covmat <- c0-var.mod$model(xx.dist,var.mod$parameters)
}
# xt.lm <- lm(formula,xtrafo)
# x.glsfit<-glsfit.solve(formula,as.data.frame(args[[i]]),weight=covmat)
# pmean <- c(pmean, x.glsfit$beta)
# pvar <- c(pvar, x.glsfit$vbeta)
#nur ok, kein gls
pmean[,i] <- c(mean(args[[i]][,namz]))
pvar[,i] <- c(matrix(var(args[[i]][,namz])/nx[i],1,1))
}
list(mean = pmean, var = pvar, dim=nx)
}
"subjective.prior" <-
function (x, bounds)
{
}
R/zzz.R 0 → 100644
View file @ b022db9c
.First.lib <- function(lib, pkg) {
library.dynam("baykrig", pkg, lib)
require(sgeostat)
}
if(version$minor < "62") {
library.dynam("baykrig")
require(sgeostat)
}
src/bk.f 0 → 100644
View file @ b022db9c
SUBROUTINE BKTILE(XSW,YSW,XNE,YNE,ANGLE,NX,NY,NZ,NTX,NTY,NT,
. DX,DY,ITX,ITY,IPT,XG,YG,ZG,VARG,DOG,
. XGWORK,YGWORK,
. LON,LAT,Z,EXTRAP,N,
. COVTYPE,COVPAR,COV,C0,COV0,TREND,
. NTREND,RSEARCH,NSEARCH,NSMIN,NSMAX,FWORK,F0WORK,
. DIST,INDSNB,INDSNA,INDSRT,KWORK,NKWORK,RHSWORK,
. IPIV,MODE,MU,LAMBDA,X0,Y0,Z0,DO0,inddo,VAR0,
. IERR)
INTEGER NX,NY,NZ,NTX,NTY,NT,ITX,ITY,N,COVTYPE,TREND,NTREND,
. NSEARCH,NSMIN,NSMAX,MODE,IERR,INDSNB(*),INDSNA(*),
. INDSRT(*),IPIV(*),IPT,EXTRAP,DOG(NX,*),DO0(*),inddo(*)
DOUBLE PRECISION XSW,YSW,XNE,YNE,ANGLE,DX,DY,
. XG(*),YG(*),ZG(NX,*),VARG(NX,*),
. LON(*),LAT(*),Z(*),COV(N,*),C0(*),COV0(N,*),
. RSEARCH,FWORK(N,*),F0WORK(NTREND,*),DIST(*),
. KWORK(NKWORK,*),RHSWORK(NKWORK,*),MU(NTREND,*),
. LAMBDA(N,*),X0(*),Y0(*),Z0(*),VAR0(*),
. XGWORK(NX,*),YGWORK(NX,*),COVPAR(*)
c subroutine for kriging prediction on tiles of a grid
c
c This subroutine takes the grid specification, forms tiles
c (i.e. rectangular subregions of the grid) and calls KRIGE
c (see below) on these tiles. The idea is to reduce computational
c burden by collecting some neigbouring "krige"-systems and
c forming a combined krige system with multiple right hand sides
c to be solved by DGESV simultaneously.
c
c
c
c [1,1] ne
c +-----------+----------+--------+
c | o o o | o o o| o o | o -- grid points (i,j)
c |(1,1) | | x |
c | o o o | o o o| o o | +---+
c | | | | | o | -- grid tiles [i,j]
c | o o o | o o o| o o | +---+
c +-----------+----------+--------+
c 1 | o o o | o o o| o o | x -- user specified