xyz.plot <- function (data, z, x, y, residual.plot = FALSE, lower.limit, upper.limit, origin, maxsize = 0.075, seg.col = "red", seg.lty = par("lty"), seg.lwd = par("lwd"), seg.azm = 0, pch.pnt = 16, plot.points = FALSE, equal.scale = TRUE, messages.screen = FALSE, ... ) { # purpose: function plots the x,y coordinates of a point and # draws line segments through the points whereby the # length of the segments is linearly related to z # # # arguments: # # data dataframe with x, y, and z-components or # list of class "geodata" from which missing # x, y, z are taken # z value of target variable # x x-coordinate of sampling location # y y-coordinate of sampling location # residual.plot TRUE: residual type plot, # sign(z-lower.limit) 45 degrees slanted segments, # FALSE: normal type plot, vertical or horizontal segments # lower.limit value for lower.limit of target variable, # default value: lower.limit=mean(z), # there are three possibilities for drawing segments # 1) lower.limit < min(z): normal type plot, # segment length proportional to (z - lower.limit) # 2) lower.limit in range(z) and residual.plot = TRUE: residual type plot; # segment length proportional to abs(z - lower.limit) # 3) else: normal type plot, # segment length proportional to (z - min(z)) # upper.limit value for upper.limit of target varible # default value: upper.limit=max(z) # upper.limit > max(z) then the segment length is given by # (z - lower.limit) / (upper.limit - lower.limit) # origin the same as lower.limit (compatibility) # maxsize maximum length of segment in y-units # equal.scale TRUE: calls eqscplot() of library MASS # FALSE: calls basic plot() # seg.col color attribute for segments # seg.lty line type attribute for segments # seg.lwd line width attribute for segments # seg.azm azimuth angle of segments in degrees (angle clockwise from north) # pch.pnt plotting character for plot() # plot.points TRUE: sampling locations are shown # FALSE: sampling locations are not shown # messages.screen write diagnostics to stdout # ... further arguments passed to plot(), eqscplot(), points() # # author: A. Papritz # date: 2003-04-01 # setup x, y, z vectors if ( !missing(data) ) { if( missing(x) ) { if( is.data.frame(data) ) { if( any( names(data) == 'x' ) ) x <- data$x else stop('dataframe without "x"-component ') } else { if( is.list(data) && class(data) == "geodata") x <- data$coords[,1] else stop('object given for argument "data" not supported') } } if( missing(y) ) { if( is.data.frame(data) ){ if( any( names(data) == 'y' ) ) y <- data$y else stop('dataframe without "y"-component ') } else { if( is.list(data) && class(data) == "geodata") y <- data$coords[,2] else stop('object given for argument "data" not supported') } } if( missing(z) ) { if( is.data.frame(data) ){ if( any( names(data) == 'z' ) ) z <- data$z else stop('dataframe without "z"-component ') } else { if( is.list(data) && class(data) == "geodata") z <- data$data else stop('object given for argument "data" not supported') } } } else { if ( missing(x) ) stop('"x" is a required argument') if ( missing(y) ) stop('"y" is a required argument') if ( missing(z) ) stop('"z" is a required argument') } # check mode of important arguments if( !is.numeric( x ) ) stop( "x-coordinate must be of mode numeric" ) if( !is.numeric( y ) ) stop( "y-coordinate must be of mode numeric" ) if( !is.numeric( z ) ) stop( "z-value must be of mode numeric" ) if( !( is.character( seg.col ) | is.numeric( seg.col ) ) ) stop("seq.col must be either of mode character or numeric") if( !is.numeric( seg.azm )) stop( "seg.azm must be of mode numeric" ) # omit NAs in x, y or z n <- min( length(x), length(y), length(z) ) if( length(seg.col) == 1 ) (seg.col <- rep(seg.col, n) ) if( length(seg.azm) == 1 ) (seg.azm <- rep(seg.azm, n) ) df <- data.frame(x = x[1:n], y = y[1:n], z = z[1:n], seg.col = I(seg.col[1:n]), seg.azm = seg.azm[1:n] ) df <- na.omit(df) x <- df$x; y <- df$y; z <- df$z; seg.col <- df$seg.col; seg.azm <- df$seg.azm if( messages.screen ) cat( ' number of omitted cases with NAs: ', n-length(x), '\n' ) rm(df) # set default value for lower.limit if ( missing(lower.limit) ) { if (missing(origin)) lower.limit <- mean(z) else lower.limit <- origin } # plot sample locations if( plot.points ) pt <- "p" else pt <- "n" if( equal.scale ) t.asp <- 1 else t.asp <- NA plot(x, y, pch=pch.pnt, type = pt, asp = t.asp, ... ) # compute length of line segments wdim <- par('usr'); wx <- diff(wdim[1:2]); wy <- diff(wdim[3:4]) if( messages.screen ) cat(' extremes of user coordinates in plotting window ', format(wdim), '\n') maxsize <- wy * maxsize * 0.5 rz <- range(z) if( messages.screen ) cat( ' range of z: ', rz, '\n' ) if ( missing(upper.limit) ) upper.limit <- rz[2] else { if (upper.limit <= rz[2]) upper.limit <- rz[2] } if ( lower.limit < rz[1] ) { dy <- (z - lower.limit) / (upper.limit - lower.limit) if( messages.screen ) { cat( ' value used for lower.limit: ', lower.limit, '\n' ) cat( ' length of segments proportional to ( z - lower.limit )\n' ) } } else { if ( lower.limit < rz[2] && residual.plot ) { dy <- (z - lower.limit) / max( abs( c( rz[1]-lower.limit, rz[2]-lower.limit ) ) ) if( messages.screen ) { cat( ' value used for lower.limit: ', lower.limit, '\n' ) cat( ' length of segments proportional to ( z - lower.limit )\n' ) } } else { dy <- (z-rz[1]) / (upper.limit - rz[1]) if( messages.screen ) cat( ' length of segments proportional to ( z - min(z) )\n' ) } } dy <- dy * maxsize # draw segments oldpar <- par(xpd=TRUE) pwin <- par('pin') # print(dy) # print(seg.col) # print(seg.azm) if ( all( dy >= 0 ) ) { dx <- sin(seg.azm/180*pi) * dy dy <- cos(seg.azm/180*pi) * dy segments( x-dx, y-dy, x+dx, y+dy, col=seg.col, lty=seg.lty, lwd=seg.lwd ) } else segments( x-abs(dy/wy*pwin[2]*wx/pwin[1]), y-dy, x+abs(dy/wy*pwin[2]*wx/pwin[1]), y+dy, col=seg.col, lty=seg.lty, lwd=seg.lwd ) par(oldpar) } azimuth.angle <- function(x, y) { # purpose: function computes azimuth angle after centering # on (center.x, center.y) # # arguments: # # x x-coordinates # y y-coordinates # # author: A. Papritz # date: 2001-06-05 x <- x[!is.na(x)] y <- y[!is.na(y)] if(length(x) != length(y)) stop('lengths of "x" and "y" differ') angle <- -Arg( complex( real=x, imag=y )) angle[angle<0] <- 2*pi + angle[angle<0] angle <- angle + pi/2 angle[angle>2*pi] <- angle[angle>2*pi] - 2*pi angle }