Merge pull request #3844 from neutrinoceros/manual_bp_3820

Manual backport #3820 to yt-4.0.x (BUG: fix plot bounds for slices in cylindrical coordinates)

yt/geometry/coordinates/coordinate_handler.py

@@ -1,5 +1,6 @@
 import weakref
 from numbers import Number
+from typing import Tuple
 
 import numpy as np
 from packaging.version import Version
@@ -309,3 +310,51 @@ def cylindrical_to_cartesian(coord, center=(0, 0, 0)):
     c2[..., 1] = np.sin(coord[..., 0]) * coord[..., 1] + center[1]
     c2[..., 2] = coord[..., 2]
     return c2
+
+
+def _get_polar_bounds(self: CoordinateHandler, axes: Tuple[str, str]):
+    # a small helper function that is needed by two unrelated classes
+    ri = self.axis_id[axes[0]]
+    pi = self.axis_id[axes[1]]
+    rmin = self.ds.domain_left_edge[ri]
+    rmax = self.ds.domain_right_edge[ri]
+    phimin = self.ds.domain_left_edge[pi]
+    phimax = self.ds.domain_right_edge[pi]
+    corners = [
+        (rmin, phimin),
+        (rmin, phimax),
+        (rmax, phimin),
+        (rmax, phimax),
+    ]
+
+    def to_polar_plane(r, phi):
+        x = r * np.cos(phi)
+        y = r * np.sin(phi)
+        return x, y
+
+    conic_corner_coords = [to_polar_plane(*corner) for corner in corners]
+
+    phimin = phimin.d
+    phimax = phimax.d
+
+    if phimin <= np.pi <= phimax:
+        xxmin = -rmax
+    else:
+        xxmin = min(xx for xx, yy in conic_corner_coords)
+
+    if phimin <= 0 <= phimax:
+        xxmax = rmax
+    else:
+        xxmax = max(xx for xx, yy in conic_corner_coords)
+
+    if phimin <= 3 * np.pi / 2 <= phimax:
+        yymin = -rmax
+    else:
+        yymin = min(yy for xx, yy in conic_corner_coords)
+
+    if phimin <= np.pi / 2 <= phimax:
+        yymax = rmax
+    else:
+        yymax = max(yy for xx, yy in conic_corner_coords)
+
+    return xxmin, xxmax, yymin, yymax

yt/geometry/coordinates/cylindrical_coordinates.py

@@ -1,10 +1,18 @@
+import sys
+
 import numpy as np
 
 from yt.utilities.lib.pixelization_routines import pixelize_cartesian, pixelize_cylinder
 
+if sys.version_info >= (3, 8):
+    from functools import cached_property
+else:
+    from yt._maintenance.backports import cached_property
+
 from .coordinate_handler import (
     CoordinateHandler,
     _get_coord_fields,
+    _get_polar_bounds,
     _unknown_coord,
     cartesian_to_cylindrical,
     cylindrical_to_cartesian,
@@ -160,6 +168,13 @@ def pixelize(
                 data_source, field, bounds, size, antialias, dimension, periodic
             )
         elif ax_name == "z":
+            # This is admittedly a very hacky way to resolve a bug
+            # it's very likely that the *right* fix would have to
+            # be applied upstream of this function, *but* this case
+            # has never worked properly so maybe it's still preferable to
+            # not having a solution ?
+            # see https://github.com/yt-project/yt/pull/3533
+            bounds = (*bounds[2:4], *bounds[:2])
             return self._cyl_pixelize(data_source, field, bounds, size, antialias)
         else:
             # Pixelizing along a cylindrical surface is a bit tricky
@@ -215,7 +230,7 @@ def image_axis_name(self):
         # non-Cartesian coordinates, we usually want to override these for
         # Cartesian coordinates, since we transform them.
         rv = {
-            self.axis_id["r"]: ("theta", "z"),
+            self.axis_id["r"]: ("\\theta", "z"),
             self.axis_id["z"]: ("x", "y"),
             self.axis_id["theta"]: ("r", "z"),
         }
@@ -247,6 +262,10 @@ def convert_from_spherical(self, coord):
     def period(self):
         return np.array([0.0, 0.0, 2.0 * np.pi])
 
+    @cached_property
+    def _polar_bounds(self):
+        return _get_polar_bounds(self, axes=("r", "theta"))
+
     def sanitize_center(self, center, axis):
         center, display_center = super().sanitize_center(center, axis)
         display_center = [
@@ -265,6 +284,11 @@ def sanitize_center(self, center, axis):
             # use existing center value
             for idx in (r_ax, z_ax):
                 display_center[idx] = center[idx]
+        elif ax_name == "z":
+            xxmin, xxmax, yymin, yymax = self._polar_bounds
+            xc = (xxmin + xxmax) / 2
+            yc = (yymin + yymax) / 2
+            display_center = (xc, yc, 0 * xc)
         return center, display_center
 
     def sanitize_width(self, axis, width, depth):
@@ -277,12 +301,12 @@ def sanitize_width(self, axis, width, depth):
         # Note: regardless of axes, these are set up to give consistent plots
         # when plotted, which is not strictly a "right hand rule" for axes.
         elif name == "r":  # soup can label
-            width = [2.0 * np.pi * self.ds.domain_width.uq, self.ds.domain_width[z_ax]]
+            width = [self.ds.domain_width[theta_ax], self.ds.domain_width[z_ax]]
         elif name == "theta":
             width = [self.ds.domain_width[r_ax], self.ds.domain_width[z_ax]]
         elif name == "z":
-            width = [
-                2.0 * self.ds.domain_right_edge[r_ax],
-                2.0 * self.ds.domain_right_edge[r_ax],
-            ]
+            xxmin, xxmax, yymin, yymax = self._polar_bounds
+            xw = xxmax - xxmin
+            yw = yymax - yymin
+            width = [xw, yw]
         return width

yt/geometry/coordinates/spherical_coordinates.py

@@ -4,7 +4,12 @@
 
 from yt.utilities.lib.pixelization_routines import pixelize_aitoff, pixelize_cylinder
 
-from .coordinate_handler import CoordinateHandler, _get_coord_fields, _unknown_coord
+from .coordinate_handler import (
+    CoordinateHandler,
+    _get_coord_fields,
+    _get_polar_bounds,
+    _unknown_coord,
+)
 
 if sys.version_info >= (3, 8):
     from functools import cached_property
@@ -343,50 +348,7 @@ def to_poloidal_plane(r, theta):
 
     @cached_property
     def _conic_bounds(self):
-        ri = self.axis_id["r"]
-        pi = self.axis_id["phi"]
-        rmin = self.ds.domain_left_edge[ri]
-        rmax = self.ds.domain_right_edge[ri]
-        phimin = self.ds.domain_left_edge[pi]
-        phimax = self.ds.domain_right_edge[pi]
-        corners = [
-            (rmin, phimin),
-            (rmin, phimax),
-            (rmax, phimin),
-            (rmax, phimax),
-        ]
-
-        def to_conic_plane(r, phi):
-            x = r * np.cos(phi)
-            y = r * np.sin(phi)
-            return x, y
-
-        conic_corner_coords = [to_conic_plane(*corner) for corner in corners]
-
-        phimin = phimin.d
-        phimax = phimax.d
-
-        if phimin <= np.pi <= phimax:
-            xxmin = -rmax
-        else:
-            xxmin = min(xx for xx, yy in conic_corner_coords)
-
-        if phimin <= 0 <= phimax:
-            xxmax = rmax
-        else:
-            xxmax = max(xx for xx, yy in conic_corner_coords)
-
-        if phimin <= 3 * np.pi / 2 <= phimax:
-            yymin = -rmax
-        else:
-            yymin = min(yy for xx, yy in conic_corner_coords)
-
-        if phimin <= np.pi / 2 <= phimax:
-            yymax = rmax
-        else:
-            yymax = max(yy for xx, yy in conic_corner_coords)
-
-        return xxmin, xxmax, yymin, yymax
+        return _get_polar_bounds(self, axes=("r", "phi"))
 
     @cached_property
     def _aitoff_bounds(self):

yt/visualization/plot_window.py

@@ -1643,6 +1643,7 @@ def __init__(
             "cylindrical",
             "geographic",
             "internal_geographic",
+            "polar",
         ):
             mylog.info("Setting origin='native' for %s geometry.", ds.geometry)
             origin = "native"
@@ -1856,6 +1857,7 @@ def __init__(
             "cylindrical",
             "geographic",
             "internal_geographic",
+            "polar",
         ):
             mylog.info("Setting origin='native' for %s geometry.", ds.geometry)
             origin = "native"