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