R community:

I have been creating code for plotting nomographs, or multiple, overlain contour plots of z-variables on a common x- and y- variable.  My input has been a matrix with observed x, y, and multiple z variables; I then create a trend surface using trmat for each z-variable.  So far so good.

One application I have for these, requires shading a portion of the nomogram that meets criteria for some of the z-variables (i.e., z[1] must be between 20 and 30, z[2] must be less than 40, etc.).  My solution was to use a logical comparison on each contour surface provided by trmat, sum the "logical surfaces" up and see if they were less than the total number of criteria.  It works, but it is quite inefficient even if I vectorize the code somewhat; for example if I specify a gridsize of 200 in trmat, have 5 z variables, and 1 criteria for each, I will have well over 200,000 comparisons to make!  So I am looking for hints or maybe an entirely different approach to speed this up.

I attached the crit.region function below along with my write up on how it works.  Can somebody give me some ideas on how to proceed?

Thanks,
Mike

Mike R. Saunders
Forest Biometrician
Cooperative Forest Research Unit
University of Maine
5755 Nutting Hall
Orono, ME  04469-5755

207-581-2763 (O)
207-581-2833 (F)


# The following function selects a region that meets a set of 
# criteria defined in terms of z-variables in a list from nomogram
# or a similarly formatted list.  This function basically is a set
# of logical comparisons on z-values at each xy-coordinate.  As such,
# the function is rasterized and can take considerable time when
# each z-variable matrix is quite large.  Parameters for the 
# function are:
# 
#   1) x        (Required)  Either a list consisting of a vector
#                           of gridded x-coordinates, a vector of 
#                           gridded y-coordinates and matrices of
#                           each z-variable, or a vector of just 
#                           the gridded x-coordinates.
#   2) y        (Optional)  A vector of gridded y-coordinates.
#   3) z        (Optional)  A matrix or data.frame of z-variates
#                           that correspond to the gridded 
#                           xy-coordinates.
#   4) critmat  (Required)  A matrix or data.frame with rows equal
#                           to the number of z-variables and 2 
#                           columns.  The first column corresponds
#                           to the minimum value allowed for each
#                           z-variable, the second to the maximum
#                           value.  If there is no minimum or
#                           maximum for a variable, NA should be 
#                           used in the appropriate row and column.
# 
# This function returns the critical area as a matrix of NA and 1 
# with dimension equal to a z-variable matrix.  The function also
# returns a message if there is no critical area solution.
# 
# [Future versions of this function will try to improve its
#  computational speed.]
# 
crit.region<-function(x,y=NULL,z=NULL,critmat) {
    if(all(missing(y),missing(z))) {
        stopifnot(class(x)=="list",sum(lapply(x,class)[1:2]!="numeric")==0,sum(sapply(x,class)[3:length(x)]!="matrix")==0,length(x[[1]])==dim(x[[3]])[1],length(x[[2]])==dim(x[[3]])[2],length(x)>4)
        y<-x[[2]]
        z<-x[c(3:length(x))]
        x<-x[[1]]
    } else if(any(missing(y),missing(z))) {
        stop("y and z are both required unless x is properly formatted list")
    } else stopifnot(class(y)=="numeric",class(z)=="list",length(x)==dim(z[[1]])[1],length(y)==dim(z[[1]])[2],sum(sapply(z,class)!="matrix")==0)
    w<-length(z)
    zrange<-sapply(z,range,na.rm=T)
    stopifnot(class(critmat)%in%c("matrix","data.frame"),dim(critmat)==c(w,2))
    critarea<-matrix(data=0,nrow=dim(z[[1]])[1],ncol=dim(z[[1]])[2])
    for(i in 1:w) {
        minz<-ifelse(is.na(critmat[i,1]),zrange[1,i],critmat[i,1])
        maxz<-ifelse(is.na(critmat[i,2]),zrange[2,i],critmat[i,2])
        critarea<-critarea+apply(z[[i]],c(1,2), function(x) ifelse(x>minz & x<maxz,1,0))
        }     
    critarea<-apply(critarea,c(1,2), function(x) ifelse(x==w,1,NA))
    if(sum(critarea,na.rm=T)==0) message("Critical region is empty set!")
    return(critarea)
}



	[[alternative HTML version deleted]]


