To do precisely what you are describing, you likely want to interpolate your grid data to a set of discrete spheres with prescribed radii, then perform the spherical transform on each set of data.
Another possibility, and a generalization of the above, would be to define a spherical basis set with the spherical harmonics as the $(\theta, \phi)$-varying functions, and some sort of radial basis function defined on $r \in \left[0,a\right]$, where $a$ is the radius of the largest sphere that will fit inside your grid box. You can then solve for the coefficients of these functions via a point-matching testing procedure and inverting the coefficient matrix. This will give you a continuous function whose value coincides with your grid data at the grid points and interpolates the grid data between grid points, and which you can evaluate at each radius of interest to give you the data you want.
I don't know if there is an extant piece of software that suits your exact needs, but both approaches I described above are easily implemented in a few lines of Matlab or Python code.