R/f.polygons.preCKrige.R

In constrainedKriging: Constrained, Covariance-Matching Constrained and Universal Point or Block Kriging

###############################################################
#                                                             #
#   Berechnungen der Polygonapproximationen durch die Pixel   #
#   ch 22-02-2010                                             #
#                                                             #
###############################################################
f.polygons.preCKrige <- function(
    				newdata,
				neighbours,
				model,
				pwidth,
				pheight,
				napp = 1
			    )
{
#
# no attributes -> only for ordinary kriging
# add an empty data.frame to the Polygons
# and build Class SpatialPolygonsDataFrame
if( class( newdata ) == "SpatialPolygons" )
{
    data = as.data.frame( matrix( ncol = 0, nrow = 0 ) )
}
#
if( class( newdata ) == "SpatialPolygonsDataFrame" )
{
    data = newdata@data
}
#
class(model) <- "list"
### covmodel of measurement free error process
model.me.free <- model[unlist(lapply(1:length(model), function(i,m){m[[i]]$model != "mev"},m = model))]
### Polygons as calss SpatialPolygons
newdata.polygons = SpatialPolygons( newdata@polygons )
#
# Ziel polygone haben keien Nachbarn
if( missing( neighbours ) )
{
    neighbours <- lapply( as.list( 1:length(newdata@polygons) ), function( x ){ return( integer(0) ) } )
}
stopifnot( is.list(neighbours) )
#
#
pixgrid<- f.pixelgrid(
    			polygons = newdata.polygons,
			neighbours = neighbours,
    			pixel.x.width = pwidth,
			pixel.y.width = pheight
		    )
#
## Kovarianzmatrix der Pixel
pixcm <- f.pixelcovmat(
              		pixgrid = pixgrid,
              		model = model.me.free
		    )

#
t.n.poly <- as.list( 1:length( newdata.polygons@polygons ) )
for( i in 1:napp ){
# pc.tmp = pixconfig.temporary
pc.tmp <- f.pixconfig(
                    		polygons = newdata.polygons,
                                neighbours = neighbours,
                                pixgrid = pixgrid,
                                n = napp
			    )
if( i == 1 )
{
    pixconfig = pc.tmp
}
else
{
t2 <- proc.time()[3]
pixconfig <- lapply( pixconfig,
    function( pc, pc.tmp )
    {
	pc$pixcenter <- cbind( pc$pixcenter,
	    pc.tmp[[ pc$posindex[1] ]]$pixcenter )
	pc$pix.in.poly <- cbind( pc$pix.in.poly, pc.tmp[[ pc$posindex[1] ]]$pix.in.poly )
	pc$sa.polygons  <- c(pc$sa.polygons,
	    pc.tmp[[ pc$posindex[1] ]]$sa.polygons)
	return(pc)
    },
    pc.tmp= pc.tmp
    )
}
rm( pc.tmp )
}
#
#stopifnot( 1 == 0 )
cm.list <- f.polygoncovmat(
    		pixconfig = pixconfig,
		pixcm = pixcm,
		model = model.me.free,
		n = napp
	    )

#
return( new( "preCKrigePolygons",
	covmat = cm.list$mean.bb.cov.mat, # block block Covariance Matrix
	se.covmat = cm.list$var.mean.bb.cov.mat, # standard error of block block Covariance
	pixconfig = pixconfig,
	pixcovmat = pixcm,
    	model = model,
	data = data,
	polygons = newdata@polygons
    )
)
rm( pixcm, pixgrid, newdata.polygons, pixconfig, cm.list)
}