WIP: add pixel-averaging option for SPH projections

nose_ignores.txt

@@ -48,5 +48,6 @@
 --ignore-file=test_disks\.py
 --ignore-file=test_offaxisprojection_pytestonly\.py
 --ignore-file=test_sph_pixelization_pytestonly\.py
+--ignore-file=test_off_axis_SPH\.py
 --ignore-file=test_time_series\.py
 --ignore-file=test_cf_radial_pytest\.py

tests/tests.yaml

@@ -231,5 +231,6 @@ other_tests:
      - "--ignore-file=test_offaxisprojection_pytestonly\\.py"
      - "--ignore-file=test_sph_pixelization_pytestonly\\.py"
      - "--ignore-file=test_cf_radial_pytest\\.py"
+     - "--ignore-file=test_off_axis_SPH\\.py"
   cookbook:
      - 'doc/source/cookbook/tests/test_cookbook.py'

yt/data_objects/selection_objects/data_selection_objects.py

@@ -4,6 +4,7 @@
 import uuid
 from collections import defaultdict
 from contextlib import contextmanager
+from typing import Literal
 
 import numpy as np
 from more_itertools import always_iterable
@@ -567,7 +568,15 @@ def _get_pw(self, fields, center, width, origin, plot_type):
         pw._setup_plots()
         return pw
 
-    def to_frb(self, width, resolution, center=None, height=None, periodic=False):
+    def to_frb(
+        self,
+        width,
+        resolution,
+        center=None,
+        height=None,
+        periodic=False,
+        pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave",
+    ):
         r"""This function returns a FixedResolutionBuffer generated from this
         object.
 
@@ -595,7 +604,15 @@ def to_frb(self, width, resolution, center=None, height=None, periodic=False):
         periodic : bool
             Should the returned Fixed Resolution Buffer be periodic?  (default:
             False).
+        pixelmeaning:
+            "pixelav": a pixel represents an average surface density or
+            surface-density-weighted average across a pixel.
+
+            "pencilbeam": a pixel represents a column density or
+            column-density-weighted average integrated over a pencil
+            beam through the pixel center.
 
+            Only applies to projections of SPH datasets.
         Returns
         -------
         frb : :class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer`
@@ -659,7 +676,9 @@ def to_frb(self, width, resolution, center=None, height=None, periodic=False):
             center[yax] - height * 0.5,
             center[yax] + height * 0.5,
         )
-        frb = FixedResolutionBuffer(self, bounds, resolution, periodic=periodic)
+        frb = FixedResolutionBuffer(
+            self, bounds, resolution, periodic=periodic, pixelmeaning=pixelmeaning
+        )
         return frb
 
 

yt/data_objects/selection_objects/slices.py

@@ -1,3 +1,5 @@
+from typing import Literal
+
 import numpy as np
 
 from yt.data_objects.selection_objects.data_selection_objects import (
@@ -313,7 +315,14 @@ def to_pw(self, fields=None, center="center", width=None, axes_unit=None):
         pw._setup_plots()
         return pw
 
-    def to_frb(self, width, resolution, height=None, periodic=False):
+    def to_frb(
+        self,
+        width,
+        resolution,
+        height=None,
+        periodic=False,
+        pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave",
+    ):
         r"""This function returns a FixedResolutionBuffer generated from this
         object.
 
@@ -338,6 +347,8 @@ def to_frb(self, width, resolution, height=None, periodic=False):
         periodic : boolean
             This can be true or false, and governs whether the pixelization
             will span the domain boundaries.
+        pixelmeaning: ignored
+            argument meant for SPH projections
 
         Returns
         -------

yt/geometry/coordinates/cartesian_coordinates.py

@@ -1,3 +1,5 @@
+from typing import Literal
+
 import numpy as np
 
 from yt.data_objects.index_subobjects.unstructured_mesh import SemiStructuredMesh
@@ -13,7 +15,8 @@
     pixelize_element_mesh_line,
     pixelize_off_axis_cartesian,
     pixelize_sph_kernel_cutting,
-    pixelize_sph_kernel_projection,
+    pixelize_sph_kernel_projection_pencilbeam,
+    pixelize_sph_kernel_projection_pixelave,
     pixelize_sph_kernel_slice,
 )
 from yt.utilities.math_utils import compute_stddev_image
@@ -170,13 +173,24 @@ def pixelize(
         size,
         antialias=True,
         periodic=True,
+        pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave",
         *,
         return_mask=False,
     ):
         """
         Method for pixelizing datasets in preparation for
         two-dimensional image plots. Relies on several sampling
         routines written in cython
+
+        pixelmeaning:
+            "pixelave": a pixel represents an average surface density or
+            surface-density-weighted average across a pixel.
+
+            "pencilbeam": a pixel represents a column density or
+            column-density-weighted average integrated over a pencil
+            beam through the pixel center.
+
+            Only applies to SPH datasets.
         """
         index = data_source.ds.index
         if hasattr(index, "meshes") and not isinstance(
@@ -236,11 +250,18 @@ def pixelize(
 
         elif self.axis_id.get(dimension, dimension) is not None:
             buff, mask = self._ortho_pixelize(
-                data_source, field, bounds, size, antialias, dimension, periodic
+                data_source,
+                field,
+                bounds,
+                size,
+                antialias,
+                dimension,
+                periodic,
+                pixelmeaning,
             )
         else:
             buff, mask = self._oblique_pixelize(
-                data_source, field, bounds, size, antialias
+                data_source, field, bounds, size, antialias, pixelmeaning
             )
 
         if return_mask:
@@ -314,7 +335,7 @@ def pixelize_line(self, field, start_point, end_point, npoints):
         return arc_length, plot_values
 
     def _ortho_pixelize(
-        self, data_source, field, bounds, size, antialias, dim, periodic
+        self, data_source, field, bounds, size, antialias, dim, periodic, pixelmeaning
     ):
         from yt.data_objects.construction_data_containers import YTParticleProj
         from yt.data_objects.selection_objects.slices import YTSlice
@@ -396,6 +417,20 @@ def _ortho_pixelize(
                 kernel_name = "cubic"
 
             if isinstance(data_source, YTParticleProj):  # projection
+                # pick out projection function. _pencilbeam does have
+                # one extra optional parameter (minimum number of
+                # subsampling pixels), but we just use the default value
+                # for now.
+                if pixelmeaning == "pencilbeam":
+                    pixelize_sph_kernel_projection = (
+                        pixelize_sph_kernel_projection_pencilbeam
+                    )
+                elif pixelmeaning == "pixelave":
+                    pixelize_sph_kernel_projection = (
+                        pixelize_sph_kernel_projection_pixelave
+                    )
+                else:
+                    raise NotImplementedError(f"No pixelmeaning option {pixelmeaning}")
                 weight = data_source.weight_field
                 moment = data_source.moment
                 le, re = data_source.data_source.get_bbox()
@@ -667,7 +702,9 @@ def _ortho_pixelize(
         assert mask is None or mask.dtype == bool
         return buff, mask
 
-    def _oblique_pixelize(self, data_source, field, bounds, size, antialias):
+    def _oblique_pixelize(
+        self, data_source, field, bounds, size, antialias, pixelmeaning
+    ):
         from yt.data_objects.selection_objects.slices import YTCuttingPlane
         from yt.frontends.sph.data_structures import ParticleDataset
         from yt.frontends.stream.data_structures import StreamParticlesDataset

yt/geometry/coordinates/coordinate_handler.py

@@ -156,6 +156,7 @@ def pixelize(
         size,
         antialias=True,
         periodic=True,
+        pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave",
         *,
         return_mask: Literal[False],
     ) -> "np.ndarray[Any, np.dtype[np.float64]]": ...
@@ -170,6 +171,7 @@ def pixelize(
         size,
         antialias=True,
         periodic=True,
+        pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave",
         *,
         return_mask: Literal[True],
     ) -> tuple[
@@ -186,6 +188,7 @@ def pixelize(
         size,
         antialias=True,
         periodic=True,
+        pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave",
         *,
         return_mask=False,
     ):

yt/geometry/coordinates/cylindrical_coordinates.py

@@ -1,4 +1,5 @@
 from functools import cached_property
+from typing import Literal
 
 import numpy as np
 
@@ -87,9 +88,16 @@ def pixelize(
         size,
         antialias=True,
         periodic=False,
+        pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave",
         *,
         return_mask=False,
     ):
+        """
+        Parameters
+        ----------
+        pixelmeaning: ignored, argument meant for cartesian SPH data
+
+        """
         # Note that above, we set periodic by default to be *false*.  This is
         # because our pixelizers, at present, do not handle periodicity
         # correctly, and if you change the "width" of a cylindrical plot, it

yt/geometry/coordinates/geographic_coordinates.py

@@ -1,3 +1,5 @@
+from typing import Literal
+
 import numpy as np
 import unyt
 
@@ -228,9 +230,16 @@ def pixelize(
         size,
         antialias=True,
         periodic=True,
+        pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave",
         *,
         return_mask=False,
     ):
+        """
+        Parameters
+        ----------
+        pixelmeaning: ignored, argument meant for cartesian SPH data
+
+        """
         if self.axis_name[dimension] in ("latitude", "longitude"):
             buff, mask = self._cyl_pixelize(
                 data_source, field, bounds, size, antialias, dimension

yt/geometry/coordinates/spherical_coordinates.py

@@ -1,4 +1,5 @@
 from functools import cached_property
+from typing import Literal
 
 import numpy as np
 
@@ -90,9 +91,16 @@ def pixelize(
         size,
         antialias=True,
         periodic=True,
+        pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave",
         *,
         return_mask=False,
     ):
+        """
+        Parameters
+        ----------
+        pixelmeaning: ignored, argument meant for cartesian SPH data
+
+        """
         self.period
         name = self.axis_name[dimension]
         if name == "r":

yt/geometry/coordinates/tests/test_sph_pixelization.py

@@ -66,6 +66,7 @@ def test_sph_projection_basic1():
         height=(2.5, "cm"),
         center=np.array([1.5, 1.5, 1.5]),
         periodic=False,
+        pixelmeaning="pencilbeam",
     )
     out = frb.get_image(("gas", "density"))
 
@@ -96,6 +97,7 @@ def test_sph_projection_basic2():
         height=(2.5, "cm"),
         center=np.array([1.375, 1.375, 1.5]),
         periodic=False,
+        pixelmeaning="pencilbeam",
     )
     out = frb.get_image(("gas", "density"))
 
@@ -163,6 +165,7 @@ def getdata_test_gridproj2():
             height=(2.5, "cm"),
             center=np.array([1.5, 1.5, 1.5]),
             periodic=False,
+            pixelmeaning="pencilbeam",
         )
         out = frb.get_image(("gas", "density"))
         outlist.append(out)
@@ -205,6 +208,7 @@ def test_sph_gridproj_reseffect2():
             height=(3.0 - 2.0 * margin, "cm"),
             center=np.array([1.5, 1.5, 1.5]),
             periodic=False,
+            pixelmeaning="pencilbeam",
         )
         out = frb.get_image(("gas", "density"))
         imgs[rl] = out