R/DataObjects.R

Defines functions hist.SpatialVx thresholder.SpatialVx thresholder.default thresholder datagrabber.matched datagrabber.features datagrabber.SpatialVx plot.SpatialVxNoMap plot.SpatialVxMap plot.SpatialVx summary.SpatialVx print.SpatialVx make.SpatialVx

Documented in datagrabber.features datagrabber.matched datagrabber.SpatialVx hist.SpatialVx make.SpatialVx plot.SpatialVx plot.SpatialVxMap plot.SpatialVxNoMap print.SpatialVx summary.SpatialVx thresholder thresholder.default thresholder.SpatialVx

make.SpatialVx <- function(X, Xhat, thresholds=NULL, loc=NULL, projection=FALSE, subset=NULL, time.vals = NULL,
			    reg.grid = TRUE, map = FALSE, loc.byrow = FALSE, field.type="", units="", data.name = "",
			    obs.name = "X", model.name = "Xhat",
			    q=c(0, 0.1, 0.25, 0.33, 0.5, 0.66, 0.75, 0.9, 0.95), qs=NULL) {

    if( is.list( X ) ) {

	nobs <- length( X )
	xdim <- dim( X[[ 1 ]] )

    } else {

	nobs <- 1
	xdim <- dim( X )

    }

    if(is.list(Xhat)) {

	nforecast <- length(Xhat)
	ydim <- dim(Xhat[[1]])

    } else {

	nforecast <- 1
	ydim <- dim( Xhat )

    }

    if(all(xdim != ydim)) {

	stopmsg <- paste("make.SpatialVx: dim of X (", xdim[ 1 ], " by ", xdim[ 2 ],
		") must be the same as dim of (each component of) Xhat (", ydim[ 1 ], " by ",
		ydim[ 2 ], ")", sep = "" )

	stop( stopmsg )
	# stop("make.SpatialVx: dim of X must be the same as dim of (each component of) Xhat")

    } # end of stop bc dims are not equal stmt.

    out <- list(X, Xhat)

    if( is.null( thresholds ) ) {

	if( nobs > 1 ) {

	    othresh <- numeric(0)
	    for( i in 1:nobs ) othresh <- cbind( othresh, quantile( c( X[[ i ]] ), probs = q, na.rm = TRUE ) )

	} else othresh <- quantile( c( X ), probs = q, na.rm = TRUE )

	if( nforecast > 1 ) {

	    fthresh <- numeric(0)
	
	    for( i in 1:nforecast ) cbind( fthresh, quantile( c( Xhat[[ i ]] ), probs = q, na.rm = TRUE ) )

	} else fthresh <- quantile( c( Xhat ), probs = q, na.rm = TRUE )

	qs <- as.character(q)
        # qs <- paste( as.character(q), "quantile" )
	# rownames( othresh ) <- qs
	# rownames( fthresh ) <- qs

	thresholds <- list( X = othresh, Xhat = fthresh )

    } else if( is.list( thresholds ) ) {

	threshnames <- names( thresholds )
	if( !( "X" %in% threshnames ) || !( "Xhat" %in% threshnames ) ) {

	    stop( "make.SpatialVx: invalid thresholds argument.  List must have components named X and Xhat." )

	} 

	odim <- dim( thresholds$X )
 	fdim <- dim( thresholds$Xhat )

	if( odim[ 1 ] != fdim[ 1 ] ) stop( "make.SpatialVx: invalid thresholds argument.  X must have same number of thresholds (rows) as Xhat." )

	nth <- odim[ 1 ]

	if( is.null( odim ) ) thresholds$X <- matrix( thresholds$X, length( thresholds$X ), nobs )
	else if( odim[ 2 ] == 1 && nobs > 1 ) thresholds$X <- matrix( thresholds$X, odim[ 1 ], nobs )
	else if( odim[ 2 ] > 1 && odim[ 2 ] != nobs ) stop( "make.SpatialVx: invalid thresholds argument." )

	if( is.null( fdim ) ) thresholds$Xhat <- matrix( thresholds$Xhat, length( thresholds$Xhat ), nforecast )
        else if( fdim[ 2 ] == 1 && nforecast > 1 ) thresholds$X <- matrix( thresholds$Xhat, fdim[ 1 ], nforecast )
        else if( fdim[ 2 ] > 1 && fdim[ 2 ] != nobs ) stop( "make.SpatialVx: invalid thresholds argument." )

	qs <- paste("threshold", 1:nth )

    } else if( is.vector(thresholds) ) {

        qs <- as.character(thresholds)
        # thresholds <- cbind(thresholds, thresholds)
        nth <- length( thresholds )
        thresholds <- list( X = matrix( thresholds, nth, nobs ),
                            Xhat = matrix( thresholds, nth, nforecast ) )

    } else stop( "make.SpatialVx: invalid thresholds argument.  Must be a vector or list." )

    # udim <- dim(thresholds)
    # if( udim[2] == 2) colnames(thresholds) <- c("X", "Xhat")
    # else if( udim[2] > 2) colnames(thresholds) <- c("X", paste("Xhat", 1:(udim[2] - 1), sep=""))

    # rownames(thresholds) <- qs

    if( is.null( loc ) ) loc <- cbind( rep( 1:xdim[ 1 ], xdim[ 2 ] ), rep( 1:xdim[ 2 ], each = xdim[ 1 ] ) )

    if(is.null(time.vals)) {

	if(length(xdim) == 3) time.vals <- 1:xdim[3]
	else time.vals <- 1

    }

    if( !missing( data.name ) ) msg <- data.name
    else msg <- ""

    if(field.type != "" && units != "") msg <- paste(msg, "\n", field.type, " (", units, ")", sep="")
    else if(field.type != "") msg <- paste(msg, "\n", field.type, sep="")
    else if(units != "") msg <- paste(msg, "\n", " (", units, ")", sep="")


    class(out) <- "SpatialVx"
    attr(out, "xdim") <- xdim
    # attr(out, "udim") <- udim
    attr(out, "time") <- time.vals
    attr(out, "thresholds") <- thresholds
    attr(out, "loc") <- loc
    attr(out, "loc.byrow") <- loc.byrow
    attr(out, "subset") <- subset
    attr(out, "data.name") <- data.name
    attr(out, "obs.name" ) <- obs.name
    attr(out, "model.name" ) <- model.name
    attr(out, "nobs" ) <- nobs
    attr(out, "nforecast") <- nforecast
    attr(out, "field.type") <- field.type
    attr(out, "units") <- units
    attr(out, "projection") <- projection
    attr(out, "reg.grid") <- reg.grid
    attr(out, "map") <- map
    attr(out, "qs") <- qs
    attr(out, "msg") <- msg

    return(out)

} # end of 'make.SpatialVx' function.

print.SpatialVx <- function(x, ...) {
    tmp <- attributes(x)
    cat("SpatialVx data object: ", tmp$data.name, "\n")
    cat(tmp$msg, "\n")

    cat("\n\nObservations: ", tmp$obs.name, "\n\n" )
    cat("\n\nModels: ", tmp$model.name, "\n\n" )
    # if(tmp$field.type != "") cat("Field type: ", tmp$field.type)
    # if(tmp$units != "") cat(" (", tmp$units, ")", "\n")

    cat("Field dimensions: ", tmp$xdim, "\n")

    if(tmp$nforecast > 1) cat(tmp$nforecast, " forecasts to be evaluated/verified.\n")
    else cat("1 forecast to be evaluated/verified.\n")

    u <- tmp$thresholds
    cat("\nThresholds\n")
    print(u)

    if(!is.null(tmp$levels)) {
        cat("\n", "Neighborhood levels to be applied for neighborhood methods:\n")
        print(tmp$levels)
        cat("\n", "Smoothing function to be applied for neighborhood methods:\n")
        print(tmp$smooth.fun)
        if(!is.null(tmp$smooth.params)) {
            cat("\n", "Smoothing function parameters:\n")
            print(tmp$smooth.params)
        }
    }

    invisible()

} # end of 'print.SpatialVx' function.

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

    out <- list()

    x <- object
    print(x)
    tmp <- attributes(x)

    print( out$X <- summary(x[[1]]))

    if(tmp$nforecast > 1) print(out$Xhat <- lapply(x[[2]], summary))
    else print(out$Xhat <- summary(x[[2]]))

    invisible(out)
} # end of 'summary.SpatialVx' function.

plot.SpatialVx <- function( x, ..., time.point = 1, obs = 1, model = 1, col = c( "gray", tim.colors( 64 ) ), zlim, mfrow = c(1, 2) ) {

    a <- attributes( x )

    # if( a$map ) class( x ) <- "SpatialVxMap"
    # else class( x ) <- "SpatialVxNoMap"

    if( a$map ) {

	if( !missing( zlim ) ) plot.SpatialVxMap( x = x, ..., time.point = time.point, obs = obs, model = model, col = col, zlim = zlim, mfrow = mfrow )
	else plot.SpatialVxMap( x = x, ..., time.point = time.point, obs = obs, model = model, col = col, mfrow = mfrow )

    } else {

	if( !missing( zlim ) ) plot.SpatialVxNoMap( x = x, ..., time.point = time.point, obs = obs, model = model, col = col, zlim = zlim, mfrow = mfrow )
	else plot.SpatialVxNoMap(  x = x, ..., time.point = time.point, obs = obs, model = model, col = col, mfrow = mfrow )

    }

    # UseMethod( "plot", x )

}

plot.SpatialVxMap <- function( x, ..., time.point = 1, obs = 1, model = 1, col = c( "gray", tim.colors( 64 ) ), zlim, mfrow = c(1, 2) ) {

    if( !is.null( mfrow ) ) {

        op <- par()

        par( mfrow = mfrow, oma = c(0, 0, 2, 0) )

    }

    a <- attributes( x )

    loc <- a$loc
    xdim <- a$xdim

    l <- list( x = matrix( loc[, 1 ], xdim[ 1 ], xdim[ 2 ], byrow = a$loc.byrow ),
		y = matrix( loc[, 2 ], xdim[ 1 ], xdim[ 2 ], byrow = a$loc.byrow ) )

    lr <- apply( loc, 2, range, finite = TRUE )

    dat <- datagrabber( x, ..., time.point = time.point, obs = obs, model = model )
    X <- dat$X
    Xhat <- dat$Xhat

    if( missing( zlim ) ) zlim <- range( c( c(X), c(Xhat) ), finite = TRUE )

    if( is.function( time.point ) ) tp <- as.character( substitute( time.point ) )
    else tp <- time.point

    msg <- a$data.name
    msgX <- paste( a$obs.name[ obs ], "\n", "(", a$field.type, ", ", a$units, ", time = ", tp, ")", sep = "" )
    msgXhat <- paste( a$model.name[ model ], "\n", "(", a$field.type, ", ", a$units, ", time = ", tp, ")", sep = "" )

    map( xlim = lr[,1], ylim = lr[,2], type = "n" )
    title( msgX )
    poly.image( l$x, l$y, X, col = col, zlim = zlim, add = TRUE, ... )
    map( add = TRUE )
    map( database = "state", add = TRUE )
    map.axes()

    map( xlim = lr[,1], ylim = lr[,2], type = "n" )
    title( msgXhat )
    poly.image( l$x, l$y, Xhat, col = col, zlim = zlim, add = TRUE, ... )
    map( add = TRUE )
    map( database = "state", add = TRUE )
    map.axes()

    image.plot( X, zlim = zlim, col = col, legend.only = TRUE, ... )

    if( !is.null( mfrow ) ) {

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

        par( mfrow = op$mfrow )

    }

    invisible()

} # end of 'plot.SpatialVxMap' function.

plot.SpatialVxNoMap <- function( x, ..., time.point = 1, obs = 1, model = 1, col = c( "gray", tim.colors( 64 ) ), zlim, mfrow = c(1, 2) ) {

    if( !is.null( mfrow ) ) {

        op <- par()

        par( mfrow = mfrow, oma = c(0, 0, 2, 0) )

    }

    a <- attributes( x )

    dat <- datagrabber( x, ..., time.point = time.point, obs = obs, model = model )
    X <- dat$X
    Xhat <- dat$Xhat

    if( missing( zlim ) ) zlim <- range( c( c(X), c(Xhat) ), finite = TRUE )

    if( is.function( time.point ) ) tp <- as.character( substitute( time.point ) )
    else tp <- time.point

    msg <- a$data.name
    msgX <- paste( a$obs.name[ obs ], "\n", "(", a$field.type, ", ", a$units, ", time = ", tp, ")", sep = "" )
    msgXhat <- paste( a$model.name[ model ], "\n", "(", a$field.type, ", ", a$units, ", time = ", tp, ")", sep = "" )

    image( X, col = col, zlim = zlim, main = msgX, ... )
    image( Xhat, col = col, zlim = zlim, main = msgXhat, ... )
    image.plot( X, zlim = zlim, col = col, legend.only = TRUE, ... )

    if( !is.null( mfrow ) ) {

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

        par( mfrow = op$mfrow )

    }

    invisible()

} # end of 'plot.SpatialVxNoMap' function.

datagrabber.SpatialVx <- function(x, ..., time.point = 1, obs = 1, model = 1) {

    tmp <- attributes(x)

    if(!missing(time.point)) {

	if( !is.function( time.point ) ) {

            if(length(time.point) > 1) stop("datagrabber: length of time.point must be one.")
            tiid <- tmp$time
            Nt <- length(tiid)
            if(!is.numeric(time.point)) time.point <- (1:Nt)[ tiid == time.point ]
            if(is.na(time.point)) stop("datagrabber: invalid time.point argument.")

	}

    }

    if( !is.numeric( obs ) ) stop( "datagrabber: invalid obs argument.  Must be numeric." )
    if( !is.numeric( model ) ) stop( "datagrabber: invalid model argument.  Must be numeric." )

    xdim <- tmp$xdim
    nobs <- tmp$nobs
    nf <- tmp$nforecast

    Vx <- x[[1]]
    if( nobs > 1 ) Vx <- Vx[[ obs ]]
    
    Fcst <- x[[2]]
    if(nf > 1) Fcst <- Fcst[[ model ]]

    if(length(xdim) == 3) {

	if( is.function( time.point ) ) {

	    afun <- match.fun( time.point )

	    Vx <- apply( Vx, 1:2, afun, ... )
	    Fcst <- apply( Fcst, 1:2, afun, ... )

	} else {

	    Vx <- Vx[,,time.point]
	    Fcst <- Fcst[,,time.point]

	} # end of if take a single time point or a function over time.

    } # if more than one time point present stmts.

    out <- list( X = Vx, Xhat = Fcst )

    return( out ) 

} # end of 'datagrabber.SpatialVx' function.

datagrabber.features <- function(x, ...) {

    return(list(X = x$X, Xhat = x$Xhat))

} # end of 'data.grabber.features' function.

datagrabber.matched <- function(x, ...) {

    return(list(X = x$X, Xhat = x$Xhat))

} # end of 'data.grabber.matched' function.

thresholder <- function( x, type = c( "binary", "replace.below" ), th, rule = ">=", replace.with = 0, ... ) {

    if( !is.element( rule, c( ">=", ">", "<", "<=" ) ) ) stop("thresholder: invalid rule argument.")

    UseMethod( "thresholder" )

} # end of 'thresholder' function.

thresholder.default <- function( x, type = c( "binary", "replace.below" ), th, rule = ">=", replace.with = 0, ... ) {

    type <- match.arg( type )

    if( rule == ">=" ) id <- x >= th
    else if( rule == ">" ) id <- x > th
    else if( rule == "<=" ) id <- x <= th
    else if( rule == "<" ) id <- x < th

    xdim <- dim( x )

    if( is.null( xdim ) ) out <- numeric( length( x ) ) + replace.with
    else out <- matrix( replace.with, xdim[ 1 ], xdim[ 2 ] )

    if( type == "binary" ) out[ id ] <- 1 
    else out[ id ] <- x[ id ]

    return( out )

} # end of 'thresholder.default' function.

thresholder.SpatialVx <- function( x, type = c( "binary", "replace.below" ), th, rule = ">=", replace.with = 0, ...,
    time.point = 1, obs = 1, model = 1 ) {

    type <- match.arg( type )

    dat <- datagrabber( x, time.point = time.point, obs = obs, model = model )

    X <- dat$X
    Xhat <- dat$Xhat

    a <- attributes( x )

    u <- a$thresholds

    if( is.list( u ) ) {

	u1 <- u$X[ th, obs ]
	u2 <- u$Xhat[ th, model ]

    } else {

	u1 <- u2 <- u[ th ]

    }

    out1 <- thresholder( X, type = type, th = u1, rule = rule, replace.with = replace.with )
    out2 <- thresholder( Xhat, type = type, th = u2, rule = rule, replace.with = replace.with )

    return( list( X = out1, Xhat = out2 ) )

} # end of 'thresholder.SpatialVx' function.

hist.SpatialVx <- function( x, ..., time.point = 1, obs = 1, model = 1, threshold.num = NULL ) {

   tmp <- attributes( x )

    if( is.null( threshold.num ) ) {

	dat <- datagrabber( x, time.point = time.point, obs = obs, model = model )
	qs <- ""

    } else {

	dat <- thresholder( x, type = "replace.below", th = threshold.num, time.point = time.point, obs = obs, model = model )
	qs <- paste( ">= ", tmp$qs[ threshold.num ], sep = "" )

    }

   X <- dat$X
   Xhat <- dat$Xhat

   dn <- tmp$data.name
    vxname <- tmp$obs.name[ obs ]
    fcname <- tmp$model.name[ model ]
   nobs <- tmp$nobs
   nf <- tmp$nforecast

    if( !missing( obs ) ) if( length( obs ) > 1 ) stop("hist.SpatialVx: length of obs argument must be one." )

   if(!missing(model)) if(length(model) > 1) stop("hist.SpatialVx: length of model argument must be one.")

    h1 <- hist( X, main = dn, xlab = paste( vxname, "(", tmp$field.type, ", ", tmp$units, ")", sep = "" ), ... )
    h2 <- hist( Xhat, main = dn, xlab = paste( fcname, "(", tmp$field.type, ", ", tmp$units, ")", sep = "" ), ... )

   invisible( list( X = h1, Xhat = h2 ) )

} # end of 'hist.SpatialVx' 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.