## Authors ## Martin Schlather, schlather@math.uni-mannheim.de ## ## ## Copyright (C) 2012 -- 2014 Alexander Malinowski & Martin Schlather ## 2015 -- 2017 Martin Schlather ## ## This program is free software; you can redistribute it and/or ## modify it under the terms of the GNU General Public License ## as published by the Free Software Foundation; either version 3 ## of the License, or (at your option) any later version. ## ## This program is distributed in the hope that it will be useful, ## but WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with this program; if not, write to the Free Software ## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. ## Methods for classes 'RFempVario' and 'RFfit' ####################### summary_RFempVariog <- function(OP, object, ...) { variogram <- do.call(OP, list(object, "empirical")) if (is.array(variogram)) { dims <- dim(variogram) dims <- dims[dims > 1] dim(variogram) <- dims } l <- list(centers=do.call(OP, list(object, "centers")), variogram=variogram) obj <- do.call(OP, list(object, "sd")) if (length(obj) > 0) { l$sd <- obj if (is.array(variogram)) dim(l$sd) <- dims } obj <- do.call(OP, list(object, "phi.centers")) if (length(obj) > 1) l$phi <- obj obj <- do.call(OP, list(object, "theta.centers")) if (length(obj) > 1) l$theta <- obj obj <- do.call(OP, list(object, "T")) if (length(obj)>0) l$T <- obj obj <- do.call(OP, list(object, "vdim")) if (obj > 1) l$vdim <- obj class(l) <- "summary.RFempVariog" l } summary.RFempVariog <- function(object, ...) summary_RFempVariog("@",object,...) summary.RF_empVariog <- function(object,...) summary_RFempVariog("$",object,...) print.summary.RFempVariog <- function(x, ...) { cat("Object of class 'RFempVariog'\n") str(x, no.list=TRUE) # invisible(x) } print.RFempVariog <- function(x, ...) { print.summary.RFempVariog(summary.RFempVariog(x, ...)) } print.RF_empVariog <- function(x, ...) { print.summary.RFempVariog(summary.RF_empVariog(x, ...)) } setMethod(f="show", signature="RFempVariog", definition=function(object) print.RFempVariog(object)) # ## coercion methods setAs("RFempVariog", "list", def=function(from) { li <- list() for (name in slotNames(from)) { li[[name]] <- eval(parse(text=paste("from@", name, sep=""))) } return(li) }) setAs("RFfit", "RFempVariog", def=function(from) from@ev) list2RFempVariog <- function(li) { RFopt.coords <- RFoptions()$coords return(new("RFempVariog", centers=li$centers, empirical=li$empirical, var=var, sd=li$sd, n.bin=li$n.bin, phi.centers=li$phi.centers, theta.centers=li$theta.centers, T=li$Tbins, vdim = li$vdim, coordunits = RFopt.coords$coordunits, varunits = RFopt.coords$varunits, call=li$call, method=li$method)) } ## plot method plotRFempVariogUnbinned <- function(x, coordunits, varunits, varnames, legend, ..., plotmethod="plot") { stop("currently not used ") # !! graphics <- RFoptions()$graphics dots = mergeWithGlobal(list(...)) dotnames <- names(dots) coords <- GridTopology2gridVectors(cbind(x@centers$x, x@centers$T)) if (length(coords)>1) { coords[[1]] <- sort(unique((coords[[1]] - min(coords[[1]])) * rep(c(-1, 1), each=length(coords[[1]])))) lab.names <- dimnames(x@centers$x)[[2]] if (!(x@centers$spatialdim==1 && x@centers$has.time.comp) && !(x@centers$x[3,2]==dim(x@empirical)[2])) coords[[2]] <- sort(unique((coords[[2]] - min(coords[[2]])) * rep(c(-1, 1), each=length(coords[[2]])))) else lab.names[2] <- "T" } else { x@centers <- coords[[1]]-min(coords[[1]]) do.call(graphics::plot, args=c(dots, list(x=x, plot.nbin=FALSE))) return() } if (!("main" %in% dotnames)) { main <- "Variogram image plot" if (length(varnames)>0) main <- paste(main, "for", varnames) dots$main <- main } lab.names <- paste(lab.names, "-distance", sep="") idx <- lab.names != "T-distance" if (any(idx) && all(coordunits[idx] != "")) lab.names[idx] <- paste(lab.names[idx], " [", coordunits[idx], "]", sep="") if (!all(idx) && all(coordunits[!idx] != "")) lab.names[!idx] <- paste(lab.names[!idx], " [", coordunits[!idx], "]", sep="") if (!("xlab" %in% dotnames)) dots$xlab <- lab.names[1] if (!("ylab" %in% dotnames)) dots$ylab <- lab.names[2] if (!("xlim" %in% dotnames)) dots$xlim <- range(coords[[1]]) if (!("ylim" %in% dotnames)) dots$ylim <- range(coords[[2]]) idx1 <- coords[[1]] >= dots$xlim[1] & coords[[1]] <= dots$xlim[2] idx2 <- coords[[2]] >= dots$ylim[1] & coords[[2]] <= dots$ylim[2] coords[[1]] <- coords[[1]][idx1] coords[[2]] <- coords[[2]][idx2] dims <- dim(x@empirical) ev.plot <- matrix(x@empirical, nrow=dims[1], ncol=dims[2]) ev.plot <- ev.plot[idx1, idx2] range.ev <- range(ev.plot) #x@empirical col <- (if ("col" %in% dotnames) dots$col else default.image.par(NULL, NULL)$default.col) if (graphics$split.screen && legend) { scr <- split.screen( rbind(c(0,.85,0,1), c(.85,1,0,1))) screen(scr[1]) par(mar=c(4,4,3,1)) } else { scr <- NULL par(mar=c(4,4,3,1)) } dots$type <- NULL do.call(plotmethod, args=c(dots, list(x=coords[[1]], y=coords[[2]], z=ev.plot))) # xlab=lab.names[1], # ylab=lab.names[2], # main="Variogram image plot" if (legend) { stop("das ist doch quark?!") if (graphics$split.screen) { screen(scr[2]) par(mar=c(4,0,3,3)) z.legend <- seq(range.ev[1], range.ev[2], length=length(col)) image(y=z.legend, x=c(0,1), z=rbind(z.legend, z.legend), axes=F, col=col, xlab="") axis(4) box() } else { my.legend(min(coords[[1]]), max(coords[[2]]), range(ev.plot), col=col, bg="white") } } if (graphics$close_screen) { close.screen(scr) scr <- NULL } return(invisible(scr)) } RFplotEmpVariogram <- function(x, model = NULL, nmax.phi = NA, nmax.theta = NA, nmax.T = NA, plot.nbin = TRUE, plot.sd=FALSE, method = "ml", variogram=TRUE, boundaries = TRUE, ...) { if (!variogram) stop("plot of estimated covariance functions not programmed yet.") graphics <- RFoptions()$graphics newx <- list() methodnames <- double() if (is(x, "RFfit") || is(x, "RF_fit")) { OP <- c("$", "@")[1 + is(x, "RFfit")] if(length(do.call(OP, list(x, "ev")))==0) stop("The fit does not contain an empirical variogram.") newx$autostart <- do.call(OP, list(x, "autostart")) newx$self <- do.call(OP, list(x, "self")) newx$plain <- do.call(OP, list(x, "plain")) newx$sqrt.nr <- do.call(OP, list(x, "sqrt.nr")) newx$sd.inv <- do.call(OP, list(x, "sd.inv")) newx$internal1 <- do.call(OP, list(x, "internal1")) newx$internal2 <- do.call(OP, list(x, "internal2")) newx$internal3 <- do.call(OP, list(x, "internal3")) newx$ml <- do.call(OP, list(x, "ml")) x <- do.call(OP, list(x, "ev")) if(is.null(method)) method <- if (length(newx$ml@name)>0) "ml" else "plain" methodidx<-names(newx[unlist(lapply(newx, function(x) is(x, CLASS_FIT)))]) methodidx <- method %in% methodidx methodnames <- method[methodidx] nomethodnames <- method[!methodidx] if(!all(methodidx)) warning( paste("The following method does not exist: ", nomethodnames)) } else { if (!is(x, "RFempVariog") && !is(x, "RF_empVariog")) stop("method only for objects of class 'RFempVariog' or 'RFfit'") OP <- c("$", "@")[1 + is(x, "RFempVariog")] } ev <- do.call(OP, list(x, "empirical")) sd <- do.call(OP, list(x, "sd")) centers <- do.call(OP, list(x, "centers")) phi.centers <- do.call(OP, list(x, "phi.centers")) theta.centers <- do.call(OP, list(x, "theta.centers")) Time <- do.call(OP, list(x, "T")) n.bin <- do.call(OP, list(x, "n.bin")) coordunits <- do.call(OP, list(x, "coordunits")) empmethod <- do.call(OP, list(x, "method")) varnames <- if (is.matrix(ev)) dimnames(ev)[[2]][1] else names(ev)[1] if(empmethod > 2 && !is.null(model)) warning("theoretical madogram model is not available at the moment") dots = list(...) dotnames <- names(dots) if (!("type" %in% dotnames)) dots$type <- "b" cex <- if ("cex" %in% dotnames) dots$cex else .8 par(cex=cex, xaxs="i") dots$cex <- NULL if (!("pch" %in% dotnames)) dots$pch <- 19 if(!("xlim" %in% dotnames)) dots$xlim <- range(centers) ylim.not.in.dotnames <- !("ylim" %in% dotnames) xlab <- if ("xlab" %in% dotnames) dots$xlab else "distance" dots$xlab <- NULL ylab.ev <- if ("ylab" %in% dotnames) { dots$ylab } else if(empmethod == 2) { "covariance" } else if(empmethod < 2) { "semivariance" } else { "madogram" } dots$ylab <- NULL main0 <- if ("main" %in% dotnames) { dots$main } else if(empmethod == 2) { "Covariance plot" } else if(empmethod < 2) { "Variogram plot" } else { "Madogram plot" } dots$main <- NULL if (!is.null(main0)) oma.top <- 2 else oma.top <- 0 has.sd <- !is.null(sd) if(!is.null(model)) { if (!is.list(model)) model <- list(model) if (!all(unlist(lapply(model, FUN=function(x) is(x, CLASS_CLIST))))) stop("model must be (a list of elements) of class 'CLASS_CLIST'") modelnames <- if(length(names(model)) > 0) names(model) else paste("model", 1:length(model)) methodnames <- c(methodnames, modelnames) names(model) <- modelnames newx <- c(newx, model) } n.methods <- length(methodnames) n.phi <- min(nmax.phi, l.phi <- max(1,length(phi.centers)), na.rm=TRUE) n.theta <- min(nmax.theta, l.theta <- max(1,length(theta.centers)), na.rm=TRUE) n.T <- min(nmax.T, l.T <- max(1,length(Time)), na.rm=TRUE) vdim <- dim(ev) vdim <- vdim[length(vdim)] if (n.phi > 6 || n.theta > 3 || n.T > 3) message("'If you feel the picture is overloaded, set the parameters 'nmax.phi', 'nmax.theta' and 'nmax.T'") halfphi.sector <- pi/(2*l.phi) halftheta.sector <- pi/(2*l.theta) phi.angles <- c(halfphi.sector, 0, -halfphi.sector) theta.angles <- seq(-halftheta.sector, halftheta.sector, len=5) # len=odd! if (n.phi>1 && boundaries) { phi.angles <- phi.angles * 0.96 theta.angles <- theta.angles * 0.96 } TandV <- (n.T > 1 && vdim > 1) && graphics$split_screen if (vdim>1 && length(varnames)==0) varnames <- paste("v", 1:vdim, sep="") range.nbin <- range(c(0, n.bin), na.rm=TRUE) ylim.nbin <- range.nbin * c(1,1.2) col.v <- col <- if ("col" %in% dotnames) rep(dots$col, len=n.phi) else 1:max(n.phi) dots$col <- NULL if (n.methods > 0){ dotsRFfit <- dots dotsRFfit$type <- "l" dotsRFfit$lwd <- 2 ltyRFfit.v <- 1:n.methods dotsRFfit$lty <- NULL } dir2vario <- function(dir.vec, x.eval, x.time, method.model, v1, v2){ x.space <- as.matrix(x.eval) %*% t(dir.vec) for (j in (1:4)) { ## orig 20 if(empmethod == 2) { # covariance plot vario.vals <- try(RFcov(x = cbind(x.space, x.time), model = method.model, grid = FALSE, internal.examples_reduced=FALSE), silent = TRUE) } else if(empmethod < 2) { vario.vals <- try(RFvariogram(x = cbind(x.space, x.time), model = method.model, grid = FALSE, internal.examples_reduced=FALSE), silent = TRUE) } else { break; } if(!is(vario.vals, "try-error")) { if (is.array(vario.vals)) { return(vario.vals[, v1, v2]) } else { return(vario.vals) } } x.space <- cbind(x.space, 0) } return(NA) } oma.left <- 6 Screens <- if (TandV) c(n.T, n.theta) else c(n.T * vdim * vdim, n.theta) ArrangeDevice(graphics, Screens) all.scr <- scr <- (if (!graphics$split_screen) NULL else split.screen(if (TandV) c(vdim, vdim) else Screens)) for (v1 in 1:vdim) { for (v2 in 1:vdim) { if (TandV) { scr <- split.screen(Screens, all.scr[v1 + (v2 - 1) * vdim]) all.scr <- c(all.scr, scr) } if (vdim == 1) { main <- if (!is.null(main0) && length(varnames)>0) paste(main0, "for", varnames) else main0 } else { main <- if (!TandV) main0 else paste(main0, "for", varnames[v1], "vs.", varnames[v2]) } if (ylim.not.in.dotnames) dots$ylim <- range(ev[,,,, v1, v2], na.rm=TRUE) for (iT in 1:n.T) { for (ith in 1:n.theta) { ## plot n.bin if (plot.nbin) { screen(scr[1]) par(oma=c(4,oma.left,oma.top,0)) scr2 <- split.screen(rbind(c(0,1,.2,1), c(0,1,0,.2)), screen=scr[1]) all.scr <- c(all.scr, scr2) screen(scr2[2]) par(mar=c(0,.5,0,.5)) for (iph in 1:n.phi) { if (n.phi > 1) col <- col.v[iph] if (iph==1) { lab <- xylabs("bin centers", NULL, units=coordunits) plot(centers, n.bin[ ,iph, ith, iT, v1, v2], xlim=dots$xlim, ylim=ylim.nbin, type=if (n.phi>1) "p" else "h", col =if (n.phi>1) col else "darkgray", lwd=4, pch=16, axes=FALSE, ylab="", xlab = lab$x) box() at <- seq(range.nbin[1], range.nbin[2], len=3) if (ith==1) axis(2, las=1, at=at, labels=format(at, scient=-1, dig=2), outer=TRUE) else axis(2, las=1, at=at, labels=FALSE) if (iT==n.T && (n.T > 1 || (v1==vdim && v2==vdim))) axis(1) if (ith==1) title(ylab="n.bin", line=5, outer=TRUE, adj=0) } else { points(centers, n.bin[ ,iph, ith, iT, v1, v2], type="p", col=col, pch=16) } } screen(scr2[1]) } else { screen(scr[1]) par(oma=c(4,oma.left,oma.top,0)) } ## plot empirical ##if (ith==1) par(mar=c(0,6,1,1)) else par(mar=c(0,.5,1,.5)) plotted.meth <- NULL # needed for legend for (iph in 1:n.phi) { if (n.phi>1) col <- col.v[iph] if (iph==1) { do.call(graphics::plot, c(dots, list( x=centers, y=ev[,iph,ith,iT,v1,v2], #ylim=ylim.ev, type=type, pch=19, col=col, axes=FALSE, ylab="")) ) box() axis(2, las=1, labels=(ith==1), outer=(ith==1)) if (!plot.nbin) axis(1) if (l.theta > 1 || l.T > 1 || vdim > 1) { L <- character(3) if (!TandV && vdim > 1) L[1] <- paste(varnames[c(v1,v2)], collapse=":") if (l.T>1) L[2] <- paste("T=",signif(Time[iT], 3), " ", sep="") if (l.theta>1) L[3] <- paste(sep="", "theta=", signif(theta.centers[ith],3)) legend("topleft", legend=paste(L[L!=""], collapse=",")) } if (ith == 1) title(ylab=ylab.ev, line=5, outer=TRUE) if (has.sd && plot.sd) legend("topright", bty="n", #col=NA, lwd=NA, legend="arrows represent +-1 sd intervals") } else { ## iph > 1 do.call(graphics::points, c(dots, list(x=centers, y=ev[,iph,ith,iT,v1,v2], col=col)) ) #type=type, pch=19 } ## for iph k <- 1 if (n.methods > 0) { for(i in 1:n.methods) { method.model <- newx[[methodnames[i]]] if(length(method.model@name) == 0){ warning("The method '", methodnames[i], "' cannot be fitted.") next } if (!is.null(phi.centers)) { x.radial <- cbind(cos(phi.centers[iph]+phi.angles), sin(phi.centers[iph]+phi.angles)) if (!is.null(phi.centers)) x.radial <- cbind(x.radial, cos(rep(theta.angles, each = length(phi.angles)) + theta.centers[ith])) } else x.radial <- matrix(1, nrow=1, ncol=1) x.time <- NULL if (!is.null(Time)) x.time <- Time[iT] x.eval <- seq(from = max(dotsRFfit$xlim[1],1e-3), to = dotsRFfit$xlim[2], len = 150) ## sehr genau Abschaetzung, indem mehrere (3) ## Winkel angeschaut werden und dann der mittlere ## Wert der Variogramme angeschaut wird, da ja auch ## das emp. Variogramm ueber ein Winkelintervall gemittelt wird dummy.vals <- as.matrix(apply(x.radial, 1, FUN = dir2vario, x.eval=x.eval, x.time=x.time, method.model=method.model, v1=v1, v2=v2)) ## Print(i, x.eval, dummy.vals);print.RMmodelFit(method.model) if(length(dummy.vals) > 0 && !all(is.na(dummy.vals))){ # Print(n.phi, boundaries) if (n.phi>1 && boundaries) { do.call(graphics::matplot, args=c(dotsRFfit, list(x=x.eval, y= if (ncol(dummy.vals) == 1) dummy.vals else t(apply(dummy.vals, 1, range)), add = TRUE, col=col, lty = 3))) } else { do.call(graphics::points, args=c(dotsRFfit, list( x=x.eval, y = rowMeans(dummy.vals) , col=col, lty = ltyRFfit.v[k]))) } k <- k+1 if(iph == 1) plotted.meth <- c(plotted.meth, methodnames[i]) } } ## for i in n.methods } ## if n.methods > 0 if (has.sd && plot.sd) { sdnot0 <- sd[ ,iph, ith, iT] != 0 arrows(centers[sdnot0], ev[sdnot0 ,iph, ith,iT] - sd[sdnot0,iph,ith,iT], centers[sdnot0], ev[sdnot0 ,iph, ith,iT] + sd[sdnot0,iph,ith,iT], code=2, angle=90, length=0.05, col=col) } } # iph in nphi pifactor <- signif((phi.centers[1:n.phi]) / pi, 2) len.mnames <- length(plotted.meth) string.emp <- "empirical" if(len.mnames > 0) { labels <- if (n.phi>1) paste("\"phi=\",", rep(pifactor, each = len.mnames+1), ", pi, ", "\", ") else "\"" labels <- paste("c(", paste("expression(paste(",labels, rep(c(string.emp, plotted.meth), l.phi), "\"", "))", collapse=","), ")") labels <- eval(parse(text=labels)) } else { labels <- if (n.phi > 1) eval(parse(text=paste( "c(", paste("expression(paste(\"phi=\",", pifactor,", pi))", collapse=","), ")"))) } if (l.phi > 1 || len.mnames > 0)# && iT==1 && ith==1 # auskommentiert auf Sebs wunsch legend("bottomright", col=rep(col.v, each = len.mnames+1), lwd=1, pch=rep(c(19,rep(NA, len.mnames)), l.phi), bty="n", legend=labels, lty = rep(c(1, if(len.mnames==0) NULL else ltyRFfit.v[1:len.mnames]), l.phi)) scr <- scr[-1] } # n.theta } # T } # vdim 1 } # vdim 2 dots$type <- NULL if (!is.null(main)) do.call(graphics::title, args=c(dots, main=main, outer=TRUE)) if (!is.null(xlab)) do.call(graphics::title, args=c(dots, xlab=xlab, outer=TRUE)) if (graphics$close_screen) { close.screen(all.scr) all.scr <- NULL } return(invisible(all.scr)) } setMethod(f="plot", signature(x="RFempVariog", y="missing"), function(x, y, ...) RFplotEmpVariogram(x, ...))