Axes labels always in code for spherical datasets, independent of "unit_system"

[n-claes]

Bug report

Bug summary

@neutrinoceros discovered this when we went over #2728, what is happening is that the axes units are always in code units, independent of unit_system. It turns out that this is only the case for spherical datasets, and that it is not amrvac frontend-related. I managed to reproduce this issue with the KeplerianDisk dataset on the yt Hub as well, which is the athena_pp frontend.

Code for reproduction

blast = "amrvac/bw_2d0000.dat"
disk = "KeplerianDisk/disk.out1.00000.athdf"
for file, name in zip((blast, disk), ("amrvac", "athenapp")):
    ds = yt.load(file, unit_system="cgs")
    p = yt.SlicePlot(ds, "phi", "density")
    p.save(f"slice_cgs_{name}.png")

Actual outcome
AMRVAC dataset:

AthenaPP dataset:

Expected outcome
R and z axis labels in "cm".
For reference, this is a 3D Cartesian dataset ("amrvac/bw_3d0000.dat") sliced along y, same code as above:

which behaves as it should.

[neutrinoceros]

Reporting that the images produced by the example script have gotten better since #3533

... but this issue is still relevant.