Making a Stacked 3d Plot (and some Interpolation!)
I have a backlog of little bits of code or test projects I’ve helped folks with over the last few (or more) years. In trying to write more and get some more posts up that may be helpful, I’m going to try and post some shorter notes with some code, mainly as a reference for myself, and maybe as help for someone else.
This post is going to focus on the following:
- How do you make a nice interactive 3d plot in R?
- How do you interpolate a surface for a 3d plot?
- How might you stack that surface (like a bunch of pancakes) and plot it?
While this may seem like an odd conundrum, let me explain why it came about. Another researcher Jorge Andres Morande, asked how he might be able to interpolate a surface from drone data flown over agricultural fields at different heights. In particular, could you visualize these surfaces for different variables which were measured at different heights, like relative humidity, evapotranspiration, temperature, etc. See a short report here if you are interested in more details.
It’s an interesting question, and one that seemed fun to tackle.
Plotting with plotly
I opted to try using a plotly plot, because it is flexible, easy to play with, and seemed one of the better 3d plotting options in R. The drone data provided us with X/Y coordinates, an altitude (height), and some variables like temperature and relative humidity (rh). For this exercise, I’m just going to make up some data, assuming a spiral flight, and some rnorm points for the weather metric.
Generating Data
This code generates data for flight paths at 3 separate “levels” or heights, with random “relative humidity”" data sampled at each level.
suppressPackageStartupMessages(library(tidyverse)) # dplyr and other tools
# make fake data
t <- seq(0, 10, by=0.01)
x <- t*cos(t)
y <- t*sin(t)
z1 <- rnorm(n = 1001, mean = 5, sd = .5)
rh1 <- rnorm(n=1001, mean=34, sd=1)
z2 <- rnorm(n = 1001, mean = 10, sd = .5)
rh2 <- rnorm(n=1001, mean=33.9, sd=1)
z3 <- rnorm(n = 1001, mean = 15, sd = .5)
rh3 <- rnorm(n=1001, mean=33.8, sd=1)
dfw1 <- data.frame(x, y, z=z1, rh=rh1, lev=1)
dfw2 <- data.frame(x, y, z=z2, rh=rh2, lev=2)
dfw3 <- data.frame(x, y, z=z3, rh=rh3, lev=3)
dfw <- bind_rows(dfw1, dfw2, dfw3)Plotting in 3d
Turns out it’s not very difficult…a quick google search turned up a help page, which then helped me create the following. I show what the 2d plot looks like, then the 3D plot.
library(plotly)
library(viridis)
# first this is what it looks like in 2d
# 2d plot
ggplot(data=dfw, aes(x=x,y=y, group=lev))+
geom_path() + theme_bw() +
labs(title = "A Spiral Flight Path in 2D")
# 3d plot
plot_ly(dfw, x=~x, y=~y, z=~z) %>%
add_markers(marker=list(color=~rh,
colorscale=viridis(n = 12),
colorbar=list(title='Relative. \n Humidity'))) %>%
layout(title="Simulated Spiral Flight Path at 3 Heights",
scene = list(zaxis=list(range=c(0,20))))Interpolation
I’m not an expert. There are many reasons and methods for interpolation, and entire books written on the subject. There are multiple ways to do this, I’m showing two that worked, linear interpolation, and spline interpolation.
Linear Interpolation
The following uses akima and fields for the interpolation.
library(akima) # interpolation
library(fields) # plotting/interp
# pull out variables for one level
x1 <- dfw1$x
y1 <- dfw1$y
z1 <- dfw1$rh
# using the interp linear method:
interp1 <- interp(x1, y1, z=dfw1$rh, linear = T, nx = 90, ny = 90) # output grid 90x90
glimpse(interp1) # quick check of data## List of 3
## $ x: num [1:90] -9.48 -9.3 -9.12 -8.94 -8.77 ...
## $ y: num [1:90] -5.44 -5.29 -5.14 -4.99 -4.84 ...
## $ z: num [1:90, 1:90] NA NA NA NA NA NA NA NA NA NA ...We can then view our newly interpolated layer using the raster package, and convert the data to a raster, and extract the XY points back out from the interpolated data.
library(raster) # for converting to raster
# Raster plot
image.plot(interp1, nlevel = 20, col = viridis(n=20))
points(x1, y1, bg="white", col="gray80",pch=21, cex=0.5) #add the original points
# Contour plot: Nice, but can't overlay points/layers
filled.contour(interp1$x,interp1$y,interp1$z, color.palette = viridis)