BUG: fix buffer masking in cartesian geometry

yt-project · Mar 21, 2022 · 014a10f · 014a10f
1 parent ef08989
commit 014a10f
Show file tree

Hide file tree

Showing 5 changed files with 11 additions and 6 deletions.
diff --git a/yt/geometry/coordinates/cartesian_coordinates.py b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -306,7 +306,7 @@ def _ortho_pixelize(
         if hasattr(period, "in_units"):
             period = period.in_units("code_length").d
 
-        buff = np.zeros((size[1], size[0]), dtype="f8")
+        buff = np.full((size[1], size[0]), np.nan, dtype="float64")
         particle_datasets = (ParticleDataset, StreamParticlesDataset)
         is_sph_field = finfo.is_sph_field
 
@@ -555,7 +555,7 @@ def _oblique_pixelize(self, data_source, field, bounds, size, antialias):
         from yt.frontends.ytdata.data_structures import YTSpatialPlotDataset
 
         indices = np.argsort(data_source["pdx"])[::-1].astype(np.int_)
-        buff = np.zeros((size[1], size[0]), dtype="f8")
+        buff = np.full((size[1], size[0]), np.nan, dtype="float64")
         ftype = "index"
         if isinstance(data_source.ds, YTSpatialPlotDataset):
             ftype = "gas"

diff --git a/yt/geometry/coordinates/cylindrical_coordinates.py b/yt/geometry/coordinates/cylindrical_coordinates.py
@@ -125,7 +125,7 @@ def _ortho_pixelize(
         period[1] = self.period[self.y_axis[dim]]
         if hasattr(period, "in_units"):
             period = period.in_units("code_length").d
-        buff = np.zeros(size, dtype="f8")
+        buff = np.full(size, np.nan, dtype="float64")
         pixelize_cartesian(
             buff,
             data_source["px"],

diff --git a/yt/geometry/coordinates/geographic_coordinates.py b/yt/geometry/coordinates/geographic_coordinates.py
@@ -246,7 +246,7 @@ def _ortho_pixelize(
         pdx = data_source["pdx"]
         py = data_source["py"]
         pdy = data_source["pdy"]
-        buff = np.zeros((size[1], size[0]), dtype="f8")
+        buff = np.full((size[1], size[0]), np.nan, dtype="float64")
         pixelize_cartesian(
             buff,
             px,

diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx
@@ -256,6 +256,8 @@ def pixelize_cartesian(np.float64_t[:,:] buff,
                                 # This will reduce artifacts if we ever move to
                                 # compositing instead of replacing bitmaps.
                                 if overlap1 * overlap2 < 1.e-6: continue
+                                # make sure pixel value is not a NaN before incrementing it
+                                if buff[i,j] != buff[i,j]: buff[i,j] = 0.0
                                 buff[i,j] += (dsp * overlap1) * overlap2
                             else:
                                 buff[i,j] = dsp
@@ -506,6 +508,8 @@ def pixelize_off_axis_cartesian(
                        fabs(zsp - cz) * 0.99 > dzsp:
                         continue
                     mask[i, j] += 1
+                    # make sure pixel value is not a NaN before incrementing it
+                    if buff[i,j] != buff[i,j]: buff[i,j] = 0.0
                     buff[i, j] += dsp
     for i in range(buff.shape[0]):
         for j in range(buff.shape[1]):

diff --git a/yt/visualization/tests/test_offaxisprojection.py b/yt/visualization/tests/test_offaxisprojection.py
@@ -3,6 +3,8 @@
 import tempfile
 import unittest
 
+import numpy as np
+
 from yt.testing import (
     assert_equal,
     assert_fname,
@@ -98,5 +100,4 @@ def test_field_cut_off_axis_octree():
         (p3.frb[("gas", "density")] == p4.frb[("gas", "density")]).all(), False
     )
     p4rho = p4.frb[("gas", "density")]
-    assert_equal(p4rho.min() == 0.0, True)  # Lots of zeros
-    assert_equal(p4rho[p4rho > 0.0].min() >= 0.5, True)
+    assert_equal(np.nanmin(p4rho[p4rho > 0.0]) >= 0.5, True)