BUG: inconsistency in cartesian buffer masking

[neutrinoceros]

Bug report

Bug summary

There is a discrepancy in what fill value is used in buffers depending on the pixelisation routine. This is apparent e.g. in CylindricalCoordinateHandler:

when cartesian pixelisation is used:

yt/yt/geometry/coordinates/cylindrical_coordinates.py

Line 128 in 0d84602

buff = np.zeros(size, dtype="f8")

and when cylindrical pixelisation is used:

yt/yt/geometry/coordinates/cylindrical_coordinates.py

Line 144 in 0d84602

buff = np.full((size[1], size[0]), np.nan, dtype="f8")

See #2480 and #2517 for historical context.
In short, we switched to NaN as a fill value to resolve issues with setting background colors.

This was not yet generalised to cartesian pixelisation because usually, every single pixel is filled with data values, leaving no room for a background color.
This is however not true in the case of the very example used in our documentation to demonstrate how to setup a background color, where we deliberately zoom out of the domain bounding box.

This example happened to work until yt 4.0.2 only because the norm used there is log, but it breaks if one switches to linear scaling

import yt

ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
normal="z"
slc = yt.SlicePlot(ds, normal, ("gas", "density"), width=(1.5, "Mpc"))
slc.set_background_color(("gas", "density"), color="C0")
slc.save("/tmp/background_log")
slc.set_log(("gas", "density"), False)
slc.save("/tmp/background_lin")

On the main branch, as well as on the 4.0.x branch, having zeros in the buffer implies that yt will switch to a symlog norm, which breaks even the default case. This behaviour is intentional and happened in #3793, but it should be apparent by now that this is incompatible with using 0 as a fill value

So ! I made a simple attempt using the following patch

diff --git a/yt/geometry/coordinates/cartesian_coordinates.py b/yt/geometry/coordinates/cartesian_coordinates.py
index 5930e44a2..c83bcb9b2 100644
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -306,7 +306,7 @@ class CartesianCoordinateHandler(CoordinateHandler):
         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 @@ class CartesianCoordinateHandler(CoordinateHandler):
         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
index 4f34d62d2..a865ff1e0 100644
--- a/yt/geometry/coordinates/cylindrical_coordinates.py
+++ b/yt/geometry/coordinates/cylindrical_coordinates.py
@@ -125,7 +125,7 @@ class CylindricalCoordinateHandler(CoordinateHandler):
         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
index 38c948835..ea5371393 100644
--- a/yt/geometry/coordinates/geographic_coordinates.py
+++ b/yt/geometry/coordinates/geographic_coordinates.py
@@ -246,7 +246,7 @@ class GeographicCoordinateHandler(CoordinateHandler):
         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
index c52d22ab3..3a0309aab 100644
--- 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
index a7af54e77..cb64c2db2 100644
--- a/yt/visualization/tests/test_offaxisprojection.py
+++ b/yt/visualization/tests/test_offaxisprojection.py
@@ -3,6 +3,8 @@ import shutil
 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)

which indeed fixed the original problem with background colors, however this also has unexpected and undesirable side effects in other image comparison tests, so I want to have a dedicated PR for this problem, instead of jamming it in #3849.

My conclusion for now is this needs to be adressed for yt 4.0.3

[matthewturk]

I'm not sure I agree that this is necessary for 4.0.3 since the scaling changes in #3849 would not go in 4.0.3, but be reserved for 4.1.

That being said, I do agree we should fill with NaNs when we create the buffer.

A more robust, long-term approach would be to address this via having a centralized buffer creation step, which would have the side effect of making it possible to overplot multiple different types of data (for instance, mixed cylindrical/cartesian data), but that is out of scope.

[neutrinoceros]

I'm not sure I agree that this is necessary for 4.0.3 since the scaling changes in #3849 would not go in 4.0.3, but be reserved for 4.1.

#3849 is indeed not meant for a bugfix release, but the breaking change was in #3793, which was already backported (but not shipped yet ! there's still the option of reverting it on the backport branch if needed).

edit: to be clear there is no version of yt in which this works with linear scaling, but the main branch as well as the backport branch are arguably even more broken since #3793