scattercontour <- function(x,y,xlim=NULL,ylim=NULL,levels=NULL, xlab="",ylab="") { x.nbins <- 50 y.nbins <- 50 x.max <- ifelse(is.null(xlim) || length(xlim) != 2, max(x, na.rm=TRUE), xlim[1]) x.min <- ifelse(is.null(xlim) || length(xlim) != 2, min(x, na.rm=TRUE), xlim[2]) x.binwidth <- (x.max - x.min)/x.nbins x.binlocs <- seq(x.min+x.binwidth/2, x.max, by=x.binwidth) x.ok <- x >= x.min & x <= x.max x.bins <- factor(ifelse(x.ok, pmin(floor((x-x.min)/x.binwidth)+1,x.nbins, na.rm=TRUE), NA), levels=1:50) y.max <- ifelse(is.null(ylim) || length(ylim) != 2, max(y, na.rm=TRUE), ylim[1]) y.min <- ifelse(is.null(ylim) || length(ylim) != 2, min(y, na.rm=TRUE), ylim[2]) y.binwidth <- (y.max - y.min)/y.nbins y.binlocs <- seq(y.min+y.binwidth/2, y.max, by=y.binwidth) y.ok <- y >= y.min & y <= y.max y.bins <- factor(ifelse(y.ok, pmin(floor((y-y.min)/y.binwidth)+1,y.nbins, na.rm=TRUE), NA), levels=1:50) grid.counts <- xtabs(~ x.bins + y.bins) contour(x.binlocs, y.binlocs, grid.counts, levels=levels, xlab=xlab, ylab=ylab) for (i in 1:x.nbins) { for (j in 1:y.nbins) { if (grid.counts[i,j] > 0 & grid.counts[i,j] < 30) { points(x[x.bins == i & y.bins == j], y[x.bins == i & y.bins == j], pch=".") } } } }