In this 2011 paper (Reyes, V.M., 2011), we reported on our efforts at representing spheroprotiens using
spherical (ρ, ϕ, θ) coordinates, with the protein centroid as origin, and presented two of its potential
applications. Here, we present the Fortran program source codes we used in our procedure, starting with the
determination of the protein molecular (geometric) centroid and the transformation of its Cartesian coordinates
to spherical coordinates. We also present here the Fortran source codes for the implementation of the two
applications we presented in the aforementioned paper, namely, the separation of the protein outer layer (OL)
from its iner core (IC), and the identification of protrusions and invaginations on the protein surface.
Once you have rho, phi and theta in some matrix text file, you could perhaps import that into R and use ggplot and its stat_density_2d function to plot phi and theta along x and y axes, and use rho as the heatmap "density" value.
Lots of examples online on how to make that type of plot with sample datasets, which you'd adjust for your particular coordinates.