1import weakref
2from functools import wraps
3
4import numpy as np
5
6from yt.data_objects.image_array import ImageArray
7from yt.frontends.ytdata.utilities import save_as_dataset
8from yt.funcs import get_output_filename, iter_fields, mylog
9from yt.loaders import load_uniform_grid
10from yt.utilities.lib.api import CICDeposit_2, add_points_to_greyscale_image
11from yt.utilities.lib.pixelization_routines import pixelize_cylinder
12from yt.utilities.on_demand_imports import _h5py as h5py
13
14from .fixed_resolution_filters import apply_filter, filter_registry
15from .volume_rendering.api import off_axis_projection
16
17
18class FixedResolutionBuffer:
19    r"""
20    FixedResolutionBuffer(data_source, bounds, buff_size, antialias = True)
21
22    This accepts a 2D data object, such as a Projection or Slice, and
23    implements a protocol for generating a pixelized, fixed-resolution
24    image buffer.
25
26    yt stores 2D AMR data internally as a set of 2D coordinates and the
27    half-width of individual pixels.  Converting this to an image buffer
28    requires a deposition step, where individual variable-resolution pixels
29    are deposited into a buffer of some resolution, to create an image.
30    This object is an interface to that pixelization step: it can deposit
31    multiple fields.  It acts as a standard YTDataContainer object, such that
32    dict-style access returns an image of a given field.
33
34    Parameters
35    ----------
36    data_source : :class:`yt.data_objects.construction_data_containers.YTQuadTreeProj`
37                   or :class:`yt.data_objects.selection_data_containers.YTSlice`
38        This is the source to be pixelized, which can be a projection, slice or
39        cutting plane.
40    bounds : sequence of floats
41        Bounds are the min and max in the image plane that we want our
42        image to cover.  It's in the order of (xmin, xmax, ymin, ymax),
43        where the coordinates are all in the appropriate code units.
44    buff_size : sequence of ints
45        The size of the image to generate.
46    antialias : boolean
47        This can be true or false.  It determines whether or not sub-pixel
48        rendering is used during data deposition.
49    periodic : boolean
50        This can be true or false, and governs whether the pixelization
51        will span the domain boundaries.
52
53    Examples
54    --------
55    To make a projection and then several images, you can generate a
56    single FRB and then access multiple fields:
57
58    >>> proj = ds.proj(0, ("gas", "density"))
59    >>> frb1 = FixedResolutionBuffer(proj, (0.2, 0.3, 0.4, 0.5), (1024, 1024))
60    >>> print(frb1[("gas", "density")].max())
61    1.0914e-9 g/cm**3
62    >>> print(frb1[("gas", "temperature")].max())
63    104923.1 K
64    """
65    _exclude_fields = (
66        "pz",
67        "pdz",
68        "dx",
69        "x",
70        "y",
71        "z",
72        "r",
73        "dr",
74        "phi",
75        "dphi",
76        "theta",
77        "dtheta",
78        ("index", "dx"),
79        ("index", "x"),
80        ("index", "y"),
81        ("index", "z"),
82        ("index", "r"),
83        ("index", "dr"),
84        ("index", "phi"),
85        ("index", "dphi"),
86        ("index", "theta"),
87        ("index", "dtheta"),
88    )
89
90    def __init__(self, data_source, bounds, buff_size, antialias=True, periodic=False):
91        self.data_source = data_source
92        self.ds = data_source.ds
93        self.bounds = bounds
94        self.buff_size = (int(buff_size[0]), int(buff_size[1]))
95        self.antialias = antialias
96        self.data = {}
97        self._filters = []
98        self.axis = data_source.axis
99        self.periodic = periodic
100
101        ds = getattr(data_source, "ds", None)
102        if ds is not None:
103            ds.plots.append(weakref.proxy(self))
104
105        # Handle periodicity, just in case
106        if self.data_source.axis < 3:
107            DLE = self.ds.domain_left_edge
108            DRE = self.ds.domain_right_edge
109            DD = float(self.periodic) * (DRE - DLE)
110            axis = self.data_source.axis
111            xax = self.ds.coordinates.x_axis[axis]
112            yax = self.ds.coordinates.y_axis[axis]
113            self._period = (DD[xax], DD[yax])
114            self._edges = ((DLE[xax], DRE[xax]), (DLE[yax], DRE[yax]))
115
116        self.setup_filters()
117
118    def keys(self):
119        return self.data.keys()
120
121    def __delitem__(self, item):
122        del self.data[item]
123
124    def __getitem__(self, item):
125        if item in self.data:
126            return self.data[item]
127        mylog.info(
128            "Making a fixed resolution buffer of (%s) %d by %d",
129            item,
130            self.buff_size[0],
131            self.buff_size[1],
132        )
133        bounds = []
134        for b in self.bounds:
135            if hasattr(b, "in_units"):
136                b = float(b.in_units("code_length"))
137            bounds.append(b)
138
139        buff = self.ds.coordinates.pixelize(
140            self.data_source.axis,
141            self.data_source,
142            item,
143            bounds,
144            self.buff_size,
145            int(self.antialias),
146        )
147
148        for name, (args, kwargs) in self._filters:
149            buff = filter_registry[name](*args[1:], **kwargs).apply(buff)
150
151        # FIXME FIXME FIXME we shouldn't need to do this for projections
152        # but that will require fixing data object access for particle
153        # projections
154        try:
155            if hasattr(item, "name"):
156                it = item.name
157            else:
158                it = item
159            units = self.data_source._projected_units[it]
160        except (KeyError, AttributeError):
161            units = self.data_source[item].units
162
163        ia = ImageArray(buff, units=units, info=self._get_info(item))
164        self.data[item] = ia
165        return self.data[item]
166
167    def __setitem__(self, item, val):
168        self.data[item] = val
169
170    def _get_data_source_fields(self):
171        exclude = self.data_source._key_fields + list(self._exclude_fields)
172        fields = getattr(self.data_source, "fields", [])
173        fields += getattr(self.data_source, "field_data", {}).keys()
174        for f in fields:
175            if f not in exclude and f[0] not in self.data_source.ds.particle_types:
176                self[f]
177
178    def _get_info(self, item):
179        info = {}
180        ftype, fname = field = self.data_source._determine_fields(item)[0]
181        finfo = self.data_source.ds._get_field_info(*field)
182        info["data_source"] = self.data_source.__str__()
183        info["axis"] = self.data_source.axis
184        info["field"] = str(item)
185        info["xlim"] = self.bounds[:2]
186        info["ylim"] = self.bounds[2:]
187        info["length_unit"] = self.data_source.ds.length_unit
188        info["length_to_cm"] = info["length_unit"].in_cgs().to_ndarray()
189        info["center"] = self.data_source.center
190
191        try:
192            info["coord"] = self.data_source.coord
193        except AttributeError:
194            pass
195
196        try:
197            info["weight_field"] = self.data_source.weight_field
198        except AttributeError:
199            pass
200
201        info["label"] = finfo.get_latex_display_name()
202
203        return info
204
205    def convert_to_pixel(self, coords):
206        r"""This function converts coordinates in code-space to pixel-space.
207
208        Parameters
209        ----------
210        coords : sequence of array_like
211            This is (x_coord, y_coord).  Because of the way the math is done,
212            these can both be arrays.
213
214        Returns
215        -------
216        output : sequence of array_like
217            This returns px_coord, py_coord
218
219        """
220        dpx = (self.bounds[1] - self.bounds[0]) / self.buff_size[0]
221        dpy = (self.bounds[3] - self.bounds[2]) / self.buff_size[1]
222        px = (coords[0] - self.bounds[0]) / dpx
223        py = (coords[1] - self.bounds[2]) / dpy
224        return (px, py)
225
226    def convert_distance_x(self, distance):
227        r"""This function converts code-space distance into pixel-space
228        distance in the x-coordinate.
229
230        Parameters
231        ----------
232        distance : array_like
233            This is x-distance in code-space you would like to convert.
234
235        Returns
236        -------
237        output : array_like
238            The return value is the distance in the y-pixel coordinates.
239
240        """
241        dpx = (self.bounds[1] - self.bounds[0]) / self.buff_size[0]
242        return distance / dpx
243
244    def convert_distance_y(self, distance):
245        r"""This function converts code-space distance into pixel-space
246        distance in the y-coordinate.
247
248        Parameters
249        ----------
250        distance : array_like
251            This is y-distance in code-space you would like to convert.
252
253        Returns
254        -------
255        output : array_like
256            The return value is the distance in the x-pixel coordinates.
257
258        """
259        dpy = (self.bounds[3] - self.bounds[2]) / self.buff_size[1]
260        return distance / dpy
261
262    def set_unit(self, field, unit, equivalency=None, equivalency_kwargs=None):
263        """Sets a new unit for the requested field
264
265        parameters
266        ----------
267        field : string or field tuple
268           The name of the field that is to be changed.
269
270        unit : string or Unit object
271           The name of the new unit.
272
273        equivalency : string, optional
274           If set, the equivalency to use to convert the current units to
275           the new requested unit. If None, the unit conversion will be done
276           without an equivalency
277
278        equivalency_kwargs : string, optional
279           Keyword arguments to be passed to the equivalency. Only used if
280           ``equivalency`` is set.
281        """
282        if equivalency_kwargs is None:
283            equivalency_kwargs = {}
284        field = self.data_source._determine_fields(field)[0]
285        if equivalency is None:
286            self[field].convert_to_units(unit)
287        else:
288            equiv_array = self[field].to_equivalent(
289                unit, equivalency, **equivalency_kwargs
290            )
291            # equiv_array isn't necessarily an ImageArray. This is an issue
292            # inherent to the way the unit system handles YTArray
293            # subclasses and I don't see how to modify the unit system to
294            # fix this. Instead, we paper over this issue and hard code
295            # that equiv_array is an ImageArray
296            self[field] = ImageArray(
297                equiv_array,
298                equiv_array.units,
299                equiv_array.units.registry,
300                self[field].info,
301            )
302
303    def export_hdf5(self, filename, fields=None):
304        r"""Export a set of fields to a set of HDF5 datasets.
305
306        This function will export any number of fields into datasets in a new
307        HDF5 file.
308
309        Parameters
310        ----------
311        filename : string
312            This file will be opened in "append" mode.
313        fields : list of strings
314            These fields will be pixelized and output.
315        """
316        if fields is None:
317            fields = list(self.data.keys())
318        output = h5py.File(filename, mode="a")
319        for field in fields:
320            output.create_dataset(field, data=self[field])
321        output.close()
322
323    def to_fits_data(self, fields=None, other_keys=None, length_unit=None, **kwargs):
324        r"""Export the fields in this FixedResolutionBuffer instance
325        to a FITSImageData instance.
326
327        This will export a set of FITS images of either the fields specified
328        or all the fields already in the object.
329
330        Parameters
331        ----------
332        fields : list of strings
333            These fields will be pixelized and output. If "None", the keys of
334            the FRB will be used.
335        other_keys : dictionary, optional
336            A set of header keys and values to write into the FITS header.
337        length_unit : string, optional
338            the length units that the coordinates are written in. The default
339            is to use the default length unit of the dataset.
340        """
341        from yt.visualization.fits_image import FITSImageData
342
343        if length_unit is None:
344            length_unit = self.ds.length_unit
345
346        if fields is None:
347            fields = list(self.data.keys())
348        else:
349            fields = list(iter_fields(fields))
350
351        if len(fields) == 0:
352            raise RuntimeError(
353                "No fields to export. Either pass a field or list of fields to "
354                "to_fits_data or access a field from the FixedResolutionBuffer "
355                "object."
356            )
357
358        fid = FITSImageData(self, fields=fields, length_unit=length_unit)
359        if other_keys is not None:
360            for k, v in other_keys.items():
361                fid.update_all_headers(k, v)
362        return fid
363
364    def export_dataset(self, fields=None, nprocs=1):
365        r"""Export a set of pixelized fields to an in-memory dataset that can be
366        analyzed as any other in yt. Unit information and other parameters (e.g.,
367        geometry, current_time, etc.) will be taken from the parent dataset.
368
369        Parameters
370        ----------
371        fields : list of strings, optional
372            These fields will be pixelized and output. If "None", the keys of the
373            FRB will be used.
374        nprocs: integer, optional
375            If greater than 1, will create this number of subarrays out of data
376
377        Examples
378        --------
379        >>> import yt
380        >>> ds = yt.load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150")
381        >>> slc = ds.slice(2, 0.0)
382        >>> frb = slc.to_frb((500.0, "kpc"), 500)
383        >>> ds2 = frb.export_dataset(
384        ...     fields=[("gas", "density"), ("gas", "temperature")], nprocs=32
385        ... )
386        """
387        nx, ny = self.buff_size
388        data = {}
389        if fields is None:
390            fields = list(self.keys())
391        for field in fields:
392            arr = self[field]
393            data[field] = (arr.d.T.reshape(nx, ny, 1), str(arr.units))
394        bounds = [b.in_units("code_length").v for b in self.bounds]
395        bbox = np.array([[bounds[0], bounds[1]], [bounds[2], bounds[3]], [0.0, 1.0]])
396        return load_uniform_grid(
397            data,
398            [nx, ny, 1],
399            length_unit=self.ds.length_unit,
400            bbox=bbox,
401            sim_time=self.ds.current_time.in_units("s").v,
402            mass_unit=self.ds.mass_unit,
403            time_unit=self.ds.time_unit,
404            velocity_unit=self.ds.velocity_unit,
405            magnetic_unit=self.ds.magnetic_unit,
406            periodicity=(False, False, False),
407            geometry=self.ds.geometry,
408            nprocs=nprocs,
409        )
410
411    def save_as_dataset(self, filename=None, fields=None):
412        r"""Export a fixed resolution buffer to a reloadable yt dataset.
413
414        This function will take a fixed resolution buffer and output a
415        dataset containing either the fields presently existing or fields
416        given in the ``fields`` list.  The resulting dataset can be
417        reloaded as a yt dataset.
418
419        Parameters
420        ----------
421        filename : str, optional
422            The name of the file to be written.  If None, the name
423            will be a combination of the original dataset and the type
424            of data container.
425        fields : list of strings or tuples, optional
426            If this is supplied, it is the list of fields to be saved to
427            disk.  If not supplied, all the fields that have been queried
428            will be saved.
429
430        Returns
431        -------
432        filename : str
433            The name of the file that has been created.
434
435        Examples
436        --------
437
438        >>> import yt
439        >>> ds = yt.load("enzo_tiny_cosmology/DD0046/DD0046")
440        >>> proj = ds.proj(("gas", "density"), "x", weight_field=("gas", "density"))
441        >>> frb = proj.to_frb(1.0, (800, 800))
442        >>> fn = frb.save_as_dataset(fields=[("gas", "density")])
443        >>> ds2 = yt.load(fn)
444        >>> print(ds2.data[("gas", "density")])
445        [[  1.25025353e-30   1.25025353e-30   1.25025353e-30 ...,   7.90820691e-31
446            7.90820691e-31   7.90820691e-31]
447         [  1.25025353e-30   1.25025353e-30   1.25025353e-30 ...,   7.90820691e-31
448            7.90820691e-31   7.90820691e-31]
449         [  1.25025353e-30   1.25025353e-30   1.25025353e-30 ...,   7.90820691e-31
450            7.90820691e-31   7.90820691e-31]
451         ...,
452         [  1.55834239e-30   1.55834239e-30   1.55834239e-30 ...,   8.51353199e-31
453            8.51353199e-31   8.51353199e-31]
454         [  1.55834239e-30   1.55834239e-30   1.55834239e-30 ...,   8.51353199e-31
455            8.51353199e-31   8.51353199e-31]
456         [  1.55834239e-30   1.55834239e-30   1.55834239e-30 ...,   8.51353199e-31
457            8.51353199e-31   8.51353199e-31]] g/cm**3
458
459        """
460
461        keyword = f"{str(self.ds)}_{self.data_source._type_name}_frb"
462        filename = get_output_filename(filename, keyword, ".h5")
463
464        data = {}
465        if fields is not None:
466            for f in self.data_source._determine_fields(fields):
467                data[f] = self[f]
468        else:
469            data.update(self.data)
470
471        ftypes = {field: "grid" for field in data}
472        extra_attrs = {
473            arg: getattr(self.data_source, arg, None)
474            for arg in self.data_source._con_args + self.data_source._tds_attrs
475        }
476        extra_attrs["con_args"] = self.data_source._con_args
477        extra_attrs["left_edge"] = self.ds.arr([self.bounds[0], self.bounds[2]])
478        extra_attrs["right_edge"] = self.ds.arr([self.bounds[1], self.bounds[3]])
479        # The data dimensions are [NY, NX] but buff_size is [NX, NY].
480        extra_attrs["ActiveDimensions"] = self.buff_size[::-1]
481        extra_attrs["level"] = 0
482        extra_attrs["data_type"] = "yt_frb"
483        extra_attrs["container_type"] = self.data_source._type_name
484        extra_attrs["dimensionality"] = self.data_source._dimensionality
485        save_as_dataset(
486            self.ds, filename, data, field_types=ftypes, extra_attrs=extra_attrs
487        )
488
489        return filename
490
491    @property
492    def limits(self):
493        rv = dict(x=None, y=None, z=None)
494        xax = self.ds.coordinates.x_axis[self.axis]
495        yax = self.ds.coordinates.y_axis[self.axis]
496        xn = self.ds.coordinates.axis_name[xax]
497        yn = self.ds.coordinates.axis_name[yax]
498        rv[xn] = (self.bounds[0], self.bounds[1])
499        rv[yn] = (self.bounds[2], self.bounds[3])
500        return rv
501
502    def setup_filters(self):
503        ignored = ["FixedResolutionBufferFilter"]
504        for key in filter_registry:
505            if key in ignored:
506                continue
507            filtername = filter_registry[key]._filter_name
508
509            # We need to wrap to create a closure so that
510            # FilterMaker is bound to the wrapped method.
511            def closure():
512                FilterMaker = filter_registry[key]
513
514                @wraps(FilterMaker)
515                def method(*args, **kwargs):
516                    # We need to also do it here as "invalidate_plot"
517                    # and "apply_callback" require the functions'
518                    # __name__ in order to work properly
519                    @wraps(FilterMaker)
520                    def cb(self, *a, **kwa):
521                        # We construct the callback method
522                        # skipping self
523                        return FilterMaker(*a, **kwa)
524
525                    # Create callback
526                    cb = apply_filter(cb)
527
528                    return cb(self, *args, **kwargs)
529
530                return method
531
532            self.__dict__["apply_" + filtername] = closure()
533
534
535class CylindricalFixedResolutionBuffer(FixedResolutionBuffer):
536    """
537    This object is a subclass of
538    :class:`yt.visualization.fixed_resolution.FixedResolutionBuffer`
539    that supports non-aligned input data objects, primarily cutting planes.
540    """
541
542    def __init__(self, data_source, radius, buff_size, antialias=True):
543
544        self.data_source = data_source
545        self.ds = data_source.ds
546        self.radius = radius
547        self.buff_size = buff_size
548        self.antialias = antialias
549        self.data = {}
550
551        ds = getattr(data_source, "ds", None)
552        if ds is not None:
553            ds.plots.append(weakref.proxy(self))
554
555    def __getitem__(self, item):
556        if item in self.data:
557            return self.data[item]
558        buff = np.zeros(self.buff_size, dtype="f8")
559        pixelize_cylinder(
560            buff,
561            self.data_source["r"],
562            self.data_source["dr"],
563            self.data_source["theta"],
564            self.data_source["dtheta"],
565            self.data_source[item].astype("float64"),
566            self.radius,
567        )
568        self[item] = buff
569        return buff
570
571
572class OffAxisProjectionFixedResolutionBuffer(FixedResolutionBuffer):
573    """
574    This object is a subclass of
575    :class:`yt.visualization.fixed_resolution.FixedResolutionBuffer`
576    that supports off axis projections.  This calls the volume renderer.
577    """
578
579    def __init__(self, data_source, bounds, buff_size, antialias=True, periodic=False):
580        self.data = {}
581        FixedResolutionBuffer.__init__(
582            self, data_source, bounds, buff_size, antialias, periodic
583        )
584
585    def __getitem__(self, item):
586        if item in self.data:
587            return self.data[item]
588        mylog.info(
589            "Making a fixed resolution buffer of (%s) %d by %d",
590            item,
591            self.buff_size[0],
592            self.buff_size[1],
593        )
594        dd = self.data_source
595        width = self.ds.arr(
596            (
597                self.bounds[1] - self.bounds[0],
598                self.bounds[3] - self.bounds[2],
599                self.bounds[5] - self.bounds[4],
600            )
601        )
602        buff = off_axis_projection(
603            dd.dd,
604            dd.center,
605            dd.normal_vector,
606            width,
607            self.buff_size,
608            item,
609            weight=dd.weight_field,
610            volume=dd.volume,
611            no_ghost=dd.no_ghost,
612            interpolated=dd.interpolated,
613            north_vector=dd.north_vector,
614            method=dd.method,
615        )
616        ia = ImageArray(buff.swapaxes(0, 1), info=self._get_info(item))
617        self[item] = ia
618        return ia
619
620
621class ParticleImageBuffer(FixedResolutionBuffer):
622    """
623
624    This object is a subclass of
625    :class:`yt.visualization.fixed_resolution.FixedResolutionBuffer`
626    that supports particle plots. It splats points onto an image
627    buffer.
628
629    """
630
631    def __init__(self, data_source, bounds, buff_size, antialias=True, periodic=False):
632        self.data = {}
633        FixedResolutionBuffer.__init__(
634            self, data_source, bounds, buff_size, antialias, periodic
635        )
636
637        # set up the axis field names
638        axis = self.axis
639        self.xax = self.ds.coordinates.x_axis[axis]
640        self.yax = self.ds.coordinates.y_axis[axis]
641        ax_field_template = "particle_position_%s"
642        self.x_field = ax_field_template % self.ds.coordinates.axis_name[self.xax]
643        self.y_field = ax_field_template % self.ds.coordinates.axis_name[self.yax]
644
645    def __getitem__(self, item):
646        if item in self.data:
647            return self.data[item]
648
649        deposition = self.data_source.deposition
650        density = self.data_source.density
651
652        mylog.info(
653            "Splatting (%s) onto a %d by %d mesh using method '%s'",
654            item,
655            self.buff_size[0],
656            self.buff_size[1],
657            deposition,
658        )
659
660        bounds = []
661        for b in self.bounds:
662            if hasattr(b, "in_units"):
663                b = float(b.in_units("code_length"))
664            bounds.append(b)
665
666        ftype = item[0]
667        x_data = self.data_source.dd[ftype, self.x_field]
668        y_data = self.data_source.dd[ftype, self.y_field]
669        data = self.data_source.dd[item]
670
671        # handle periodicity
672        dx = x_data.in_units("code_length").d - bounds[0]
673        dy = y_data.in_units("code_length").d - bounds[2]
674        if self.periodic:
675            dx %= float(self._period[0].in_units("code_length"))
676            dy %= float(self._period[1].in_units("code_length"))
677
678        # convert to pixels
679        px = dx / (bounds[1] - bounds[0])
680        py = dy / (bounds[3] - bounds[2])
681
682        # select only the particles that will actually show up in the image
683        mask = np.logical_and(
684            np.logical_and(px >= 0.0, px <= 1.0), np.logical_and(py >= 0.0, py <= 1.0)
685        )
686
687        weight_field = self.data_source.weight_field
688        if weight_field is None:
689            weight_data = np.ones_like(data.v)
690        else:
691            weight_data = self.data_source.dd[weight_field]
692        splat_vals = weight_data[mask] * data[mask]
693
694        x_bin_edges = np.linspace(0.0, 1.0, self.buff_size[0] + 1)
695        y_bin_edges = np.linspace(0.0, 1.0, self.buff_size[1] + 1)
696
697        # splat particles
698        buff = np.zeros(self.buff_size)
699        buff_mask = np.zeros_like(buff, dtype="uint8")
700        if deposition == "ngp":
701            add_points_to_greyscale_image(
702                buff, buff_mask, px[mask], py[mask], splat_vals
703            )
704        elif deposition == "cic":
705            CICDeposit_2(
706                py[mask],
707                px[mask],
708                splat_vals,
709                mask.sum(),
710                buff,
711                buff_mask,
712                x_bin_edges,
713                y_bin_edges,
714            )
715        else:
716            raise ValueError(f"Received unknown deposition method '{deposition}'")
717
718        # remove values in no-particle region
719        buff[buff_mask == 0] = np.nan
720
721        # Normalize by the surface area of the pixel or volume of pencil if
722        # requested
723        info = self._get_info(item)
724        if density:
725            width = self.data_source.width
726            norm = width[self.xax] * width[self.yax] / np.prod(self.buff_size)
727            norm = norm.in_base()
728            buff /= norm.v
729            units = data.units / norm.units
730            info["label"] = "%s $\\rm{Density}$" % info["label"]
731        else:
732            units = data.units
733
734        ia = ImageArray(buff, units=units, info=info)
735
736        # divide by the weight_field, if needed
737        if weight_field is not None:
738            weight_buff = np.zeros(self.buff_size)
739            weight_buff_mask = np.zeros(self.buff_size, dtype="uint8")
740            if deposition == "ngp":
741                add_points_to_greyscale_image(
742                    weight_buff, weight_buff_mask, px[mask], py[mask], weight_data[mask]
743                )
744            elif deposition == "cic":
745                CICDeposit_2(
746                    py[mask],
747                    px[mask],
748                    weight_data[mask],
749                    mask.sum(),
750                    weight_buff,
751                    weight_buff_mask,
752                    y_bin_edges,
753                    x_bin_edges,
754                )
755            weight_array = ImageArray(
756                weight_buff, units=weight_data.units, info=self._get_info(item)
757            )
758            # remove values in no-particle region
759            weight_buff[weight_buff_mask == 0] = np.nan
760            locs = np.where(weight_array > 0)
761            ia[locs] /= weight_array[locs]
762
763        self.data[item] = ia
764        return self.data[item]
765
766    # over-ride the base class version, since we don't want to exclude
767    # particle fields
768    def _get_data_source_fields(self):
769        exclude = self.data_source._key_fields + list(self._exclude_fields)
770        fields = getattr(self.data_source, "fields", [])
771        fields += getattr(self.data_source, "field_data", {}).keys()
772        for f in fields:
773            if f not in exclude:
774                self[f]
775