"Cut out" 3D Surface Plot in R

Question

I've created a 3d surface plot in R using persp() - though I'm not wed to that package by any means if something else would solve this problem better. My variables on the x and y axis, respectively, are temperature and precipitation. I've created a grid that extends from my min-max temp, and min-max precip, and predicted values from my function for the z-axis.

However, realistically, my data does not truly "extend that far" (i.e. I don't have data at the tails [e.g. a combination of the lowest precip and the lowest temp]) so the extrapolation is a little misleading.

I'm wondering if there's an easy way to "bound" my plot, and I'm open to many types of suggestions. I could color-code values that are out of the range of observed weather conditions, or if there's a way to "cut out" the surface, that would be great too. For example, I've found the "bounds" (e.g. the highest and lowest observed precip for each rounded temp value). I could use the bounds as a "jagged cookie cutter," but I have also fit a Gaussian-like curve using scatter.smooth() from which I think I can get a smoothed function to describe the boundary.

Please let me know if you have any suggestions; I'm finding a solution hard to come by but maybe I'm just looking in the wrong place. Thanks!

Answer 1

If you want to do your "cookie-cutter" approach, you could do something like this with the sp library and point.in.polygon function. You would need the points specifying the perimeter you were interested in to contruct the polygon.

## Make some data
f <- function(theta) (1+1.5*cos(9*theta)) * (1+sin(theta)) + (1+0.5*cos(100*theta)) + (1+0.3*cos(50*theta))
x <- seq(-pi, pi, len=100)
dat <- data.frame(x=cos(x), y=sin(x)) * f(x)  # points for edges of polygon

library(sp)
d <- function(a, b) {  # density function to get z values
    ifelse(point.in.polygon(a, b, dat$x, dat$y) > 0, 1, 0)  # only points in range
}
ps <- outer(seq(min(dat$x), max(dat$x), len=100), 
             seq(min(dat$y), max(dat$y), len=100), d)
persp(ps, zlim=c(0, 3), theta=35, xlab="x", ylab="y", phi=25)

Answer 2

Take a look at the exclude.too.far() function in the mgcv package, which was designed for the exact same plotting purpose you describe but in relation to 2-d splines.

An illustration of what this function does is given below, taken from the example of ?exclude.too.far

library("mgcv")
x <- rnorm(100)
y <- rnorm(100) # some "data"
n <- 40 # generate a grid....
mx <- seq(min(x), max(x), length = n)
my <- seq(min(y), max(y), length = n)
gx <- rep(mx, n)
gy <- rep(my, rep(n, n))
tf <- exclude.too.far(gx, gy, x, y, 0.1)
plot(gx[!tf], gy[!tf], pch = ".")
points(x, y, col=2)

In the figure, the black points are the grid locations within 0.1 distance (after translating the grid to the unit square; so think of the 0.1 as being about 10% of the range of the data) of the observations, which are shown as large red points.

In the persp() setting, you have x, y, and z. Proceed as above to define the grid (the same grid that you used to create the persp() plot) and then use exclude.too.far() to identify which grid points lie too far from the data. Set the values in z to NA where exclude.too.far() returns TRUE. Then plot the modified z using persp().