Merge pull request #4065 from neutrinoceros/manual_bps_4.0.5

REL: manual backports + prep release 4.0.5

doc/source/conf.py

@@ -79,7 +79,7 @@
 # The short X.Y version.
 version = "4.0"
 # The full version, including alpha/beta/rc tags.
-release = "4.0.4"
+release = "4.0.5"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

doc/source/visualizing/callbacks.rst

@@ -563,7 +563,7 @@ To plot only the central particles
 
    ds = yt.load("Enzo_64/DD0043/data0043")
    p = yt.ProjectionPlot(ds, "x", ("gas", "density"), center="m", width=(10, "Mpc"))
-   sp = ds.sphere([0.5, 0.5, 0.5], ds.quan(1, "Mpc"))
+   sp = ds.sphere(p.data_source.center, ds.quan(1, "Mpc"))
    p.annotate_particles((10, "Mpc"), data_source=sp)
    p.save()
 

pyproject.toml

@@ -1,15 +1,9 @@
 [build-system]
 requires = [
-  "setuptools>=19.6",
-  # see https://github.com/numpy/numpy/pull/18389
-  "wheel>=0.36.2",
-
-  # cython version is imposed by that of numpy, see release notes
-  # https://github.com/numpy/numpy/releases/tag/v1.19.2
-  # Cython 3.0 is the next version after 0.29, and a major change,
-  # we forbid it until we can properly test against it
-  "Cython>=0.26.1,<3.0; python_version=='3.6'",
-  "Cython>=0.29.21,<3.0; python_version>='3.7'",
+  "setuptools==59.6.0; python_version=='3.6'",
+  "setuptools==63.4.2; python_version>'3.6'",
+  "wheel==0.37.1",
+  "Cython==0.29.32",
   "oldest-supported-numpy",
 ]
 

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = yt
-version = 4.0.4
+version = 4.0.5
 description = An analysis and visualization toolkit for volumetric data
 long_description = file: README.md
 long_description_content_type = text/markdown

yt/__init__.py

@@ -7,7 +7,7 @@
 * Contribute: https://github.com/yt-project/yt
 
 """
-__version__ = "4.0.4"
+__version__ = "4.0.5"
 
 import yt.units as units
 import yt.utilities.physical_constants as physical_constants

yt/data_objects/construction_data_containers.py

@@ -1672,9 +1672,8 @@ def calculate_flux(self, field_x, field_y, field_z, fluxing_field=None):
 
         Returns
         -------
-        flux : float
-            The summed flux.  Note that it is not currently scaled; this is
-            simply the code-unit area times the fields.
+        flux : YTQuantity
+            The summed flux.
 
         References
         ----------
@@ -1717,7 +1716,10 @@ def _calculate_flux_in_grid(
 
         vc_data = grid.get_vertex_centered_data(vc_fields)
         if fluxing_field is None:
-            ff = np.ones_like(vc_data[self.surface_field], dtype="float64")
+            ff = self.ds.arr(
+                np.ones_like(vc_data[self.surface_field].d, dtype="float64"),
+                "dimensionless",
+            )
         else:
             ff = vc_data[fluxing_field]
         surf_vals = vc_data[self.surface_field]

yt/data_objects/image_array.py

@@ -331,10 +331,8 @@ def write_png(
             sigma_clip = clip_ratio
 
         if sigma_clip is not None:
-            nz = out[:, :, :3][out[:, :, :3].nonzero()]
-            return write_bitmap(
-                out.swapaxes(0, 1), filename, nz.mean() + sigma_clip * nz.std()
-            )
+            clip_value = self._clipping_value(sigma_clip, im=out)
+            return write_bitmap(out.swapaxes(0, 1), filename, clip_value)
         else:
             return write_bitmap(out.swapaxes(0, 1), filename)
 
@@ -442,3 +440,11 @@ def save(self, filename, png=True, hdf5=True, dataset_name=None):
             if not filename.endswith(".h5"):
                 filename = filename + ".h5"
             self.write_hdf5(filename, dataset_name)
+
+    def _clipping_value(self, sigma_clip, im=None):
+        # return the max value to clip with given a sigma_clip value. If im
+        # is None, the current instance is used
+        if im is None:
+            im = self
+        nz = im[:, :, :3][im[:, :, :3].nonzero()]
+        return nz.mean() + sigma_clip * nz.std()

yt/data_objects/tests/test_fluxes.py

@@ -24,6 +24,11 @@ def test_flux_calculation():
     )
     assert_almost_equal(flux.value, 1.0, 12)
     assert_equal(str(flux.units), "cm**2")
+    flux2 = surf.calculate_flux(
+        ("index", "ones"), ("index", "zeros"), ("index", "zeros")
+    )
+    assert_almost_equal(flux2.value, 1.0, 12)
+    assert_equal(str(flux2.units), "cm**2")
 
 
 def test_sampling():

yt/data_objects/tests/test_image_array.py

@@ -134,6 +134,18 @@ def test_write_image(self):
         im_arr.write_image("with_cmap", cmap_name="hot")
         im_arr.write_image("channel_1.png", channel=1)
 
+    def test_clipping_value(self):
+        im_arr = ImageArray(dummy_image(10.0, 4))
+        clip_val1 = im_arr._clipping_value(1)
+        clip_val2 = im_arr._clipping_value(1, im=im_arr)
+        assert clip_val2 == clip_val1
+
+        clip_val3 = im_arr._clipping_value(6)
+        assert clip_val3 > clip_val2
+
+        im_arr[:] = 1.0  # std will be 0, mean will be 1, so clip value will be 1
+        assert im_arr._clipping_value(1) == 1.0
+
     def tearDown(self):
         os.chdir(self.curdir)
         # clean up

yt/data_objects/tests/test_sph_data_objects.py

@@ -82,6 +82,21 @@ def test_region():
                 assert_equal(reg["gas", "density"].shape[0], answer)
 
 
+def test_periodic_region():
+    ds = fake_sph_grid_ds(10.0)
+    coords = [0.7, 1.4, 2.8]
+
+    for x in coords:
+        for y in coords:
+            for z in coords:
+                center = np.array([x, y, z])
+                for n, w in zip((8, 27), (1.0, 2.0)):
+                    le = center - 0.5 * w
+                    re = center + 0.5 * w
+                    box = ds.box(le, re)
+                    assert box["io", "particle_ones"].size == n
+
+
 SPHERE_ANSWERS = {
     ((0, 0, 0), 4): 7,
     ((0, 0, 0), 3): 7,

yt/frontends/enzo_e/data_structures.py

@@ -276,12 +276,15 @@ def _detect_output_fields(self):
             grid = self.grids[0]
             field_list, ptypes = self.io._read_field_names(grid)
             mylog.debug("Grid %s has: %s", grid.id, field_list)
+            sample_pfields = self.io.sample_pfields
         else:
             field_list = None
             ptypes = None
+            sample_pfields = None
         self.field_list = list(self.comm.mpi_bcast(field_list))
         self.dataset.particle_types = list(self.comm.mpi_bcast(ptypes))
         self.dataset.particle_types_raw = self.dataset.particle_types[:]
+        self.io.sample_pfields = self.comm.mpi_bcast(sample_pfields)
 
 
 class EnzoEDataset(Dataset):

yt/funcs.py

@@ -167,7 +167,8 @@ def get_memory_usage(subtract_share=False):
     status_file = f"/proc/{pid}/statm"
     if not os.path.isfile(status_file):
         return -1024
-    line = open(status_file).read()
+    with open(status_file) as fh:
+        line = fh.read()
     size, resident, share, text, library, data, dt = (int(i) for i in line.split())
     if subtract_share:
         resident -= share

yt/geometry/_selection_routines/region_selector.pxi

@@ -145,15 +145,14 @@ cdef class RegionSelector(SelectorObject):
         cdef np.float64_t p
         cdef np.float64_t r2 = radius**2
         cdef np.float64_t dmin = 0
+        cdef np.float64_t d = 0
         for i in range(3):
-            if self.periodicity[i] and self.check_period[i]:
-                p = pos[i] + self.right_edge_shift[i]
-            else:
-                p = pos[i]
-            if p < self.left_edge[i]:
-                dmin += (p - self.left_edge[i])**2
-            elif pos[i] > self.right_edge[i]:
-                dmin += (p - self.right_edge[i])**2
+            if (pos[i]+radius < self.left_edge[i] and \
+                pos[i]-radius >= self.right_edge_shift[i]):
+                d = self.periodic_difference(pos[i], self.left_edge[i], i)
+            elif pos[i]-radius > self.right_edge[i]:
+                d = self.periodic_difference(pos[i], self.right_edge[i], i)
+            dmin += d*d
         return int(dmin <= r2)
 
     @cython.boundscheck(False)

yt/geometry/coordinates/geographic_coordinates.py

@@ -382,6 +382,7 @@ def image_axis_name(self):
 
     @property
     def data_projection(self):
+        # this will control the default projection to use when displaying data
         if self._data_projection is not None:
             return self._data_projection
         dpj = {}
@@ -397,6 +398,7 @@ def data_projection(self):
 
     @property
     def data_transform(self):
+        # this is the coordinate system on which the data is defined (the crs).
         if self._data_transform is not None:
             return self._data_transform
         dtx = {}

yt/utilities/lib/image_samplers.pyx

@@ -267,7 +267,7 @@ cdef class ImageSampler:
                     continue
 
                 for i in range(Nch):
-                    idata.rgba[i] = 0
+                    idata.rgba[i] = self.image[vi, vj, i]
                 for i in range(8):
                     vc.mask[i] = 1
 

yt/visualization/base_plot_types.py

@@ -1,3 +1,4 @@
+import sys
 import warnings
 from io import BytesIO
 
@@ -266,33 +267,8 @@ def _init_image(self, data, cbnorm, cblinthresh, cmap, extent, aspect):
             transform = self.axes.transData
         else:
             transform = self._transform
-        if hasattr(self.axes, "set_extent"):
-            # CartoPy hangs if we do not set_extent before imshow if we are
-            # displaying a small subset of the globe.  What I believe happens is
-            # that the transform for the points on the outside results in
-            # infinities, and then the scipy.spatial cKDTree hangs trying to
-            # identify nearest points.
-            #
-            # Also, set_extent is defined by cartopy, so not all axes will have
-            # it as a method.
-            #
-            # A potential downside is that other images may change, but I believe
-            # the result of imshow is to set_extent *regardless*.  This just
-            # changes the order in which it happens.
-            #
-            # NOTE: This is currently commented out because it breaks in some
-            # instances.  It is left as a historical note because we will
-            # eventually need some form of it.
-            # self.axes.set_extent(extent)
-
-            # possibly related issue (operation order dependency)
-            # https://github.com/SciTools/cartopy/issues/1468
-
-            # in cartopy 0.19 (or 0.20), some intented behaviour changes produced an
-            # incompatibility here where default values for extent would lead to a crash.
-            # A solution is to let the transform object set the image extents internally
-            # see https://github.com/SciTools/cartopy/issues/1955
-            extent = None
+
+        self._validate_axes_extent(extent, transform)
 
         self.image = self.axes.imshow(
             data.to_ndarray(),
@@ -360,6 +336,33 @@ def _init_image(self, data, cbnorm, cblinthresh, cmap, extent, aspect):
             self.cb = self.figure.colorbar(self.image, self.cax)
         self.cax.tick_params(which="both", axis="y", direction="in")
 
+    def _validate_axes_extent(self, extent, transform):
+        # if the axes are cartopy GeoAxes, this checks that the axes extent
+        # is properly set.
+
+        if "cartopy" not in sys.modules:
+            # cartopy isn't already loaded, nothing to do here
+            return
+
+        from cartopy.mpl.geoaxes import GeoAxes
+
+        if isinstance(self.axes, GeoAxes):
+            # some projections have trouble when passing extents at or near the
+            # limits. So we only set_extent when the plot is a subset of the
+            # globe, within the tolerance of the transform.
+
+            # note that `set_extent` here is setting the extent of the axes.
+            # still need to pass the extent arg to imshow in order to
+            # ensure that it is properly scaled. also note that set_extent
+            # expects values in the coordinates of the transform: it will
+            # calculate the coordinates in the projection.
+            global_extent = transform.x_limits + transform.y_limits
+            thresh = transform.threshold
+            if all(
+                abs(extent[ie]) < (abs(global_extent[ie]) - thresh) for ie in range(4)
+            ):
+                self.axes.set_extent(extent, crs=transform)
+
     def _get_best_layout(self):
 
         # Ensure the figure size along the long axis is always equal to _figure_size

yt/visualization/plot_container.py

@@ -1030,7 +1030,7 @@ def _sanitize_units(z, _field):
                     # they can rely on for inspecting units
                     mylog.warning(
                         "%s class doesn't support zmin/zmax"
-                        " as tuples or unyt_quantitiy",
+                        " as tuples or unyt_quantity",
                         self.__class__.__name__,
                     )
                     z = z.value

yt/visualization/tests/test_geo_projections.py

@@ -1,5 +1,6 @@
 import unittest
 
+import numpy as np
 from nose.plugins.attrib import attr
 
 import yt
@@ -62,6 +63,9 @@ class TestGeoProjections(unittest.TestCase):
     @requires_module("cartopy")
     def setUp(self):
         self.ds = fake_amr_ds(geometry="geographic")
+        # switch off the log plot to avoid some unrelated matplotlib issues
+        f = self.ds._get_field_info("stream", "Density")
+        f.take_log = False
 
     @requires_module("cartopy")
     def tearDown(self):
@@ -133,6 +137,12 @@ def test_nondefault_transform(self):
         from yt.utilities.on_demand_imports import _cartopy as cartopy
 
         axis = "altitude"
+        # Note: The Miller transform has an extent of approx. +/- 180 in x,
+        # +/-132 in y (in Miller, x is longitude, y is a factor of latitude).
+        # So by changing the projection in this way, the dataset goes from
+        # covering the whole globe (+/- 90 latitude), to covering part of the
+        # globe (+/-72 latitude). Totally fine for testing that the code runs,
+        # but good to be aware of!
         self.ds.coordinates.data_transform[axis] = "Miller"
         self.slc = yt.SlicePlot(self.ds, axis, ("stream", "Density"), origin="native")
 
@@ -151,6 +161,21 @@ def test_nondefault_transform(self):
                 self.slc.plots[("stream", "Density")].axes.projection, proj_type
             )
 
+    @requires_module("cartopy")
+    def test_extent(self):
+        # checks that the axis extent is narrowed when doing a subselection
+        axis = "altitude"
+        slc = yt.SlicePlot(self.ds, axis, ("stream", "Density"), origin="native")
+        ax = slc.plots[("stream", "Density")].axes
+        full_extent = np.abs(ax.get_extent())
+
+        slc = yt.SlicePlot(
+            self.ds, axis, ("stream", "Density"), origin="native", width=(80.0, 50.0)
+        )
+        ax = slc.plots[("stream", "Density")].axes
+        extent = np.abs(ax.get_extent())
+        assert np.all(extent < full_extent)
+
 
 class TestNonGeoProjections(unittest.TestCase):
     def setUp(self):

yt/visualization/volume_rendering/scene.py

@@ -343,8 +343,10 @@ def save(self, fname=None, sigma_clip=None, render=True):
             ax = fig.add_axes([0, 0, 1, 1])
             ax.set_axis_off()
             out = self._last_render
-            nz = out[:, :, :3][out[:, :, :3].nonzero()]
-            max_val = nz.mean() + sigma_clip * nz.std()
+            if sigma_clip is not None:
+                max_val = out._clipping_value(sigma_clip)
+            else:
+                max_val = out[:, :, :3].max()
             alpha = 255 * out[:, :, 3].astype("uint8")
             out = np.clip(out[:, :, :3] / max_val, 0.0, 1.0) * 255
             out = np.concatenate([out.astype("uint8"), alpha[..., None]], axis=-1)
@@ -472,11 +474,9 @@ def _show_mpl(self, im, sigma_clip=None, dpi=100):
         ax.set_position([0, 0, 1, 1])
 
         if sigma_clip is not None:
-            nz = im[im > 0.0]
-            nim = im / (nz.mean() + sigma_clip * np.std(nz))
+            nim = im / im._clipping_value(sigma_clip)
             nim[nim > 1.0] = 1.0
             nim[nim < 0.0] = 0.0
-            del nz
         else:
             nim = im
         axim = ax.imshow(nim[:, :, :3] / nim[:, :, :3].max(), interpolation="bilinear")