BUG: plot bounds are off with ds.r[...].plot in non-cartesian geometries

[neutrinoceros]

Bug report

Bug summary

This was discovered in #4497, which I reproduce here for self-consistency

Code for reproduction

import yt
import yt.testing

ds = yt.testing.fake_amr_ds(geometry="spherical")
c = ds.domain_center
ds.r[:,c[1],:].plot("Density").save("/tmp/test.png")

Actual outcome

Expected outcome

The output should look identical to

import yt
import yt.testing

ds = yt.testing.fake_amr_ds(geometry="spherical")
c = ds.domain_center
yt.SlicePlot(ds, "theta", "Density").save("/tmp/test2.png")

[neutrinoceros]

#4500 fixes the simple cases that are illustrated here and in #4497, but actual tests are needed to close this, which will be easier to manage on the dev branch (especially considering #4362)

[neutrinoceros]

I am not so sure about what the "correct" behavior should be (the current one is still clearly incorrect). Given that YTSlice.plot is already slightly different from SlicePlot in that it uses ˋorigin="native"` by default, maybe it would make sense to also use native plot axes (as #4179 attempts at enabling).
@matthewturk, @jzuhone I would be interested in your opinions on this.