R/Vgrams.R

Defines functions plot.variographied print.variographied variographier.SpatialVx variographier.default variographier plot.structurogram structurogram plot.structurogram.matrix structurogram.matrix summary.lossdiff print.lossdiff plot.lossdiff ORSS print.flossdiff.expvg predict.flossdiff.expvg expvg flossdiff empiricalVG.lossdiff lossdiff.default lossdiff.SpatialVx lossdiff distmaploss abserrloss sqerrloss corrskill plot.griddedVgram griddedVgram

Documented in abserrloss corrskill distmaploss empiricalVG.lossdiff expvg flossdiff griddedVgram lossdiff lossdiff.default lossdiff.SpatialVx ORSS plot.griddedVgram plot.lossdiff plot.structurogram plot.structurogram.matrix plot.variographied predict.flossdiff.expvg print.flossdiff.expvg print.lossdiff print.variographied sqerrloss structurogram structurogram.matrix summary.lossdiff variographier variographier.default variographier.SpatialVx

griddedVgram <- function(object, zero.in=TRUE, zero.out=TRUE, time.point=1, obs = 1, model=1, ...) {

    a <- attributes(object)
    out <- list()
    if(!is.null(a$names)) a$names <- NULL
    attributes(out) <- a

    ## Begin: Get the data sets
    dat <- datagrabber(object, time.point=time.point, obs = obs, model=model)

    X <- dat$X
    Y <- dat$Xhat
    ## End: Get the data sets

    out$zero.in <- zero.in
    out$zero.out <- zero.out

    out$Vx.vgram.matrix <- list()
    out$Fcst.vgram.matrix <- list()

    if(zero.in) {
	out$Vx.vgram.matrix[[1]] <- vgram.matrix(dat=X, ...)
	out$Fcst.vgram.matrix[[1]] <- vgram.matrix(dat=Y, ...)
    } else {
	out$Vx.vgram.matrix[[1]] <- NULL
	out$Fcst.vgram.matrix[[1]] <- NULL
    }

    if(zero.out) {
	out$Vx.vgram.matrix[[2]] <- variogram.matrix(dat=X, zero.out=TRUE, ...)
	out$Fcst.vgram.matrix[[2]] <- variogram.matrix(dat=Y, zero.out=TRUE, ...)
    } else {
	out$Vx.vgram.matrix[[2]] <- NULL
	out$Fcst.vgram.matrix[[2]] <- NULL
    }

    attr(out, "time.point") <- time.point
    attr(out, "obs") <- obs
    attr(out, "model") <- model

    attr(out, "obs.name") <- a$obs.name[ obs ]
    attr(out, "model.name") <- a$model.name[ model ]

    class(out) <- "griddedVgram"

    return(out)

} # end of 'griddedVgram' function.

plot.griddedVgram <- function( x, ... ) {

    a <- attributes(x)

    mainX <- a$obs.name
    mainY <- a$model.name

   if(x$zero.in) {
      vgX <- x$Vx.vgram.matrix[[1]]
      vgY <- x$Fcst.vgram.matrix[[1]]
   }

   if(x$zero.out) {
      vgX.zero <- x$Vx.vgram.matrix[[2]]
      vgY.zero <- x$Fcst.vgram.matrix[[2]]
   }

    if( !is.null( a$msg ) ) {

        op <- par()
        par( oma = c(0, 0, 2, 0) )

    }

   if(x$zero.in) {

      plot(vgX$d, vgX$vgram, xlab="separation distance", ylab="variogram", main=mainX, col="darkblue")
      points(vgX$d.full, vgX$vgram.full, pch=".", cex=1.25, col="darkblue")

      plot(vgY$d, vgY$vgram, xlab="separation distance", ylab="variogram", main=mainY, col="darkblue") 
      points(vgY$d.full, vgY$vgram.full, pch=".", cex=1.25, col="darkblue")

      class(vgX) <- "vgram.matrix"
      class(vgY) <- "vgram.matrix"

      zl <- range(c(c(vgX$vgram.full), c(vgY$vgram.full)), finite=TRUE)
      plot(vgX, xlab="x separations", ylab="y separations", zlim=zl, ...)
      plot(vgY, xlab="x separations", ylab="y separations", zlim=zl, ...)

   }
   if(x$zero.out) {
      plot(vgX.zero$d, vgX.zero$vgram, xlab="separation distance", ylab="variogram", main=paste(mainX, " (non-zeros only)", sep=""), col="darkblue")
      points(vgX.zero$d.full, vgX.zero$vgram.full, pch=".", cex=1.25, col="darkblue")
      plot(vgY.zero$d, vgY.zero$vgram, xlab="separation distance", ylab="variogram", main=paste(mainY, " (non-zeros only)", sep=""), col="darkblue")
      points(vgY.zero$d.full, vgY.zero$vgram.full, pch=".", cex=1.25, col="darkblue")
      class(vgX.zero) <- "vgram.matrix"
      class(vgY.zero) <- "vgram.matrix"
      zl <- range(c(c(vgX.zero$vgram.full), c(vgY.zero$vgram.full)), finite=TRUE)
      plot(vgX.zero, xlab="x separations", ylab="y separations", zlim=zl, ...)
      plot(vgY.zero, xlab="x separations", ylab="y separations", zlim=zl, ...)
   }

    if(!is.null(a$msg)) {

	title("")
	mtext(a$msg, line=0.05, outer=TRUE)

	par( oma = op$oma )

    }

   invisible()

} # end of 'plot.griddedVgram' function.

corrskill <- function(x,y,...) {
   good <- !is.na(x) & !is.na(y)
   n <- sum(good, na.rm=TRUE)
   s1 <- sd(c(x), ...)
   s2 <- sd(c(y), ...)
   m1 <- sum(colSums(x,na.rm=TRUE),na.rm=TRUE)/n
   m2 <- sum(colSums(y,na.rm=TRUE),na.rm=TRUE)/n
   out <- ((n/(n-1))/(s1*s2))*(x-m1)*(y-m2)
   return(out)
} # end of 'corrskill' function.

sqerrloss <- function(x,y,...) return((x-y)^2)
abserrloss <- function(x,y,...) return(abs(x-y))
distmaploss <- function(x,y, threshold=0, const=Inf, ...) {
   if(length(threshold)==1) threshold <- rep(threshold,2)
   xdim <- dim(x)
   if(is.finite(const)) {
	x <- cbind(matrix(0, xdim[1], const), x, matrix(0, xdim[1], const))
	x <- rbind(matrix(0, const, xdim[2]+2*const), x, matrix(0, const, xdim[2] + 2*const))
	y <- cbind(matrix(0, xdim[1], const), y, matrix(0, xdim[1], const))
        y <- rbind(matrix(0, const, xdim[2]+2*const), y, matrix(0, const, xdim[2] + 2*const))
   }
   x <- im(x)
   y <- im(y)
   x <- solutionset(x >= threshold[2])
   y <- solutionset(y >= threshold[1])
   bb <- boundingbox(as.rectangle(x), as.rectangle(y))
   x <- rebound(x, bb)
   y <- rebound(y, bb)
   dy <- distmap(y, ...)$v
   dx <- distmap(x, ...)$v
   if(is.finite(const)) {
	dy[dy>const] <- const
	dx[dx>const] <- const
	dy <- dy[(const+1):(xdim[1]+const),(const+1):(xdim[2]+const)]
	dx <- dx[(const+1):(xdim[1]+const),(const+1):(xdim[2]+const)]
   }
   return(abs(dx - dy))
} # end of 'distmaploss' function.

lossdiff <- function(x, ...) {

    UseMethod("lossdiff", x)

} # end of 'lossdiff' function.

lossdiff.SpatialVx <- function(x, ..., time.point = 1, obs = 1, model = c(1,2), threshold = NULL,
    lossfun = "corrskill", zero.out = FALSE) {

    theCall <- match.call()

    a <- attributes(x)

    # if( missing( maxrad ) ) maxrad <- floor( max( a$xdim ) / 2 )

    # if( is.null( maxrad ) || is.na( maxrad ) || !is.finite( maxrad ) || !is.numeric( maxrad ) ) maxrad <- 20

    ## Begin: Get the data sets
    dat <- datagrabber(x, time.point=time.point, obs = obs, model=model[1])
   X <- dat$X
    Xhat <- dat$Xhat

    dat <- datagrabber(x, time.point=time.point, obs = obs, model=model[2])
    Xhat2 <- dat$Xhat
    ## End: Get the data sets

    out <- lossdiff.default(x = X, ..., xhat1 = Xhat, xhat2 = Xhat2, threshold = threshold, lossfun = lossfun, 
			    loc = a$loc, zero.out = zero.out)

    out$data.name <- c( a$obs.name[ obs ], a$model.name[ model ] )
    out$call <- theCall

    attr(out, "time.point") <- time.point
    attr( out, "obs" ) <- obs
    attr(out, "model") <- model

    attr(out, "msg") <- a$msg
    attr(out, "projection") <- a$projection
    attr(out, "map") <- a$map
    attr(out, "loc") <- a$loc
    attr(out, "loc.byrow") <- a$loc.byrow
    attr(out, "xdim") <- a$xdim

    return(out)

} # end of 'lossdiff.SpatialVx' function.

lossdiff.default <- function(x, ..., xhat1, xhat2, threshold = NULL, lossfun = "corrskill", loc = NULL, zero.out = FALSE) {

    theCall <- match.call()

    out <- list()
    attributes(out) <- atmp <- list(...)
    if(is.null(atmp$loc.byrow)) attr(out, "loc.byrow") <- FALSE

    data.name <- c(as.character(substitute(x)),as.character(substitute(xhat1)),as.character(substitute(xhat2)))
    names(data.name) <- c("verification","model1", "model2")
    out$data.name <- data.name
    xdim <- dim(x)

    out$maxdim <- floor( sqrt(xdim[ 1 ]^2 + xdim[2]^2 ) )
    out$lossfun <- lossfun
    out$lossfun.args <- list(...)

    if( !is.null( threshold ) ) {

	numthresh <- length( threshold )
	if( numthresh == 1 ) threshold <- rep( threshold, 3 )
	else if( numthresh > 3 ) stop("lossdiff: invalid threshold argument.  Must be of length one, two or three.")
	else if( numthresh == 2 ) threshold <- c( threshold[ 1 ], rep( threshold[ 2], 2 ) )

	x[ x < threshold[ 1 ] ] <- 0
	xhat1[ xhat1 < threshold[ 2 ] ] <- 0
	xhat2[ xhat2 < threshold[ 3 ] ] <- 0

	out$threshold <- threshold

    } # end of if 'threshold' stmt.
 
    g1 <- do.call(lossfun, args=c(list(x=x,y=xhat1),list(...)))
    g2 <- do.call(lossfun, args=c(list(x=x,y=xhat2),list(...)))
    d <- matrix(g1-g2, xdim[1], xdim[2])

    if(zero.out) {

	# zeros <- d == 0
	zeros <- (x==0) & (xhat1==0) & (xhat2==0)
	beta <- mean(!zeros, na.rm=TRUE)
	out$zeros <- zeros
	out$beta <- beta

    }

#     if(trend=="ols") {
 
 	if(is.null(loc)) loc <- cbind(rep(1:xdim[1],xdim[2]), rep(1:xdim[2],each=xdim[1]))
 	# else out$loc <- character(substitute(loc))
 	if( zero.out ) dat <- data.frame(y = c(d)[ !zeros ], x1 = loc[ !zeros, 1], x2 = loc[ !zeros, 2])
	else dat <- data.frame(y = c(d), x1 = loc[, 1], x2 = loc[, 2] )
 	fit <- try( lm(y~x1+x2, data=dat), silent = TRUE )
 	# tr <- matrix(predict(fit),xdim[1],xdim[2])
 	# if(zero.out) tr[zeros] <- 0
 	# d <- d - tr
 	out$trend.fit <- fit
# 
#     } else if(is.numeric(trend)) {
#
#	if(zeros) warning("lossdiff: zero.out is TRUE, but trend provided.  Not setting original zeros back after removing trend.")
# 	d <- d - trend
# 	if( zero.out ) d[ zero.out ] <- 0
# 	out$trend <- trend
#
#   }

    out$call <- theCall
    out$d <- d
    class(out) <- "lossdiff"

    return(out)

} # end of 'lossdiff.default' function.

empiricalVG.lossdiff <- function( x, trend = 0, maxrad, dx = 1, dy = 1 ) {

    theCall <- match.call()

    a <- attributes( x )

    out <- x

    if( missing( maxrad ) ) maxrad <- floor( max( a$xdim ) / 2 )

    if( is.null( maxrad ) || is.na( maxrad ) || !is.finite( maxrad ) || !is.numeric( maxrad ) ) maxrad <- 20

    out$trend <- trend
    out$vgram.args <- list( maxrad = maxrad, dx = dx, dy = dy )

    if( !is.null( x$zeros ) ) zero.out <- TRUE
    else zero.out <- FALSE

    d <- x$d - trend

    if( zero.out ) d[ x$zeros ] <- 0

    out$loss.differential.detrended <- d

    if( !zero.out ) vg <- vgram.matrix(dat = d, R = maxrad, dx = dx, dy = dy)
    else vg <- variogram.matrix(dat = d, R = maxrad, dx = dx, dy = dy, zero.out = zero.out)

    out$lossdiff.vgram <- vg

    out$call <- list( out$call, theCall )

    return( out )

} # end of 'empiricalVG.lossdiff' function.

flossdiff <- function(object, vgmodel = "expvg", ... ) {

    out <- object
    a <- attributes( object )
    vg <- object$lossdiff.vgram

    p <- c( sqrt( vg$vgram[ 1 ] ), ifelse( vg$vgram[ 2 ] > 0, -(vg$vgram[ 2 ] - vg$vgram[ 1 ]), -(0 - vg$vgram[ 1 ]) ) )

    fit <- try( nlminb( p, ORSS, vg = vg, model = vgmodel, ..., lower = c(0, 0), upper = c(Inf, Inf) ) )

    hold <- list( model = vgmodel, p = fit$par, objective = fit$objective, maxdim = object$maxdim,
		convergence = fit$convergence, message = fit$message, iterations = fit$iterations,
		evaluations = fit$evaluations )

    class( hold ) <- paste("flossdiff", vgmodel, sep = ".")

    out$vgmodel <- hold

    # class( out ) <- c("lossdiff", "flossdiff")

    return(out)

} # end of 'flossdiff' function.

expvg <- function( p, vg, ... ) {

    return( p[ 1 ] * ( 1 - exp( - vg$d * p[ 2 ] ) ) )

} # end of 'expvg' function.

predict.flossdiff.expvg <- function( object, newdata, ... ) {

    if( missing( newdata ) ) {

	h <- seq(0, object$maxdim - 1, by = 1 )

    } else h <- newdata

    p <- object$p

    return( p[ 1 ] * ( 1 - exp( - h * p[ 2 ] ) ) )

} # end of 'predict.expvg' function.

print.flossdiff.expvg <- function( x, ... ) {

    cat("\n\nExponential variogram fit\n")

    look <- x$p
    names( look ) <- c("scale", "range")

    print( look )

    invisible()

} # end of 'print.flossdiff.expvg' function.

ORSS <- function( p, vg, model = "expvg", ... ) {

    fit <- do.call( model, c( list( p = p, vg = vg ), list( ... ) ) )

    out <- sum( ( fit - vg$vgram )^2, na.rm = TRUE )

    return( out )

} # end of 'ORSS' function.

# TO DO: make a WRSS function too for weighted RSS (maybe not necessary).

plot.lossdiff <- function(x, ..., icol=c("gray", tim.colors(64))) {

    tmp <- attributes(x)
    loc.byrow <- tmp$loc.byrow

    if(is.null(tmp$msg)) msg <- paste("\n", x$lossfun, ": ", x$data.name[2], " vs ", x$data.name[3], " (", x$data.name[1], ")", sep="")
    else msg <- paste(tmp$msg, "\n", x$lossfun, ": ", x$data.name[2], " vs ", x$data.name[3], " (", x$data.name[1], ")", sep="")
    
    par(mfrow=c(2,2), mar=rep(4.1,4), oma=c(0,0,2,0))

    if(is.null(tmp$projection)) proj <- FALSE
    else proj <- tmp$projection

    if(is.null(tmp$map)) domap <- FALSE
    else domap <- tmp$map

    if(is.null(tmp$xdim)) proj <- domap <- FALSE
    else xd <- tmp$xdim

    if(proj) {
	loc <- list(x=matrix(tmp$loc[,1], xd[1], xd[2], byrow=loc.byrow),
		    y=matrix(tmp$loc[,2], xd[1], xd[2], byrow=loc.byrow))
    }

    if(domap) {
	locr <- apply(tmp$loc, 2, range, finite=TRUE)
	ax <- list(x=pretty(round(tmp$loc[,1], digits=2)), y=pretty(round(tmp$loc[,2], digits=2)))
    }

    if(is.null(x$zeros)) Im <- x$d
    else {
        Im <- x$d
        Im[x$zeros] <- NA
    }

    if(domap) {
	map(xlim=locr[,1], ylim=locr[,2], type="n")
	# axis(1, at=ax$x, labels=ax$x)
	# axis(2, at=ax$y, labels=ax$y)
        if(proj) image.plot(loc$x, loc$y, Im, col=icol, add=TRUE, ...)
        else image.plot(Im, col=icol, add=TRUE, ...)
        map(add=TRUE, lwd=1.5)
        map(add=TRUE, database="state")
	map.axes()
    } else {
	if(proj) image.plot(loc$x, loc$y, Im, col=icol, ...)
        else image.plot(Im, col=icol, ...)
    }
    title("")
    title("Loss Differential Field")

    hist(x$d, breaks="FD", xlab="Mean Loss Differential", col="darkblue", freq=FALSE,
	main="Histogram of\nLoss Differential Field")
    a <- x$lossdiff.vgram
    plot(a$d, a$vgram, col="darkblue", xlab="separation distance", ylab="variogram")
    if(!is.null(x$vgmodel)) {
	lines(a$d, predict(x$vgmodel, newdata = a$d), col="darkorange", lwd=1.5)
   	legend("bottomright", legend=c("Empirical", "Model"), pch=c("o", ""), col=c("darkblue","darkorange"), lty=c(0,1), lwd=1.5, bty="n")
    }
    plot.vgram.matrix(a, main="variogram by direction")

    title("")
    mtext(msg, line=0.05, outer=TRUE)

    invisible()

} # end of 'plot.lossdiff' function.

print.lossdiff <- function( x, ... ) {

    print( x$call )

    print( x$data.name )

    if( !is.null( x$trend ) ) print( summary( x$trend ) )

    print( x$lossfun )

    if( !is.null( x$vgmodel ) ) {

	cat("\n", "Fitted variogram model:\n")

	print( x$vgmodel )

    } # end of if model has been fit function.

    invisible()

} # end of 'print.lossdiff' function.

summary.lossdiff <- function(object, ...) {

   out <- object
   out$summary.call <- match.call()
   msg <- paste(object$lossfun, ": ", object$data.name[2], " vs ", object$data.name[3], " (against verification: ", object$data.name[1], ")", sep="")
   print(msg)
   cat("Fitted trend information (note: not used by subsequent functions, information only):\n\n")
   print( summary( object$trend.fit ) )
   d <- object$d
   if(zero.out <- !is.null(object$beta)) {
	cat("\n", "Estimate of beta present, so calculating Dbar and test statistic over non-zero entries only.\n")
	cat("\n", "Frequency of non-zero loss differential (beta) is: ", object$beta, "\n")
	d[d==0] <- NA
   }
   good <- !is.na(d)
   n <- length(d[good])
   cat("\n", "number of non-zero loss differential points is: ", n, "\n")
   cat("\n", "Mean Loss Differential: \n")
   Dbar <- sum(colSums(d,na.rm=TRUE),na.rm=TRUE)/n
   out$Dbar <- Dbar
   cat("\n", Dbar, "\n", "\n")
   a <- object$lossdiff.vgram
   cat("\n", "Summary of empirical variogram values:\n")
   print(stats(a))
   if(!is.null(object$vgmodel)) {
      # co <- coef(object$vgmodel)
      # sig2 <- co[1]^2
      # r <- co[2]
      # b <- sig2*exp(-a$d.full/r)
      # if(zero.out) denom <- sqrt(mean(object$beta^2*b,na.rm=TRUE))
      # else 
      b <- predict(object$vgmodel)
      denom <- sqrt(mean(b,na.rm=TRUE))
      SV <- Dbar/denom
      out$test.statistic <- SV
      cat("Test Statistic for null hypothesis of equal predictive ability on average\n")
      print(SV)
      pval <- c( 2 * pnorm( -abs(SV) ), pnorm( SV ), pnorm( SV, lower.tail=FALSE ) )
      names(pval) <- c("two.sided", "less", "greater")
      cat("p-value for two-sided alternative hypothesis is: ", pval[1], "\n")
      cat("p-value for (one-sided) alternative hypothesis that mu(D) < 0 is: ", pval[2], "\n")
      cat("p-value for (one-sided) alternative hypothesis that mu(D) > 0 is: ", pval[3], "\n")
      out$p.value <- pval
   }

   invisible(out)

} # end of 'summary.lossdiff' function.

variogram.matrix <- function (dat, R = 5, dx = 1, dy = 1, zero.out = FALSE) 
{
    SI <- function(ntemp, delta) {
        n1 <- 1:ntemp
        n2 <- n1 + delta
        good <- (n2 >= 1) & (n2 <= ntemp)
        cbind(n1[good], n2[good])
    }
    if(zero.out) dat[dat==0] <- NA
    N <- ncol(dat)
    M <- nrow(dat)
    m <- min(c(round(R/dx), M))
    n <- min(c(round(R/dy), N))
    ind <- rbind(as.matrix(expand.grid(0, 1:n)), as.matrix(expand.grid(1:m, 
        0)), as.matrix(expand.grid(c(-(m:1), 1:m), 1:n)))
    d <- sqrt((dx * ind[, 1])^2 + (dy * ind[, 2])^2)
    good <- (d > 0) & (d <= R)
    ind <- ind[good, ]
    d <- d[good]
    ind <- ind[order(d), ]
    d <- sort(d)
    nbin <- nrow(ind)
    holdVG <- rep(NA, nbin)
    holdN <- rep(NA, nbin)
    for (k in 1:nbin) {
        MM <- SI(M, ind[k, 1])
        NN <- SI(N, ind[k, 2])
	numNA <- sum(is.na(dat[MM[,1],NN[,1]]) | is.na(dat[MM[,2],NN[,2]]),na.rm=TRUE)
        holdN[k] <- length(MM) * length(NN) - numNA
        BigDiff <- (dat[MM[, 1], NN[, 1]] - dat[MM[, 2], NN[,2]])
        holdVG[k] <- mean(0.5 * (BigDiff)^2, na.rm=TRUE)
    }
    top <- tapply(holdVG * holdN, d, FUN = "sum")
    bottom <- tapply(holdN, d, FUN = "sum")
    dcollapsed <- as.numeric(names(bottom))
    vgram <- top/bottom
    dimnames(vgram) <- NULL
    out <- list(vgram = vgram, d = dcollapsed, ind = ind, d.full = d, 
        vgram.full = holdVG, N = holdN, dx = dx, dy = dy)
    class(out) <- "vgram.matrix"
    return(out)
} # end of 'variogram.matrix' function.

structurogram.matrix <- function(dat, q=2, R=5, dx=1, dy=1, zero.out=FALSE) {
   SI <- function(ntemp, delta) {
        n1 <- 1:ntemp
        n2 <- n1 + delta
        good <- (n2 >= 1) & (n2 <= ntemp)
        cbind(n1[good], n2[good])
    }
    if(zero.out) dat[dat==0] <- NA
    N <- ncol(dat)
    M <- nrow(dat)
    m <- min(c(round(R/dx), M))
    n <- min(c(round(R/dy), N))
    ind <- rbind(as.matrix(expand.grid(0, 1:n)), as.matrix(expand.grid(1:m,
        0)), as.matrix(expand.grid(c(-(m:1), 1:m), 1:n)))
    d <- sqrt(abs(dx * ind[, 1])^2 + abs(dy * ind[, 2])^2)
    good <- (d > 0) & (d <= R)
    ind <- ind[good, ]
    d <- d[good]
    ind <- ind[order(d), ]
    d <- sort(d)
    nbin <- nrow(ind)
    holdVG <- rep(NA, nbin)
    holdN <- rep(NA, nbin)
    for (k in 1:nbin) {
        MM <- SI(M, ind[k, 1])
        NN <- SI(N, ind[k, 2])
        numNA <- sum(is.na(dat[MM[,1],NN[,1]]) | is.na(dat[MM[,2],NN[,2]]),na.rm=TRUE)
        holdN[k] <- length(MM) * length(NN) - numNA
        BigDiff <- (dat[MM[, 1], NN[, 1]] - dat[MM[, 2], NN[,2]])
        holdVG[k] <- mean(0.5 * (BigDiff)^q, na.rm=TRUE)
    }
    top <- tapply(holdVG * holdN, d, FUN = "sum")
    bottom <- tapply(holdN, d, FUN = "sum")
    dcollapsed <- as.numeric(names(bottom))
    vgram <- top/bottom
    dimnames(vgram) <- NULL
    out <- list(vgram = vgram, d = dcollapsed, ind = ind, d.full = d,
        vgram.full = holdVG, N = holdN, dx = dx, dy = dy, q=q)
    class(out) <- "structurogram.matrix"
    return(out)
} # end of 'structurogram.matrix' function.

plot.structurogram.matrix <- function(x,...) {
   par(mfrow=c(1,2), bg="beige")
   plot( x$d, x$vgram, xlab="separation distance", ylab=paste("structure (q=", x$q, ")", sep=""), ...)
   points( x$d.full, x$vgram.full, pch=".")
   plot.vgram.matrix(x, main="Structure by direction")
   invisible()
} # end of 'plot.structurogram.matrix' function.

structurogram <- function(loc, y, q=2, id = NULL, d = NULL, lon.lat = FALSE, dmax = NULL, N = NULL, breaks = NULL) 
{
    y <- cbind(y)
    if (is.null(id)) {
        n <- nrow(loc)
        ind <- rep(1:n, n) > rep(1:n, rep(n, n))
        id <- cbind(rep(1:n, n), rep(1:n, rep(n, n)))[ind, ]
    }
    if (is.null(d)) {
        loc <- as.matrix(loc)
        if (lon.lat) {
            d <- rdist.earth(loc)[id]
        }
        else {
            d <- rdist(loc, loc)[id]
        }
    }
    vg <- 0.5 * rowMeans(cbind(abs(y[id[, 1], ] - y[id[, 2], ])^q), 
        na.rm = TRUE)
    call <- match.call()
    if (is.null(dmax)) {
        dmax <- max(d)
    }
    od <- order(d)
    d <- d[od]
    vg <- vg[od]
    ind <- d <= dmax & !is.na(vg)
    out <- list(d = d[ind], val = vg[ind], call = call, q=q)
    if (!is.null(breaks) | !is.null(N)) {
        out <- c(out, stats.bin(d[ind], vg[ind], N = N, breaks = breaks))
    }
    class(out) <- "structurogram"
    return(out)
} # end of 'structurogram' function.

plot.structurogram <- function(x,...) {
   plot(x$d, x$val, xlab="separation distance", ylab=paste("structure (q=",x$q,")", sep=""), ...)
   lines(x$centers, x$stats["mean",], col="darkblue")
   invisible()
} # end of 'plot.structurogram' function.

variographier <- function( x, init, zero.out = FALSE, ... ) {

    UseMethod( "variographier", x )

} # end of 'variographier' function.

variographier.default <- function( x, init, zero.out = FALSE, ..., y ) {

    if( !zero.out ) {

	vgx <- try( vgram.matrix( dat = x, ... ) )
	vgy <- try( vgram.matrix( dat = y, ... ) )

    } else {

	vgx <- try( variogram.matrix( dat = x, zero.out = zero.out, ... ) )
	vgy <- try( variogram.matrix( dat = y, zero.out = zero.out, ... ) )

    }

    if( class( vgx ) == "try-error" || class( vgy ) == "try-error" ) return( NA )

    if( missing( init ) ) {

        px <- c( sqrt( vgx$vgram[ 1 ] ), ifelse( vgx$vgram[ 2 ] > 0, -(vgx$vgram[ 2 ] - vgx$vgram[ 1 ]), -(0 - vgx$vgram[ 1 ]) ) )
	py <- c( sqrt( vgy$vgram[ 1 ] ), ifelse( vgy$vgram[ 2 ] > 0, -(vgy$vgram[ 2 ] - vgy$vgram[ 1 ]), -(0 - vgy$vgram[ 1 ]) ) )

    } else {

	px <- init$px
	py <- init$py

    }

    fitx <- try( nlminb( px, ORSS, vg = vgx, model = "expvg", ..., lower = c(0, 0), upper = c(Inf, Inf) ) )
    if( class( fitx ) == "try-error" ) return( NA )
    fity <- try( nlminb( py, ORSS, vg = vgy, model = "expvg", ..., lower = c(0, 0), upper = c(Inf, Inf) ) )
    if( class( fity ) == "try-error" ) return( NA )

    res <- 1 / sqrt( fitx$par[ 1 ]^2 + fity$par[ 1 ]^2 + ( 3 / fitx$par[ 2 ] - 3 / fity$par[ 2 ] )^2 )

    out <- list( obs.vg = vgx, mod.vg = vgy, obs.parvg = fitx, mod.parvg = fity, variography = res )

    class( out ) <- "variographied"

    return( out )

} # end of 'variographier.default' function.

variographier.SpatialVx <- function( x, init, zero.out = FALSE, ..., obs = 1, model = 1, time.point = 1 ) {

    a <- attributes( x )

    dat <- datagrabber( x, time.point = time.point, obs = obs, model = model)
    Obs <- dat$X
    Fcst <- dat$Xhat
    mainname <- a$data.name
    vxname <- a$obs.name[obs]

    if( !missing( init ) ) out <- variographier( x = Obs, y = Fcst, init = init, zero.out = zero.out, ... )
    else out <- variographier( x = Obs, y = Fcst, zero.out = zero.out, ... )

    attr( out, "time.point") <- time.point
    attr( out, "model") <- model
    attr( out, "model.name" ) <- a$model.name[ model ]
    attr( out, "obs.name" ) <- a$obs.name[ obs ]
    attr( out, "data.name" ) <- c( mainname, vxname, a$model.name[ model ] )

    return( out )

} # end of 'variographier.SpatialVx' function.

print.variographied <- function( x, ... ) {

    a <- attributes( x )

    if( !is.null( x$data.name ) ) print( a$data.name )

    cat( "\nEstimated Variography Value = ", x$variography, "\n\n" )

    invisible()

} # end of 'print.variographied' function.

plot.variographied <- function( x, ... ) {

    yl <- range( c( x$obs.vg$vgram, x$mod.vg$vgram ), finite = TRUE )

    vg <- x$obs.vg
    plot( vg$d, vg$vgram, col = "darkblue", xlab = "distance (grid squares)", ylab = "variogram", ylim = yl, ... )
    vg <- x$mod.vg
    points( vg$d, vg$vgram, col = "darkred", ... )

    p <- x$obs.parvg$par
    modvg <- p[1] * ( 1 - exp( -vg$d * p[2] ) )
    lines( vg$d, modvg, col = "darkblue" )
    p <- x$mod.parvg$par
    modvg <- p[1] * (1 - exp(-vg$d * p[2] ) )
    lines( vg$d, modvg, col = "darkred", lty = 2 )

    a <- attributes( x )
    if( !is.null( a$obs.name ) ) xname <- a$obs.name
    else xname <- "observation"
    if( !is.null( a$mod.name ) ) yname <- a$mod.name
    else yname <- "model"

    zname <- paste( "variography = ", x$variography )

    legend( "topleft", legend = c( xname, yname, zname ), pch = c( 1, 1, 0 ), col = c( "darkblue", "darkred", "white" ), bty = "n" )

    invisible()

} # end of 'plot.variographied' function

Try the SpatialVx package in your browser

Any scripts or data that you put into this service are public.

SpatialVx documentation built on March 28, 2021, 1:10 a.m.