c2x and Spherical Plots
In addition to cylindrical averaging, c2x from version 2.40d can also perform spherical averaging.
This example considers a hydrogen atom in a cubic box. It does use a pseudopotential. Although there are no core electrons to pseudise away, the pseudopotential will still make the Coulomb potential broader and less peaked in order to reduce its high frequency Fourier components, so the result will not be the same as the analytic result for an electron in a Coulomb potential.
$ c2x -cv -P='(0.5,0.5,0.5):R0:100' --gnu H.castep_bin H.gnu $ gnuplot H.gnu

The above plot is the volumetric charge density. Unsurprisingly, it peaks at the point where the H atom is.
The syntax to the plot command is the usual '-P=' for a line plot, the centre point, here '(0.5,0.5,0.5)', a capital R to denote a spherical radial plot, a length of zero to denote the maximum radius possible given the cell, here a 12A cube, and finally 100 points.
If the quantity plotted is a density, as it is here, sometimes it is helpful to multiply it by 4πr2 to produce a linear density. C2x can produce such a weighted plot. One simply adds a suffix of "w" after the number of points. Here we have also specified the centre as being the position of the unique H atom, rather than writing out (0.5,0.5,0.5).
$ c2x -cv -P=H:R0:100w --gnu H.castep_bin H.gnu $ gnuplot H.gnu

Now one can see that the most likely radius at which to find the electron is about 0.5A, and the density appears to extend significantly further from the nucleus. The units of this plot are now e/A, not e/A3.
Finally c2x can also plot the cumulative density by integrating after weighting. For this one replaces the "w" by an "a".

The total enclosed charge is seen to be one, unsurprisingly, and 97% of the charge is found within 2A of the atom.
One can also use "w" and "a" with cylindrical plots, where the weighting is then 2πr, not 4πr2.
Note that a potential source of error is the weighting of a density in e/Bohr3 by a radius in Angstroms. C2x attempts to store densities in A-3 internally, but it may become confused, particularly after reading files which lack a strong convention for length units. If used with -v, c2x always reports the integral of the quantity over the volume used. For densities, this can be a useful sanity check.
Details
The data are kept on whatever grid has been provided, perhaps after interpolation by an -i option. The spherical surface is sampled by points which are approximately as far apart as the samples output in the radial direction. A Fibonacci spiral is used to sample each spherical shell, as this performs unbiased sampling and can be constructed for an arbitrary number of points. Linear interpolation is then used to find the density at each sampling point. This can be slow, and note that the number of samples required scales as the cube of the number of points in the final output. No attempt to do any intelligent adaptive sampling is made.