TL;DR There are many options. Three are explained below, including one with plotly
as requested by the OP. plotly
is probably not the best. Some data preparation may be required as is the case in the example data given by both the OP and the answer by @rar.
Disclaimer: Unfortunately 3d surface plots are a bit of a pain in R. Python is a possible alternative; I find python-matplotlib
is extremely easy, straightforward, and intuitive compared to all of the many options in R.
Data preparation 1: You should actually have 3d data
Plotting 2d data as 3d is a common mistake. The data in the OP, for instance, forms nearly linear curve in 3d. While it may be possible to treat these points as independent observations of z-axis values in particular locations in the x-y-plain and to interpolate a 3d-surface from that, I doubt that the resulting surface plot would be very useful or informative.
In any case for plotting of 3d-surfaces, you will almost always need a grid of observations over the x-y-plain. To interpolate values necessary for such a grid in R, you may use akima::interp
like so
library(akima)
df <- data.frame(matrix(rnorm(60), nrow=20))
plot_data <- interp(df$X1,df$X2,df$X3)
However, akima::interp
is rather fragile and will crash your R interpreter shell in certain cases (e.g. for the example of the data_grid
variable below).
On the plus side, it creates a data structure with tow vectors and a matrix for you, the way all (?) 3d surface plotting functions in R expect the data.
Data preparation 2: You cannot have data in long format
Say, you have data for 9 points in a 3x3 grid (to keep this example simple):
data_grid <- data.frame(data_col = c(45.62151, 60.30996, 66.01667, 45.48701,
60.39519, 65.42441, 45.48208, 60.39041, 65.52165),
axis_one=c(10000, 10000, 10000, 1000000, 1000000, 1000000, 100000000,
100000000, 100000000),
axis_two=c(1, 100, 10000,1, 100, 10000,1, 100, 10000))
data_grid
# data_col axis_one axis_two
#1 45.62151 1e+04 1
#2 60.30996 1e+04 100
#3 66.01667 1e+04 10000
#4 45.48701 1e+06 1
#5 60.39519 1e+06 100
#6 65.42441 1e+06 10000
#7 45.48208 1e+08 1
#8 60.39041 1e+08 100
#9 65.52165 1e+08 10000
I was not able to find any 3d-plotting package or function in R that would accept a data structure like this. Most replies regarding 3d plotting in R on stackoverflow assume this, probably because in most cases it is assumed that this is not empirical data but values pulled from a 3d function. Examples here and here and here or the (for 3d plotting examples) frequently used volcano
data set included in base R.
So what do these packages expect? They want two vectors (for the grid coordinates in the base plane) and a matrix with entries corresponding to the meshgrid created by these two vectors for the z axis. In other words you need to rearrange any long format data to wide format, i.e. a matrix with the x and y coordinates as row and column indices/labels.
This can be done with reshape2::acast
:
library(reshape2)
plot_matrix <- t(acast(data_grid, axis_one~axis_two, value.var="data_col"))
plot_matrix
# 10000 1e+06 1e+08
#1 45.62151 45.48701 45.48208
#100 60.30996 60.39519 60.39041
#10000 66.01667 65.42441 65.52165
If you have data in wide format already: brilliant, you don't have to do this.
Note that the other answer by @rar seems to uses example data structured in long format but plots it the wrong way, such that R interprets columns cs
and w
not as base plain coordinates but as additional observations. It is difficult to tell if this is indeed a mistake or if the answer intended it that way; unfortunately it includes very little text.
Plotting option 1: Plotting it with persp
Plotting is straightforward enough. Unfortunately, we have modified the x and y data as well now. So we have to pull these back out of the column and row names of the matrix. There may be other ways to save these before transforming the data to a wide format matrix but this may be a source of errors if data for different axes is taken from different objects.
persp(x = as.numeric(colnames(plot_matrix)),
y = as.numeric(rownames(plot_matrix)),
z = plot_matrix,
xlab = "Axis one",
ylab = "Axis two",
zlab = "Data",
ticktype ='detailed',
theta = 310,
phi = 20,
col = "green", shade = 0.5)
A major drawback of persp
is that it does not seem to be able to have log scale axes. This is inconvenient; so we have to do it manually.
persp(x = log(as.numeric(colnames(plot_matrix))),
y = log(as.numeric(rownames(plot_matrix))),
z = plot_matrix,
xlab = "log Axis one",
ylab = "log Axis two",
zlab = "Data",
ticktype ='detailed',
theta = 310,
phi = 20,
col = "green", shade = 0.5)
We save the plot as usual (albeit quite counterintuitive compared to, say, Python)
pdf(file="out.pdf", width=5, height=5)
persp(x = log(as.numeric(colnames(plot_matrix))),
y = log(as.numeric(rownames(plot_matrix))),
z = plot_matrix,
xlab = "log Axis one",
ylab = "log Axis two",
zlab = "Data",
ticktype ='detailed',
theta = 310,
phi = 20,
col = "green", shade = 0.5)
dev.off()
...or equivalently for png etc. etc.

Plotting option 2: Plotting it with plotly
The OP asked about plotly instead. So here we go:
Plotly has more functionality and looks nicer, but retains some of the same drawbacks (requires wide format).
library(plotly)
plot_ly(
x = as.numeric(colnames(plot_matrix)),
y = as.numeric(rownames(plot_matrix)),
z = plot_matrix
) %>%
add_surface() %>%
layout(
title = "",
scene = list(
xaxis = list(type = "log", title = "Total observations"),
yaxis = list(type = "log", title = "Firm size"),
zaxis = list(title = "Median"),
camera = list(eye = list(x = 1.95, y = -1.25, z = 1.25))
))
What may be near impossible with plotly, however, is saving the figure. In the past, plotly seems to have required subscription to a paid service to create a figure online on their server with plotly::export
. Today it recommends using a package called orca
which has to be installed. The installation options are described on orca's github repository and invariably require you to hack your system by installing either untrusted third party packages or a second package manager conda
which may also want root access or an entire second Python distribution.
If you are using their tool stack already, e.g., if you have an Anaconda distribution for Python or so, you may be able to install orca
without compromising your system's integrity. Otherwise, don't bother. in this case, you can save the plot (apparently; I can't test it) by doing
plotly_figure <- plot_ly(
x = as.numeric(colnames(plot_matrix)),
y = as.numeric(rownames(plot_matrix)),
z = plot_matrix
) %>%
add_surface() %>%
layout(
title = "",
scene = list(
xaxis = list(type = "log", title = "Total observations"),
yaxis = list(type = "log", title = "Firm size"),
zaxis = list(title = "Median"),
camera = list(eye = list(x = 1.95, y = -1.25, z = 1.25))
))
orca(plotly_figure, file="out.pdf")
As a side note, their website's employment of third party scripts that will break browser add-ons (privacy enhancing ones in particular) is also a bit of a red flag.
The figure is, of course more like a screenshot as it won't permit me to save anything.
Plotting option 3+: Other packages
There are some more packages. E.g. plot3D
, which has some functions including plot3D::persp3D
, the one here appropriate. Compared to persp
, this looks a bit more fancy and adds colors highlighting surface elevation.
plot3D::persp3D(x = log(as.numeric(colnames(plot_matrix))),
y = log(as.numeric(rownames(plot_matrix))),
z = plot_matrix,
xlab = "log Axis one",
ylab = "log Axis two",
zlab = "Data",
ticktype ='detailed',
theta = 310,
phi = 20)
It looks more interesting for more dense data like the data created above with akima::interp
:
plot3D::persp3D(x = plot_data$x,
y = plot_data$y,
z = plot_data$z,
xlab = "log Axis one",
ylab = "log Axis two",
zlab = "Data",
ticktype ='detailed',
theta = 310,
phi = 20)
Saving the figures should work like in persp
above.

There is another one rgl:plot3d
; if I understand it correctly, this requires extremely dense data and otherwise falls back to a 3d scatter plot.