We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
BUG: wrong axis labels representation in SlicePlot for spherical data
### Bug report

**Bug summary**

There's an inconsistency in math representation (or lack thereof) of axis labels when performing SlicePlot along the "r" direction for spherical data.

**Code for reproduction**


```python
import yt
ds = yt.load_sample("bw_spherical_2d")
p = yt.SlicePlot(ds, "r", "density")
p.save("/tmp/bug.png")
```

**Actual outcome**

![bug](https://user-images.githubusercontent.com/14075922/134900776-771c7787-44cb-4751-9a16-b06432be46f3.png)

**Expected outcome**

I expect axis labels to use `\theta` and `\phi` latex glyphs instead of "theta" and "phi" strings, consistently with what happens when slicing along the "theta" direction, for instance
![bug_theta](https://user-images.githubusercontent.com/14075922/134901098-a37df2fa-3b1e-4d8f-aebb-a9c3c5197864.png)

(here the frame is wrong, this was already reported as #3529)


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `yt/geometry/coordinates/spherical_coordinates.py`
Content:
```
1 import numpy as np
2 
3 from yt.utilities.lib.pixelization_routines import pixelize_aitoff, pixelize_cylinder
4 
5 from .coordinate_handler import CoordinateHandler, _get_coord_fields, _unknown_coord
6 
7 
8 class SphericalCoordinateHandler(CoordinateHandler):
9     name = "spherical"
10 
11     def __init__(self, ds, ordering=("r", "theta", "phi")):
12         super().__init__(ds, ordering)
13         # Generate
14         self.image_units = {}
15         self.image_units[self.axis_id["r"]] = ("rad", "rad")
16         self.image_units[self.axis_id["theta"]] = (None, None)
17         self.image_units[self.axis_id["phi"]] = (None, None)
18 
19     def setup_fields(self, registry):
20         # return the fields for r, z, theta
21         registry.add_field(
22             ("index", "dx"), sampling_type="cell", function=_unknown_coord
23         )
24 
25         registry.add_field(
26             ("index", "dy"), sampling_type="cell", function=_unknown_coord
27         )
28 
29         registry.add_field(
30             ("index", "dz"), sampling_type="cell", function=_unknown_coord
31         )
32 
33         registry.add_field(
34             ("index", "x"), sampling_type="cell", function=_unknown_coord
35         )
36 
37         registry.add_field(
38             ("index", "y"), sampling_type="cell", function=_unknown_coord
39         )
40 
41         registry.add_field(
42             ("index", "z"), sampling_type="cell", function=_unknown_coord
43         )
44 
45         f1, f2 = _get_coord_fields(self.axis_id["r"])
46         registry.add_field(
47             ("index", "dr"),
48             sampling_type="cell",
49             function=f1,
50             display_field=False,
51             units="code_length",
52         )
53 
54         registry.add_field(
55             ("index", "r"),
56             sampling_type="cell",
57             function=f2,
58             display_field=False,
59             units="code_length",
60         )
61 
62         f1, f2 = _get_coord_fields(self.axis_id["theta"], "")
63         registry.add_field(
64             ("index", "dtheta"),
65             sampling_type="cell",
66             function=f1,
67             display_field=False,
68             units="",
69         )
70 
71         registry.add_field(
72             ("index", "theta"),
73             sampling_type="cell",
74             function=f2,
75             display_field=False,
76             units="",
77         )
78 
79         f1, f2 = _get_coord_fields(self.axis_id["phi"], "")
80         registry.add_field(
81             ("index", "dphi"),
82             sampling_type="cell",
83             function=f1,
84             display_field=False,
85             units="",
86         )
87 
88         registry.add_field(
89             ("index", "phi"),
90             sampling_type="cell",
91             function=f2,
92             display_field=False,
93             units="",
94         )
95 
96         def _SphericalVolume(field, data):
97             # Here we compute the spherical volume element exactly
98             r = data["index", "r"]
99             dr = data["index", "dr"]
100             theta = data["index", "theta"]
101             dtheta = data["index", "dtheta"]
102             vol = ((r + 0.5 * dr) ** 3 - (r - 0.5 * dr) ** 3) / 3.0
103             vol *= np.cos(theta - 0.5 * dtheta) - np.cos(theta + 0.5 * dtheta)
104             vol *= data["index", "dphi"]
105             return vol
106 
107         registry.add_field(
108             ("index", "cell_volume"),
109             sampling_type="cell",
110             function=_SphericalVolume,
111             units="code_length**3",
112         )
113         registry.alias(("index", "volume"), ("index", "cell_volume"))
114 
115         def _path_r(field, data):
116             return data["index", "dr"]
117 
118         registry.add_field(
119             ("index", "path_element_r"),
120             sampling_type="cell",
121             function=_path_r,
122             units="code_length",
123         )
124 
125         def _path_theta(field, data):
126             # Note: this already assumes cell-centered
127             return data["index", "r"] * data["index", "dtheta"]
128 
129         registry.add_field(
130             ("index", "path_element_theta"),
131             sampling_type="cell",
132             function=_path_theta,
133             units="code_length",
134         )
135 
136         def _path_phi(field, data):
137             # Note: this already assumes cell-centered
138             return (
139                 data["index", "r"]
140                 * data["index", "dphi"]
141                 * np.sin(data["index", "theta"])
142             )
143 
144         registry.add_field(
145             ("index", "path_element_phi"),
146             sampling_type="cell",
147             function=_path_phi,
148             units="code_length",
149         )
150 
151     def pixelize(
152         self, dimension, data_source, field, bounds, size, antialias=True, periodic=True
153     ):
154         self.period
155         name = self.axis_name[dimension]
156         if name == "r":
157             return self._ortho_pixelize(
158                 data_source, field, bounds, size, antialias, dimension, periodic
159             )
160         elif name in ("theta", "phi"):
161             return self._cyl_pixelize(
162                 data_source, field, bounds, size, antialias, dimension
163             )
164         else:
165             raise NotImplementedError
166 
167     def pixelize_line(self, field, start_point, end_point, npoints):
168         raise NotImplementedError
169 
170     def _ortho_pixelize(
171         self, data_source, field, bounds, size, antialias, dim, periodic
172     ):
173         buff = pixelize_aitoff(
174             data_source["py"],
175             data_source["pdy"],
176             data_source["px"],
177             data_source["pdx"],
178             size,
179             data_source[field],
180             None,
181             None,
182             theta_offset=0,
183             phi_offset=0,
184         ).transpose()
185         return buff
186 
187     def _cyl_pixelize(self, data_source, field, bounds, size, antialias, dimension):
188         name = self.axis_name[dimension]
189         buff = np.full((size[1], size[0]), np.nan, dtype="f8")
190         if name == "theta":
191             pixelize_cylinder(
192                 buff,
193                 data_source["px"],
194                 data_source["pdx"],
195                 data_source["py"],
196                 data_source["pdy"],
197                 data_source[field],
198                 bounds,
199             )
200         elif name == "phi":
201             # Note that we feed in buff.T here
202             pixelize_cylinder(
203                 buff.T,
204                 data_source["px"],
205                 data_source["pdx"],
206                 data_source["py"],
207                 data_source["pdy"],
208                 data_source[field],
209                 bounds,
210             )
211         else:
212             raise RuntimeError
213         self.sanitize_buffer_fill_values(buff)
214         return buff
215 
216     def convert_from_cartesian(self, coord):
217         raise NotImplementedError
218 
219     def convert_to_cartesian(self, coord):
220         if isinstance(coord, np.ndarray) and len(coord.shape) > 1:
221             ri = self.axis_id["r"]
222             thetai = self.axis_id["theta"]
223             phii = self.axis_id["phi"]
224             r = coord[:, ri]
225             theta = coord[:, thetai]
226             phi = coord[:, phii]
227             nc = np.zeros_like(coord)
228             # r, theta, phi
229             nc[:, ri] = np.cos(phi) * np.sin(theta) * r
230             nc[:, thetai] = np.sin(phi) * np.sin(theta) * r
231             nc[:, phii] = np.cos(theta) * r
232         else:
233             r, theta, phi = coord
234             nc = (
235                 np.cos(phi) * np.sin(theta) * r,
236                 np.sin(phi) * np.sin(theta) * r,
237                 np.cos(theta) * r,
238             )
239         return nc
240 
241     def convert_to_cylindrical(self, coord):
242         raise NotImplementedError
243 
244     def convert_from_cylindrical(self, coord):
245         raise NotImplementedError
246 
247     def convert_to_spherical(self, coord):
248         raise NotImplementedError
249 
250     def convert_from_spherical(self, coord):
251         raise NotImplementedError
252 
253     _image_axis_name = None
254 
255     @property
256     def image_axis_name(self):
257         if self._image_axis_name is not None:
258             return self._image_axis_name
259         # This is the x and y axes labels that get displayed.  For
260         # non-Cartesian coordinates, we usually want to override these for
261         # Cartesian coordinates, since we transform them.
262         rv = {
263             self.axis_id["r"]: ("theta", "phi"),
264             self.axis_id["theta"]: ("x / \\sin(\\theta)", "y / \\sin(\\theta)"),
265             self.axis_id["phi"]: ("R", "z"),
266         }
267         for i in list(rv.keys()):
268             rv[self.axis_name[i]] = rv[i]
269             rv[self.axis_name[i].capitalize()] = rv[i]
270         self._image_axis_name = rv
271         return rv
272 
273     _x_pairs = (("r", "theta"), ("theta", "r"), ("phi", "r"))
274     _y_pairs = (("r", "phi"), ("theta", "phi"), ("phi", "theta"))
275 
276     @property
277     def period(self):
278         return self.ds.domain_width
279 
280     def sanitize_center(self, center, axis):
281         center, display_center = super().sanitize_center(center, axis)
282         name = self.axis_name[axis]
283         if name == "r":
284             display_center = center
285         elif name == "theta":
286             display_center = (
287                 0.0 * display_center[0],
288                 0.0 * display_center[1],
289                 0.0 * display_center[2],
290             )
291         elif name == "phi":
292             display_center = [
293                 self.ds.domain_width[0] / 2.0 + self.ds.domain_left_edge[0],
294                 0.0 * display_center[1],
295                 0.0 * display_center[2],
296             ]
297             ri = self.axis_id["r"]
298             c = self.ds.domain_width[ri] / 2.0 + self.ds.domain_left_edge[ri]
299             display_center[ri] = c
300             display_center = tuple(display_center)
301         return center, display_center
302 
303     def sanitize_width(self, axis, width, depth):
304         name = self.axis_name[axis]
305         if width is not None:
306             width = super().sanitize_width(axis, width, depth)
307         elif name == "r":
308             width = [
309                 self.ds.domain_width[self.x_axis["r"]],
310                 self.ds.domain_width[self.y_axis["r"]],
311             ]
312         elif name == "theta":
313             ri = self.axis_id["r"]
314             # Remember, in spherical coordinates when we cut in theta,
315             # we create a conic section
316             width = [2.0 * self.ds.domain_width[ri], 2.0 * self.ds.domain_width[ri]]
317         elif name == "phi":
318             ri = self.axis_id["r"]
319             width = [self.ds.domain_right_edge[ri], 2.0 * self.ds.domain_width[ri]]
320         return width
321 
322     def _sanity_check(self):
323         """This prints out a handful of diagnostics that help verify the
324         dataset is well formed."""
325         # We just check a few things here.
326         dd = self.ds.all_data()
327         r0 = self.ds.domain_left_edge[self.axis_id["r"]]
328         r1 = self.ds.domain_right_edge[self.axis_id["r"]]
329         v1 = 4.0 * np.pi / 3.0 * (r1 ** 3 - r0 ** 3)
330         print(f"Total volume should be 4*pi*r**3 = {v1:0.16e}")
331         v2 = dd.quantities.total_quantity("cell_volume")
332         print(f"Actual volume is                   {v2:0.16e}")
333         print(f"Relative difference: {np.abs(v2 - v1) / (v2 + v1):0.16e}")
334
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/yt/geometry/coordinates/spherical_coordinates.py b/yt/geometry/coordinates/spherical_coordinates.py
--- a/yt/geometry/coordinates/spherical_coordinates.py
+++ b/yt/geometry/coordinates/spherical_coordinates.py
@@ -260,7 +260,7 @@
         # non-Cartesian coordinates, we usually want to override these for
         # Cartesian coordinates, since we transform them.
         rv = {
-            self.axis_id["r"]: ("theta", "phi"),
+            self.axis_id["r"]: ("\\theta", "\\phi"),
             self.axis_id["theta"]: ("x / \\sin(\\theta)", "y / \\sin(\\theta)"),
             self.axis_id["phi"]: ("R", "z"),
         }

{"golden_diff": "diff --git a/yt/geometry/coordinates/spherical_coordinates.py b/yt/geometry/coordinates/spherical_coordinates.py\n--- a/yt/geometry/coordinates/spherical_coordinates.py\n+++ b/yt/geometry/coordinates/spherical_coordinates.py\n@@ -260,7 +260,7 @@\n         # non-Cartesian coordinates, we usually want to override these for\n         # Cartesian coordinates, since we transform them.\n         rv = {\n-            self.axis_id[\"r\"]: (\"theta\", \"phi\"),\n+            self.axis_id[\"r\"]: (\"\\\\theta\", \"\\\\phi\"),\n             self.axis_id[\"theta\"]: (\"x / \\\\sin(\\\\theta)\", \"y / \\\\sin(\\\\theta)\"),\n             self.axis_id[\"phi\"]: (\"R\", \"z\"),\n         }\n", "issue": "BUG: wrong axis labels representation in SlicePlot for spherical data\n### Bug report\r\n\r\n**Bug summary**\r\n\r\nThere's an inconsistency in math representation (or lack thereof) of axis labels when performing SlicePlot along the \"r\" direction for spherical data.\r\n\r\n**Code for reproduction**\r\n\r\n\r\n```python\r\nimport yt\r\nds = yt.load_sample(\"bw_spherical_2d\")\r\np = yt.SlicePlot(ds, \"r\", \"density\")\r\np.save(\"/tmp/bug.png\")\r\n```\r\n\r\n**Actual outcome**\r\n\r\n![bug](https://user-images.githubusercontent.com/14075922/134900776-771c7787-44cb-4751-9a16-b06432be46f3.png)\r\n\r\n**Expected outcome**\r\n\r\nI expect axis labels to use `\\theta` and `\\phi` latex glyphs instead of \"theta\" and \"phi\" strings, consistently with what happens when slicing along the \"theta\" direction, for instance\r\n![bug_theta](https://user-images.githubusercontent.com/14075922/134901098-a37df2fa-3b1e-4d8f-aebb-a9c3c5197864.png)\r\n\r\n(here the frame is wrong, this was already reported as #3529)\r\n\n", "before_files": [{"content": "import numpy as np\n\nfrom yt.utilities.lib.pixelization_routines import pixelize_aitoff, pixelize_cylinder\n\nfrom .coordinate_handler import CoordinateHandler, _get_coord_fields, _unknown_coord\n\n\nclass SphericalCoordinateHandler(CoordinateHandler):\n    name = \"spherical\"\n\n    def __init__(self, ds, ordering=(\"r\", \"theta\", \"phi\")):\n        super().__init__(ds, ordering)\n        # Generate\n        self.image_units = {}\n        self.image_units[self.axis_id[\"r\"]] = (\"rad\", \"rad\")\n        self.image_units[self.axis_id[\"theta\"]] = (None, None)\n        self.image_units[self.axis_id[\"phi\"]] = (None, None)\n\n    def setup_fields(self, registry):\n        # return the fields for r, z, theta\n        registry.add_field(\n            (\"index\", \"dx\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        registry.add_field(\n            (\"index\", \"dy\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        registry.add_field(\n            (\"index\", \"dz\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        registry.add_field(\n            (\"index\", \"x\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        registry.add_field(\n            (\"index\", \"y\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        registry.add_field(\n            (\"index\", \"z\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        f1, f2 = _get_coord_fields(self.axis_id[\"r\"])\n        registry.add_field(\n            (\"index\", \"dr\"),\n            sampling_type=\"cell\",\n            function=f1,\n            display_field=False,\n            units=\"code_length\",\n        )\n\n        registry.add_field(\n            (\"index\", \"r\"),\n            sampling_type=\"cell\",\n            function=f2,\n            display_field=False,\n            units=\"code_length\",\n        )\n\n        f1, f2 = _get_coord_fields(self.axis_id[\"theta\"], \"\")\n        registry.add_field(\n            (\"index\", \"dtheta\"),\n            sampling_type=\"cell\",\n            function=f1,\n            display_field=False,\n            units=\"\",\n        )\n\n        registry.add_field(\n            (\"index\", \"theta\"),\n            sampling_type=\"cell\",\n            function=f2,\n            display_field=False,\n            units=\"\",\n        )\n\n        f1, f2 = _get_coord_fields(self.axis_id[\"phi\"], \"\")\n        registry.add_field(\n            (\"index\", \"dphi\"),\n            sampling_type=\"cell\",\n            function=f1,\n            display_field=False,\n            units=\"\",\n        )\n\n        registry.add_field(\n            (\"index\", \"phi\"),\n            sampling_type=\"cell\",\n            function=f2,\n            display_field=False,\n            units=\"\",\n        )\n\n        def _SphericalVolume(field, data):\n            # Here we compute the spherical volume element exactly\n            r = data[\"index\", \"r\"]\n            dr = data[\"index\", \"dr\"]\n            theta = data[\"index\", \"theta\"]\n            dtheta = data[\"index\", \"dtheta\"]\n            vol = ((r + 0.5 * dr) ** 3 - (r - 0.5 * dr) ** 3) / 3.0\n            vol *= np.cos(theta - 0.5 * dtheta) - np.cos(theta + 0.5 * dtheta)\n            vol *= data[\"index\", \"dphi\"]\n            return vol\n\n        registry.add_field(\n            (\"index\", \"cell_volume\"),\n            sampling_type=\"cell\",\n            function=_SphericalVolume,\n            units=\"code_length**3\",\n        )\n        registry.alias((\"index\", \"volume\"), (\"index\", \"cell_volume\"))\n\n        def _path_r(field, data):\n            return data[\"index\", \"dr\"]\n\n        registry.add_field(\n            (\"index\", \"path_element_r\"),\n            sampling_type=\"cell\",\n            function=_path_r,\n            units=\"code_length\",\n        )\n\n        def _path_theta(field, data):\n            # Note: this already assumes cell-centered\n            return data[\"index\", \"r\"] * data[\"index\", \"dtheta\"]\n\n        registry.add_field(\n            (\"index\", \"path_element_theta\"),\n            sampling_type=\"cell\",\n            function=_path_theta,\n            units=\"code_length\",\n        )\n\n        def _path_phi(field, data):\n            # Note: this already assumes cell-centered\n            return (\n                data[\"index\", \"r\"]\n                * data[\"index\", \"dphi\"]\n                * np.sin(data[\"index\", \"theta\"])\n            )\n\n        registry.add_field(\n            (\"index\", \"path_element_phi\"),\n            sampling_type=\"cell\",\n            function=_path_phi,\n            units=\"code_length\",\n        )\n\n    def pixelize(\n        self, dimension, data_source, field, bounds, size, antialias=True, periodic=True\n    ):\n        self.period\n        name = self.axis_name[dimension]\n        if name == \"r\":\n            return self._ortho_pixelize(\n                data_source, field, bounds, size, antialias, dimension, periodic\n            )\n        elif name in (\"theta\", \"phi\"):\n            return self._cyl_pixelize(\n                data_source, field, bounds, size, antialias, dimension\n            )\n        else:\n            raise NotImplementedError\n\n    def pixelize_line(self, field, start_point, end_point, npoints):\n        raise NotImplementedError\n\n    def _ortho_pixelize(\n        self, data_source, field, bounds, size, antialias, dim, periodic\n    ):\n        buff = pixelize_aitoff(\n            data_source[\"py\"],\n            data_source[\"pdy\"],\n            data_source[\"px\"],\n            data_source[\"pdx\"],\n            size,\n            data_source[field],\n            None,\n            None,\n            theta_offset=0,\n            phi_offset=0,\n        ).transpose()\n        return buff\n\n    def _cyl_pixelize(self, data_source, field, bounds, size, antialias, dimension):\n        name = self.axis_name[dimension]\n        buff = np.full((size[1], size[0]), np.nan, dtype=\"f8\")\n        if name == \"theta\":\n            pixelize_cylinder(\n                buff,\n                data_source[\"px\"],\n                data_source[\"pdx\"],\n                data_source[\"py\"],\n                data_source[\"pdy\"],\n                data_source[field],\n                bounds,\n            )\n        elif name == \"phi\":\n            # Note that we feed in buff.T here\n            pixelize_cylinder(\n                buff.T,\n                data_source[\"px\"],\n                data_source[\"pdx\"],\n                data_source[\"py\"],\n                data_source[\"pdy\"],\n                data_source[field],\n                bounds,\n            )\n        else:\n            raise RuntimeError\n        self.sanitize_buffer_fill_values(buff)\n        return buff\n\n    def convert_from_cartesian(self, coord):\n        raise NotImplementedError\n\n    def convert_to_cartesian(self, coord):\n        if isinstance(coord, np.ndarray) and len(coord.shape) > 1:\n            ri = self.axis_id[\"r\"]\n            thetai = self.axis_id[\"theta\"]\n            phii = self.axis_id[\"phi\"]\n            r = coord[:, ri]\n            theta = coord[:, thetai]\n            phi = coord[:, phii]\n            nc = np.zeros_like(coord)\n            # r, theta, phi\n            nc[:, ri] = np.cos(phi) * np.sin(theta) * r\n            nc[:, thetai] = np.sin(phi) * np.sin(theta) * r\n            nc[:, phii] = np.cos(theta) * r\n        else:\n            r, theta, phi = coord\n            nc = (\n                np.cos(phi) * np.sin(theta) * r,\n                np.sin(phi) * np.sin(theta) * r,\n                np.cos(theta) * r,\n            )\n        return nc\n\n    def convert_to_cylindrical(self, coord):\n        raise NotImplementedError\n\n    def convert_from_cylindrical(self, coord):\n        raise NotImplementedError\n\n    def convert_to_spherical(self, coord):\n        raise NotImplementedError\n\n    def convert_from_spherical(self, coord):\n        raise NotImplementedError\n\n    _image_axis_name = None\n\n    @property\n    def image_axis_name(self):\n        if self._image_axis_name is not None:\n            return self._image_axis_name\n        # This is the x and y axes labels that get displayed.  For\n        # non-Cartesian coordinates, we usually want to override these for\n        # Cartesian coordinates, since we transform them.\n        rv = {\n            self.axis_id[\"r\"]: (\"theta\", \"phi\"),\n            self.axis_id[\"theta\"]: (\"x / \\\\sin(\\\\theta)\", \"y / \\\\sin(\\\\theta)\"),\n            self.axis_id[\"phi\"]: (\"R\", \"z\"),\n        }\n        for i in list(rv.keys()):\n            rv[self.axis_name[i]] = rv[i]\n            rv[self.axis_name[i].capitalize()] = rv[i]\n        self._image_axis_name = rv\n        return rv\n\n    _x_pairs = ((\"r\", \"theta\"), (\"theta\", \"r\"), (\"phi\", \"r\"))\n    _y_pairs = ((\"r\", \"phi\"), (\"theta\", \"phi\"), (\"phi\", \"theta\"))\n\n    @property\n    def period(self):\n        return self.ds.domain_width\n\n    def sanitize_center(self, center, axis):\n        center, display_center = super().sanitize_center(center, axis)\n        name = self.axis_name[axis]\n        if name == \"r\":\n            display_center = center\n        elif name == \"theta\":\n            display_center = (\n                0.0 * display_center[0],\n                0.0 * display_center[1],\n                0.0 * display_center[2],\n            )\n        elif name == \"phi\":\n            display_center = [\n                self.ds.domain_width[0] / 2.0 + self.ds.domain_left_edge[0],\n                0.0 * display_center[1],\n                0.0 * display_center[2],\n            ]\n            ri = self.axis_id[\"r\"]\n            c = self.ds.domain_width[ri] / 2.0 + self.ds.domain_left_edge[ri]\n            display_center[ri] = c\n            display_center = tuple(display_center)\n        return center, display_center\n\n    def sanitize_width(self, axis, width, depth):\n        name = self.axis_name[axis]\n        if width is not None:\n            width = super().sanitize_width(axis, width, depth)\n        elif name == \"r\":\n            width = [\n                self.ds.domain_width[self.x_axis[\"r\"]],\n                self.ds.domain_width[self.y_axis[\"r\"]],\n            ]\n        elif name == \"theta\":\n            ri = self.axis_id[\"r\"]\n            # Remember, in spherical coordinates when we cut in theta,\n            # we create a conic section\n            width = [2.0 * self.ds.domain_width[ri], 2.0 * self.ds.domain_width[ri]]\n        elif name == \"phi\":\n            ri = self.axis_id[\"r\"]\n            width = [self.ds.domain_right_edge[ri], 2.0 * self.ds.domain_width[ri]]\n        return width\n\n    def _sanity_check(self):\n        \"\"\"This prints out a handful of diagnostics that help verify the\n        dataset is well formed.\"\"\"\n        # We just check a few things here.\n        dd = self.ds.all_data()\n        r0 = self.ds.domain_left_edge[self.axis_id[\"r\"]]\n        r1 = self.ds.domain_right_edge[self.axis_id[\"r\"]]\n        v1 = 4.0 * np.pi / 3.0 * (r1 ** 3 - r0 ** 3)\n        print(f\"Total volume should be 4*pi*r**3 = {v1:0.16e}\")\n        v2 = dd.quantities.total_quantity(\"cell_volume\")\n        print(f\"Actual volume is                   {v2:0.16e}\")\n        print(f\"Relative difference: {np.abs(v2 - v1) / (v2 + v1):0.16e}\")\n", "path": "yt/geometry/coordinates/spherical_coordinates.py"}], "after_files": [{"content": "import numpy as np\n\nfrom yt.utilities.lib.pixelization_routines import pixelize_aitoff, pixelize_cylinder\n\nfrom .coordinate_handler import CoordinateHandler, _get_coord_fields, _unknown_coord\n\n\nclass SphericalCoordinateHandler(CoordinateHandler):\n    name = \"spherical\"\n\n    def __init__(self, ds, ordering=(\"r\", \"theta\", \"phi\")):\n        super().__init__(ds, ordering)\n        # Generate\n        self.image_units = {}\n        self.image_units[self.axis_id[\"r\"]] = (\"rad\", \"rad\")\n        self.image_units[self.axis_id[\"theta\"]] = (None, None)\n        self.image_units[self.axis_id[\"phi\"]] = (None, None)\n\n    def setup_fields(self, registry):\n        # return the fields for r, z, theta\n        registry.add_field(\n            (\"index\", \"dx\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        registry.add_field(\n            (\"index\", \"dy\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        registry.add_field(\n            (\"index\", \"dz\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        registry.add_field(\n            (\"index\", \"x\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        registry.add_field(\n            (\"index\", \"y\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        registry.add_field(\n            (\"index\", \"z\"), sampling_type=\"cell\", function=_unknown_coord\n        )\n\n        f1, f2 = _get_coord_fields(self.axis_id[\"r\"])\n        registry.add_field(\n            (\"index\", \"dr\"),\n            sampling_type=\"cell\",\n            function=f1,\n            display_field=False,\n            units=\"code_length\",\n        )\n\n        registry.add_field(\n            (\"index\", \"r\"),\n            sampling_type=\"cell\",\n            function=f2,\n            display_field=False,\n            units=\"code_length\",\n        )\n\n        f1, f2 = _get_coord_fields(self.axis_id[\"theta\"], \"\")\n        registry.add_field(\n            (\"index\", \"dtheta\"),\n            sampling_type=\"cell\",\n            function=f1,\n            display_field=False,\n            units=\"\",\n        )\n\n        registry.add_field(\n            (\"index\", \"theta\"),\n            sampling_type=\"cell\",\n            function=f2,\n            display_field=False,\n            units=\"\",\n        )\n\n        f1, f2 = _get_coord_fields(self.axis_id[\"phi\"], \"\")\n        registry.add_field(\n            (\"index\", \"dphi\"),\n            sampling_type=\"cell\",\n            function=f1,\n            display_field=False,\n            units=\"\",\n        )\n\n        registry.add_field(\n            (\"index\", \"phi\"),\n            sampling_type=\"cell\",\n            function=f2,\n            display_field=False,\n            units=\"\",\n        )\n\n        def _SphericalVolume(field, data):\n            # Here we compute the spherical volume element exactly\n            r = data[\"index\", \"r\"]\n            dr = data[\"index\", \"dr\"]\n            theta = data[\"index\", \"theta\"]\n            dtheta = data[\"index\", \"dtheta\"]\n            vol = ((r + 0.5 * dr) ** 3 - (r - 0.5 * dr) ** 3) / 3.0\n            vol *= np.cos(theta - 0.5 * dtheta) - np.cos(theta + 0.5 * dtheta)\n            vol *= data[\"index\", \"dphi\"]\n            return vol\n\n        registry.add_field(\n            (\"index\", \"cell_volume\"),\n            sampling_type=\"cell\",\n            function=_SphericalVolume,\n            units=\"code_length**3\",\n        )\n        registry.alias((\"index\", \"volume\"), (\"index\", \"cell_volume\"))\n\n        def _path_r(field, data):\n            return data[\"index\", \"dr\"]\n\n        registry.add_field(\n            (\"index\", \"path_element_r\"),\n            sampling_type=\"cell\",\n            function=_path_r,\n            units=\"code_length\",\n        )\n\n        def _path_theta(field, data):\n            # Note: this already assumes cell-centered\n            return data[\"index\", \"r\"] * data[\"index\", \"dtheta\"]\n\n        registry.add_field(\n            (\"index\", \"path_element_theta\"),\n            sampling_type=\"cell\",\n            function=_path_theta,\n            units=\"code_length\",\n        )\n\n        def _path_phi(field, data):\n            # Note: this already assumes cell-centered\n            return (\n                data[\"index\", \"r\"]\n                * data[\"index\", \"dphi\"]\n                * np.sin(data[\"index\", \"theta\"])\n            )\n\n        registry.add_field(\n            (\"index\", \"path_element_phi\"),\n            sampling_type=\"cell\",\n            function=_path_phi,\n            units=\"code_length\",\n        )\n\n    def pixelize(\n        self, dimension, data_source, field, bounds, size, antialias=True, periodic=True\n    ):\n        self.period\n        name = self.axis_name[dimension]\n        if name == \"r\":\n            return self._ortho_pixelize(\n                data_source, field, bounds, size, antialias, dimension, periodic\n            )\n        elif name in (\"theta\", \"phi\"):\n            return self._cyl_pixelize(\n                data_source, field, bounds, size, antialias, dimension\n            )\n        else:\n            raise NotImplementedError\n\n    def pixelize_line(self, field, start_point, end_point, npoints):\n        raise NotImplementedError\n\n    def _ortho_pixelize(\n        self, data_source, field, bounds, size, antialias, dim, periodic\n    ):\n        buff = pixelize_aitoff(\n            data_source[\"py\"],\n            data_source[\"pdy\"],\n            data_source[\"px\"],\n            data_source[\"pdx\"],\n            size,\n            data_source[field],\n            None,\n            None,\n            theta_offset=0,\n            phi_offset=0,\n        ).transpose()\n        return buff\n\n    def _cyl_pixelize(self, data_source, field, bounds, size, antialias, dimension):\n        name = self.axis_name[dimension]\n        buff = np.full((size[1], size[0]), np.nan, dtype=\"f8\")\n        if name == \"theta\":\n            pixelize_cylinder(\n                buff,\n                data_source[\"px\"],\n                data_source[\"pdx\"],\n                data_source[\"py\"],\n                data_source[\"pdy\"],\n                data_source[field],\n                bounds,\n            )\n        elif name == \"phi\":\n            # Note that we feed in buff.T here\n            pixelize_cylinder(\n                buff.T,\n                data_source[\"px\"],\n                data_source[\"pdx\"],\n                data_source[\"py\"],\n                data_source[\"pdy\"],\n                data_source[field],\n                bounds,\n            )\n        else:\n            raise RuntimeError\n        self.sanitize_buffer_fill_values(buff)\n        return buff\n\n    def convert_from_cartesian(self, coord):\n        raise NotImplementedError\n\n    def convert_to_cartesian(self, coord):\n        if isinstance(coord, np.ndarray) and len(coord.shape) > 1:\n            ri = self.axis_id[\"r\"]\n            thetai = self.axis_id[\"theta\"]\n            phii = self.axis_id[\"phi\"]\n            r = coord[:, ri]\n            theta = coord[:, thetai]\n            phi = coord[:, phii]\n            nc = np.zeros_like(coord)\n            # r, theta, phi\n            nc[:, ri] = np.cos(phi) * np.sin(theta) * r\n            nc[:, thetai] = np.sin(phi) * np.sin(theta) * r\n            nc[:, phii] = np.cos(theta) * r\n        else:\n            r, theta, phi = coord\n            nc = (\n                np.cos(phi) * np.sin(theta) * r,\n                np.sin(phi) * np.sin(theta) * r,\n                np.cos(theta) * r,\n            )\n        return nc\n\n    def convert_to_cylindrical(self, coord):\n        raise NotImplementedError\n\n    def convert_from_cylindrical(self, coord):\n        raise NotImplementedError\n\n    def convert_to_spherical(self, coord):\n        raise NotImplementedError\n\n    def convert_from_spherical(self, coord):\n        raise NotImplementedError\n\n    _image_axis_name = None\n\n    @property\n    def image_axis_name(self):\n        if self._image_axis_name is not None:\n            return self._image_axis_name\n        # This is the x and y axes labels that get displayed.  For\n        # non-Cartesian coordinates, we usually want to override these for\n        # Cartesian coordinates, since we transform them.\n        rv = {\n            self.axis_id[\"r\"]: (\"\\\\theta\", \"\\\\phi\"),\n            self.axis_id[\"theta\"]: (\"x / \\\\sin(\\\\theta)\", \"y / \\\\sin(\\\\theta)\"),\n            self.axis_id[\"phi\"]: (\"R\", \"z\"),\n        }\n        for i in list(rv.keys()):\n            rv[self.axis_name[i]] = rv[i]\n            rv[self.axis_name[i].capitalize()] = rv[i]\n        self._image_axis_name = rv\n        return rv\n\n    _x_pairs = ((\"r\", \"theta\"), (\"theta\", \"r\"), (\"phi\", \"r\"))\n    _y_pairs = ((\"r\", \"phi\"), (\"theta\", \"phi\"), (\"phi\", \"theta\"))\n\n    @property\n    def period(self):\n        return self.ds.domain_width\n\n    def sanitize_center(self, center, axis):\n        center, display_center = super().sanitize_center(center, axis)\n        name = self.axis_name[axis]\n        if name == \"r\":\n            display_center = center\n        elif name == \"theta\":\n            display_center = (\n                0.0 * display_center[0],\n                0.0 * display_center[1],\n                0.0 * display_center[2],\n            )\n        elif name == \"phi\":\n            display_center = [\n                self.ds.domain_width[0] / 2.0 + self.ds.domain_left_edge[0],\n                0.0 * display_center[1],\n                0.0 * display_center[2],\n            ]\n            ri = self.axis_id[\"r\"]\n            c = self.ds.domain_width[ri] / 2.0 + self.ds.domain_left_edge[ri]\n            display_center[ri] = c\n            display_center = tuple(display_center)\n        return center, display_center\n\n    def sanitize_width(self, axis, width, depth):\n        name = self.axis_name[axis]\n        if width is not None:\n            width = super().sanitize_width(axis, width, depth)\n        elif name == \"r\":\n            width = [\n                self.ds.domain_width[self.x_axis[\"r\"]],\n                self.ds.domain_width[self.y_axis[\"r\"]],\n            ]\n        elif name == \"theta\":\n            ri = self.axis_id[\"r\"]\n            # Remember, in spherical coordinates when we cut in theta,\n            # we create a conic section\n            width = [2.0 * self.ds.domain_width[ri], 2.0 * self.ds.domain_width[ri]]\n        elif name == \"phi\":\n            ri = self.axis_id[\"r\"]\n            width = [self.ds.domain_right_edge[ri], 2.0 * self.ds.domain_width[ri]]\n        return width\n\n    def _sanity_check(self):\n        \"\"\"This prints out a handful of diagnostics that help verify the\n        dataset is well formed.\"\"\"\n        # We just check a few things here.\n        dd = self.ds.all_data()\n        r0 = self.ds.domain_left_edge[self.axis_id[\"r\"]]\n        r1 = self.ds.domain_right_edge[self.axis_id[\"r\"]]\n        v1 = 4.0 * np.pi / 3.0 * (r1 ** 3 - r0 ** 3)\n        print(f\"Total volume should be 4*pi*r**3 = {v1:0.16e}\")\n        v2 = dd.quantities.total_quantity(\"cell_volume\")\n        print(f\"Actual volume is                   {v2:0.16e}\")\n        print(f\"Relative difference: {np.abs(v2 - v1) / (v2 + v1):0.16e}\")\n", "path": "yt/geometry/coordinates/spherical_coordinates.py"}]}