Color Scale based on different values than z column in surface plot

I am generating a 3D Kernal Density Plot of Steph Curry’s shooting data using the following code:

library(rjson)
library(jsonlite)
library(RCurl)

url = “http://stats.nba.com/stats/shotchartdetail?CFID=33&CFPARAMS=2017-18&ContextFilter=&ContextMeasure=FGA&DateFrom=&DateTo=&GameID=&GameSegment=&LastNGames=0&LeagueID=00&Location=&MeasureType=Base&Month=0&OpponentTeamID=0&Outcome=&PaceAdjust=N&PerMode=PerGame&Period=0&PlayerID=201939&PlusMinus=N&PlayerPosition=&Rank=N&RookieYear=&Season=2017-18&SeasonSegment=&SeasonType=Regular%20Season&TeamID=0&VsConference=&VsDivision=

content = getURLContent(url, verbose = TRUE, useragent = getOption(“HTTPUserAgent”))

shot.data <- fromJSON( txt = content )

steph_curry <- as.data.frame(shot.data$resultSets$rowSet[[1]], stringsAsFactors = FALSE )

colnames(steph_curry) <- shot.data$resultSets$headers[[1]]

library(plotly)
library(MASS)

steph_curry$LOC_X = as.numeric(steph_curry$LOC_X)
steph_curry$LOC_Y = as.numeric(steph_curry$LOC_Y)
steph_curry$SHOT_MADE_FLAG = as.numeric(steph_curry$SHOT_MADE_FLAG)

den3d <- kde2d(steph_curry$LOC_X,steph_curry$LOC_Y)

plot_ly(x=den3d$x, y=den3d$y, z=den3d$z) %>% add_surface()

All of the above code works just fine. Next what I am wanting to do is color code the graph based on field goal percentage rather than the z value of the kernal density. I use the following code to generate a dataframe of these field goal percentages by create hexbins of the shooting data where each hexbins’ z value is shooting percentage, and then assigning each x, y pairing in my kernal density the z value of the hexbin that is the minimum cartesian distance away. Again, I am able to successfully do this with the following code:

playermeans = den3d$z
row.names(playermeans) = den3d$x
colnames(playermeans) = den3d$y

a = hexbin(x = steph_curry$LOC_X,y=steph_curry$LOC_Y,IDs = TRUE,xbins = 10)
b = hexTapply(a,steph_curry$SHOT_MADE_FLAG,FUN = mean)

hexmeans=NULL
hexmeans$X = a@xcm
hexmeans = as.data.frame(means)
hexmeans$Y = a@ycm
hexmeans$Z = b

playermeans = as.data.frame(playermeans)

for(i in 1:nrow(playermeans)){
for (j in 1:ncol(playermeans)){
x_coord = as.numeric(rownames(playermeans)[i])
y_coord = as.numeric(colnames(playermeans)[j])
dists <- (x_coord - hexmeans$X)^2 + (y_coord - hexmeans$Y)^2
playermeans[i,j] = hexmeans$Z[dists == min(dists)]
}
}
The only remaining thing to do is assign a color gradient based on the shooting percentages found in the playermeans variable. Assuming it is possible, how do I do this?