diff --git a/doc/source/visualizing/FITSImageData.ipynb b/doc/source/visualizing/FITSImageData.ipynb index a96677b3c05..75526640267 100644 --- a/doc/source/visualizing/FITSImageData.ipynb +++ b/doc/source/visualizing/FITSImageData.ipynb @@ -133,10 +133,7 @@ { "cell_type": "markdown", "metadata": { - "collapsed": false, - "pycharm": { - "name": "#%% md\n" - } + "collapsed": false }, "source": [ "## Making FITS images from Particle Projections" @@ -149,10 +146,7 @@ "`FITSParticleProjection`:" ], "metadata": { - "collapsed": false, - "pycharm": { - "name": "#%% md\n" - } + "collapsed": false } }, { @@ -167,10 +161,7 @@ "prjp_fits.writeto(\"prjp.fits\", overwrite=True)" ], "metadata": { - "collapsed": false, - "pycharm": { - "name": "#%%\n" - } + "collapsed": false } }, { @@ -183,10 +174,7 @@ "for example), supply the ``density`` keyword argument:" ], "metadata": { - "collapsed": false, - "pycharm": { - "name": "#%% md\n" - } + "collapsed": false } }, { @@ -200,10 +188,32 @@ "prjpd_fits.writeto(\"prjpd.fits\", overwrite=True)" ], "metadata": { - "collapsed": false, - "pycharm": { - "name": "#%%\n" - } + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "`FITSParticleOffAxisProjection` can be used to make a projection along any arbitrary sight line:" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "L = [1, -1, 1] # normal or \"line of sight\" vector\n", + "N = [0, 0, 1] # north or \"up\" vector\n", + "poff_fits = yt.FITSParticleOffAxisProjection(\n", + " dsp, L, (\"PartType1\", \"particle_mass\"), deposition=\"cic\", north_vector=N\n", + ")\n", + "poff_fits.writeto(\"poff.fits\", overwrite=True)" + ], + "metadata": { + "collapsed": false } }, { @@ -212,10 +222,7 @@ "## Using `HDUList` Methods" ], "metadata": { - "collapsed": false, - "pycharm": { - "name": "#%% md\n" - } + "collapsed": false } }, { diff --git a/doc/source/visualizing/plots.rst b/doc/source/visualizing/plots.rst index fe4c2e49d86..80978129bc8 100644 --- a/doc/source/visualizing/plots.rst +++ b/doc/source/visualizing/plots.rst @@ -2027,7 +2027,25 @@ domain: p.set_unit(("all", "particle_mass"), "Msun") p.save() -and here is an example of using the ``data_source`` argument to :class:`~yt.visualization.particle_plots.ParticlePhasePlot` +Using :class:`~yt.visualization.particle_plots.ParticleProjectionPlot`, you can also plot particles +along an off-axis direction: + +.. python-script:: + + import yt + + ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030") + + L = [1, 1, 1] # normal or "line of sight" vector + N = [0, 1, 0] # north or "up" vector + + p = yt.ParticleProjectionPlot( + ds, L, [("all", "particle_mass")], width=(0.05, 0.05), depth=0.3, north_vector=N + ) + p.set_unit(("all", "particle_mass"), "Msun") + p.save() + +Here is an example of using the ``data_source`` argument to :class:`~yt.visualization.particle_plots.ParticlePhasePlot` to only consider the particles that lie within a 50 kpc sphere around the domain center: .. python-script:: diff --git a/tests/tests.yaml b/tests/tests.yaml index 007aba334f3..bbc93148102 100644 --- a/tests/tests.yaml +++ b/tests/tests.yaml @@ -102,9 +102,10 @@ answer_tests: - yt/frontends/owls/tests/test_outputs.py:test_snapshot_033 - yt/frontends/owls/tests/test_outputs.py:test_OWLS_particlefilter - local_pw_048: # PR 4417 + local_pw_049: # PR 4440 - yt/visualization/tests/test_plotwindow.py:test_attributes - yt/visualization/tests/test_particle_plot.py:test_particle_projection_answers + - yt/visualization/tests/test_particle_plot.py:test_particle_offaxis_projection_answers - yt/visualization/tests/test_particle_plot.py:test_particle_projection_filter - yt/visualization/tests/test_particle_plot.py:test_particle_phase_answers - yt/visualization/tests/test_callbacks.py:test_axis_manipulations diff --git a/yt/__init__.py b/yt/__init__.py index d54931f7da8..53094f87b84 100644 --- a/yt/__init__.py +++ b/yt/__init__.py @@ -118,6 +118,7 @@ def run_nose(*args, **kwargs): FITSImageData, FITSOffAxisProjection, FITSOffAxisSlice, + FITSParticleOffAxisProjection, FITSParticleProjection, FITSProjection, FITSSlice, diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 4689bff6cd9..2511d9bd9be 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1781,24 +1781,13 @@ def pixelize_element_mesh_line(np.ndarray[np.float64_t, ndim=2] coords, @cython.boundscheck(False) @cython.wraparound(False) -def off_axis_projection_SPH(np.float64_t[:] px, - np.float64_t[:] py, - np.float64_t[:] pz, - np.float64_t[:] particle_masses, - np.float64_t[:] particle_densities, - np.float64_t[:] smoothing_lengths, - bounds, - center, - width, - np.float64_t[:] quantity_to_smooth, - np.float64_t[:, :] projection_array, - normal_vector, - north_vector, - weight_field=None): - # Do nothing in event of a 0 normal vector - if np.allclose(normal_vector, np.array([0., 0., 0.]), rtol=1e-09): - return - +def rotate_particle_coord(np.float64_t[:] px, + np.float64_t[:] py, + np.float64_t[:] pz, + center, + width, + normal_vector, + north_vector): # We want to do two rotations, one to first rotate our coordinates to have # the normal vector be the z-axis (i.e., the viewer's perspective), and then # another rotation to make the north-vector be the y-axis (i.e., north). @@ -1839,6 +1828,34 @@ def off_axis_projection_SPH(np.float64_t[:] px, px_rotated[i] = rotated_coordinates[0] py_rotated[i] = rotated_coordinates[1] + return px_rotated, py_rotated, rot_bounds_x0, rot_bounds_x1, rot_bounds_y0, rot_bounds_y1 + + +@cython.boundscheck(False) +@cython.wraparound(False) +def off_axis_projection_SPH(np.float64_t[:] px, + np.float64_t[:] py, + np.float64_t[:] pz, + np.float64_t[:] particle_masses, + np.float64_t[:] particle_densities, + np.float64_t[:] smoothing_lengths, + bounds, + center, + width, + np.float64_t[:] quantity_to_smooth, + np.float64_t[:, :] projection_array, + normal_vector, + north_vector, + weight_field=None): + # Do nothing in event of a 0 normal vector + if np.allclose(normal_vector, 0.): + return + + px_rotated, py_rotated, \ + rot_bounds_x0, rot_bounds_x1, \ + rot_bounds_y0, rot_bounds_y1 = rotate_particle_coord(px, py, pz, + center, width, normal_vector, north_vector) + pixelize_sph_kernel_projection(projection_array, px_rotated, py_rotated, diff --git a/yt/visualization/api.py b/yt/visualization/api.py index 12bc0e6019b..aaaa347b2f4 100644 --- a/yt/visualization/api.py +++ b/yt/visualization/api.py @@ -4,6 +4,7 @@ FITSImageData, FITSOffAxisProjection, FITSOffAxisSlice, + FITSParticleOffAxisProjection, FITSParticleProjection, FITSProjection, FITSSlice, diff --git a/yt/visualization/fits_image.py b/yt/visualization/fits_image.py index 77e35c1309a..85a098fa30e 100644 --- a/yt/visualization/fits_image.py +++ b/yt/visualization/fits_image.py @@ -18,7 +18,11 @@ from yt.utilities.on_demand_imports import _astropy from yt.utilities.parallel_tools.parallel_analysis_interface import parallel_root_only from yt.visualization.fixed_resolution import FixedResolutionBuffer, ParticleImageBuffer -from yt.visualization.particle_plots import ParticleAxisAlignedDummyDataSource +from yt.visualization.particle_plots import ( + ParticleAxisAlignedDummyDataSource, + ParticleDummyDataSource, + ParticleOffAxisDummyDataSource, +) from yt.visualization.volume_rendering.off_axis_projection import off_axis_projection @@ -906,17 +910,23 @@ def construct_image( frb = data_source.to_frb(width[0], (nx, ny), height=width[1]) else: frb = data_source.to_frb(width[0], (nx, ny), center=center, height=width[1]) - elif isinstance(data_source, ParticleAxisAlignedDummyDataSource): - axes = axis_wcs[axis] - bounds = ( - center[axes[0]] - width[0] / 2, - center[axes[0]] + width[0] / 2, - center[axes[1]] - width[1] / 2, - center[axes[1]] + width[1] / 2, - ) - frb = ParticleImageBuffer( - data_source, bounds, (nx, ny), periodic=all(ds.periodicity) - ) + elif isinstance(data_source, ParticleDummyDataSource): + if hasattr(data_source, "normal_vector"): + # If we have a normal vector, this means + # that the data source is off-axis + bounds = (-width[0] / 2, width[0] / 2, -width[1] / 2, width[1] / 2) + periodic = False + else: + # Otherwise, this is an on-axis data source + axes = axis_wcs[axis] + bounds = ( + center[axes[0]] - width[0] / 2, + center[axes[0]] + width[0] / 2, + center[axes[1]] - width[1] / 2, + center[axes[1]] + width[1] / 2, + ) + periodic = all(ds.periodicity) + frb = ParticleImageBuffer(data_source, bounds, (nx, ny), periodic=periodic) else: frb = None w = _astropy.pywcs.WCS(naxis=2) @@ -1107,8 +1117,8 @@ class FITSProjection(FITSImageData): For example, (10, 'kpc') specifies a width that is 10 kiloparsecs wide in the x and y directions, ((10,'kpc'),(15,'kpc')) specifies a - width that is 10 kiloparsecs wide along the x axis and 15 - kiloparsecs wide along the y axis. In the other two examples, code + width that is 10 kiloparsecs wide along the x-axis and 15 + kiloparsecs wide along the y-axis. In the other two examples, code units are assumed, for example (0.2, 0.3) specifies a width that has an x width of 0.2 and a y width of 0.3 in code units. weight_field : string @@ -1214,8 +1224,8 @@ class FITSParticleProjection(FITSImageData): For example, (10, 'kpc') specifies a width that is 10 kiloparsecs wide in the x and y directions, ((10,'kpc'),(15,'kpc')) specifies a - width that is 10 kiloparsecs wide along the x axis and 15 - kiloparsecs wide along the y axis. In the other two examples, code + width that is 10 kiloparsecs wide along the x-axis and 15 + kiloparsecs wide along the y-axis. In the other two examples, code units are assumed, for example (0.2, 0.3) specifies a width that has an x width of 0.2 and a y width of 0.3 in code units. depth : A tuple or a float @@ -1280,7 +1290,7 @@ def __init__( axis, width, fields, - weight_field, + weight_field=weight_field, field_parameters=field_parameters, data_source=data_source, deposition=deposition, @@ -1292,6 +1302,132 @@ def __init__( super().__init__(frb, fields=fields, length_unit=lunit, wcs=w) +class FITSParticleOffAxisProjection(FITSImageData): + r""" + Generate a FITSImageData of an off-axis projection of a + particle field. + + Parameters + ---------- + ds : :class:`~yt.data_objects.static_output.Dataset` + The dataset object. + axis : character or integer + The axis along which to project. One of "x","y","z", or 0,1,2. + fields : string or list of strings + The fields to project + image_res : an int or 2-tuple of ints + Specify the resolution of the resulting image. A single value will be + used for both axes, whereas a tuple of values will be used for the + individual axes. Default: 512 + center : A sequence of floats, a string, or a tuple. + The coordinate of the center of the image. If set to 'c', 'center' or + left blank, the plot is centered on the middle of the domain. If set + to 'max' or 'm', the center will be located at the maximum of the + ('gas', 'density') field. Centering on the max or min of a specific + field is supported by providing a tuple such as ("min","temperature") + or ("max","dark_matter_density"). Units can be specified by passing in + *center* as a tuple containing a coordinate and string unit name or by + passing in a YTArray. If a list or unitless array is supplied, code + units are assumed. + width : tuple or a float. + Width can have four different formats to support variable + x and y widths. They are: + + ================================== ======================= + format example + ================================== ======================= + (float, string) (10,'kpc') + ((float, string), (float, string)) ((10,'kpc'),(15,'kpc')) + float 0.2 + (float, float) (0.2, 0.3) + ================================== ======================= + + For example, (10, 'kpc') specifies a width that is 10 kiloparsecs + wide in the x and y directions, ((10,'kpc'),(15,'kpc')) specifies a + width that is 10 kiloparsecs wide along the x-axis and 15 + kiloparsecs wide along the y-axis. In the other two examples, code + units are assumed, for example (0.2, 0.3) specifies a width that has an + x width of 0.2 and a y width of 0.3 in code units. + depth : A tuple or a float + A tuple containing the depth to project through and the string + key of the unit: (width, 'unit'). If set to a float, code units + are assumed. Defaults to the entire domain. + weight_field : string + The field used to weight the projection. + length_unit : string, optional + the length units that the coordinates are written in. The default + is to use the default length unit of the dataset. + deposition : string, optional + Controls the order of the interpolation of the particles onto the + mesh. "ngp" is 0th-order "nearest-grid-point" method (the default), + "cic" is 1st-order "cloud-in-cell". + density : boolean, optional + If True, the quantity to be projected will be divided by the area of + the cells, to make a projected density of the quantity. Default: False + field_parameters : dictionary + A dictionary of field parameters than can be accessed by derived + fields. + data_source : yt.data_objects.data_containers.YTSelectionContainer, optional + If specified, this will be the data source used for selecting regions + to project. + origin : string + The origin of the coordinate system in the file. If "domain", then the + center coordinates will be the same as the center of the image as + defined by the *center* keyword argument. If "image", then the center + coordinates will be set to (0,0). Default: "domain" + """ + + def __init__( + self, + ds, + normal, + fields, + image_res=512, + center="c", + width=None, + depth=(1, "1"), + weight_field=None, + length_unit=None, + deposition="ngp", + density=False, + field_parameters=None, + data_source=None, + north_vector=None, + ): + fields = list(iter_fields(fields)) + center, dcenter = ds.coordinates.sanitize_center(center, None) + width = ds.coordinates.sanitize_width(normal, width, depth) + wd = tuple(el.in_units("code_length").v for el in width) + + if field_parameters is None: + field_parameters = {} + + ps = ParticleOffAxisDummyDataSource( + center, + ds, + normal, + wd, + fields, + weight_field=weight_field, + field_parameters=field_parameters, + data_source=data_source, + deposition=deposition, + density=density, + north_vector=north_vector, + ) + w, frb, lunit = construct_image( + ds, + normal, + ps, + dcenter, + image_res, + width, + length_unit, + origin="image", + ) + super().__init__(frb, fields=fields, length_unit=lunit, wcs=w) + + class FITSOffAxisSlice(FITSImageData): r""" Generate a FITSImageData of an off-axis slice. @@ -1426,8 +1562,8 @@ class FITSOffAxisProjection(FITSImageData): For example, (10, 'kpc') specifies a width that is 10 kiloparsecs wide in the x and y directions, ((10,'kpc'),(15,'kpc')) specifies a - width that is 10 kiloparsecs wide along the x axis and 15 - kiloparsecs wide along the y axis. In the other two examples, code + width that is 10 kiloparsecs wide along the x-axis and 15 + kiloparsecs wide along the y-axis. In the other two examples, code units are assumed, for example (0.2, 0.3) specifies a width that has an x width of 0.2 and a y width of 0.3 in code units. depth : A tuple or a float diff --git a/yt/visualization/fixed_resolution.py b/yt/visualization/fixed_resolution.py index 83326e3daf0..8650ccec0bc 100644 --- a/yt/visualization/fixed_resolution.py +++ b/yt/visualization/fixed_resolution.py @@ -14,7 +14,10 @@ CICDeposit_2, add_points_to_greyscale_image, ) -from yt.utilities.lib.pixelization_routines import pixelize_cylinder +from yt.utilities.lib.pixelization_routines import ( + pixelize_cylinder, + rotate_particle_coord, +) from yt.utilities.math_utils import compute_stddev_image from yt.utilities.on_demand_imports import _h5py as h5py @@ -687,11 +690,12 @@ def __init__( # set up the axis field names axis = self.axis - self.xax = self.ds.coordinates.x_axis[axis] - self.yax = self.ds.coordinates.y_axis[axis] - ax_field_template = "particle_position_%s" - self.x_field = ax_field_template % self.ds.coordinates.axis_name[self.xax] - self.y_field = ax_field_template % self.ds.coordinates.axis_name[self.yax] + if axis is not None: + self.xax = self.ds.coordinates.x_axis[axis] + self.yax = self.ds.coordinates.y_axis[axis] + axis_name = self.ds.coordinates.axis_name + self.x_field = f"particle_position_{axis_name[self.xax]}" + self.y_field = f"particle_position_{axis_name[self.yax]}" def __getitem__(self, item): if item in self.data: @@ -708,20 +712,39 @@ def __getitem__(self, item): deposition, ) - bounds = [] - for b in self.bounds: - if hasattr(b, "in_units"): - b = float(b.in_units("code_length")) - bounds.append(b) + dd = self.data_source.dd ftype = item[0] - x_data = self.data_source.dd[ftype, self.x_field] - y_data = self.data_source.dd[ftype, self.y_field] - data = self.data_source.dd[item] + if self.axis is None: + wd = [] + for w in self.data_source.width: + if hasattr(w, "to_value"): + w = w.to_value("code_length") + wd.append(w) + x_data, y_data, *bounds = rotate_particle_coord( + dd[ftype, "particle_position_x"].to_value("code_length"), + dd[ftype, "particle_position_y"].to_value("code_length"), + dd[ftype, "particle_position_z"].to_value("code_length"), + self.data_source.center.to_value("code_length"), + wd, + self.data_source.normal_vector, + self.data_source.north_vector, + ) + x_data = np.array(x_data) + y_data = np.array(y_data) + else: + bounds = [] + for b in self.bounds: + if hasattr(b, "to_value"): + b = b.to_value("code_length") + bounds.append(b) + x_data = dd[ftype, self.x_field].to_value("code_length") + y_data = dd[ftype, self.y_field].to_value("code_length") + data = dd[item] # handle periodicity - dx = x_data.in_units("code_length").d - bounds[0] - dy = y_data.in_units("code_length").d - bounds[2] + dx = x_data - bounds[0] + dy = y_data - bounds[2] if self.periodic: dx %= float(self._period[0].in_units("code_length")) dy %= float(self._period[1].in_units("code_length")) @@ -739,7 +762,7 @@ def __getitem__(self, item): if weight_field is None: weight_data = np.ones_like(data.v) else: - weight_data = self.data_source.dd[weight_field] + weight_data = dd[weight_field] splat_vals = weight_data[mask] * data[mask] x_bin_edges = np.linspace(0.0, 1.0, self.buff_size[0] + 1) diff --git a/yt/visualization/particle_plots.py b/yt/visualization/particle_plots.py index dc2fee71c49..35bdf3bcfad 100644 --- a/yt/visualization/particle_plots.py +++ b/yt/visualization/particle_plots.py @@ -2,17 +2,25 @@ import numpy as np +from yt._maintenance.deprecation import issue_deprecation_warning from yt.data_objects.profiles import create_profile from yt.data_objects.static_output import Dataset -from yt.funcs import fix_axis, iter_fields +from yt.funcs import fix_axis, iter_fields, mylog from yt.units.yt_array import YTArray +from yt.utilities.orientation import Orientation from yt.visualization.fixed_resolution import ParticleImageBuffer from yt.visualization.profile_plotter import PhasePlot -from .plot_window import PWViewerMPL, get_axes_unit, get_window_parameters +from .plot_window import ( + NormalPlot, + PWViewerMPL, + get_axes_unit, + get_oblique_window_parameters, + get_window_parameters, +) -class ParticleAxisAlignedDummyDataSource: +class ParticleDummyDataSource: _type_name = "Particle" _dimensionality = 2 _con_args = ("center", "axis", "width", "fields", "weight_field") @@ -23,20 +31,24 @@ def __init__( self, center, ds, - axis, width, fields, + dd, + *, weight_field=None, field_parameters=None, - data_source=None, deposition="ngp", density=False, ): self.center = center self.ds = ds - self.axis = axis self.width = width + self.dd = dd + + if weight_field is not None: + weight_field = self._determine_fields(weight_field)[0] self.weight_field = weight_field + self.deposition = deposition self.density = density @@ -45,6 +57,40 @@ def __init__( else: self.field_parameters = field_parameters + fields = self._determine_fields(fields) + self.fields = fields + + def _determine_fields(self, *args): + return self.dd._determine_fields(*args) + + def get_field_parameter(self, name, default=None): + """ + This is typically only used by derived field functions, but + it returns parameters used to generate fields. + """ + if name in self.field_parameters: + return self.field_parameters[name] + else: + return default + + +class ParticleAxisAlignedDummyDataSource(ParticleDummyDataSource): + def __init__( + self, + center, + ds, + axis, + width, + fields, + *, + weight_field=None, + field_parameters=None, + data_source=None, + deposition="ngp", + density=False, + ): + self.axis = axis + LE = center - 0.5 * YTArray(width) RE = center + 0.5 * YTArray(width) for ax in range(3): @@ -52,7 +98,7 @@ def __init__( LE[ax] = max(LE[ax], ds.domain_left_edge[ax]) RE[ax] = min(RE[ax], ds.domain_right_edge[ax]) - self.dd = ds.region( + dd = ds.region( center, LE, RE, @@ -61,27 +107,79 @@ def __init__( data_source=data_source, ) - fields = self.dd._determine_fields(fields) - self.fields = fields + super().__init__( + center, + ds, + width, + fields, + dd, + weight_field=weight_field, + field_parameters=field_parameters, + deposition=deposition, + density=density, + ) - def _determine_fields(self, *args): - return self.dd._determine_fields(*args) - def get_field_parameter(self, name, default=None): - """ - This is typically only used by derived field functions, but - it returns parameters used to generate fields. - """ - if name in self.field_parameters: - return self.field_parameters[name] +class ParticleOffAxisDummyDataSource(ParticleDummyDataSource): + def __init__( + self, + center, + ds, + normal_vector, + width, + fields, + *, + weight_field=None, + field_parameters=None, + data_source=None, + deposition="ngp", + density=False, + north_vector=None, + ): + self.axis = None # always true for oblique data objects + normal = np.array(normal_vector) + normal = normal / np.linalg.norm(normal) + + # If north_vector is None, we set the default here. + # This is chosen so that if normal_vector is one of the + # cartesian coordinate axes, the projection will match + # the corresponding on-axis projection. + if north_vector is None: + vecs = np.identity(3) + t = np.cross(vecs, normal).sum(axis=1) + ax = t.argmax() + east_vector = np.cross(vecs[ax, :], normal).ravel() + north = np.cross(normal, east_vector).ravel() else: - return default + north = np.array(north_vector) + north = north / np.linalg.norm(north) + self.normal_vector = normal + self.north_vector = north + + if data_source is None: + dd = ds.all_data() + else: + dd = data_source + + self.orienter = Orientation(normal_vector, north_vector=north_vector) + + super().__init__( + center, + ds, + width, + fields, + dd, + weight_field=weight_field, + field_parameters=field_parameters, + deposition=deposition, + density=density, + ) -class ParticleProjectionPlot(PWViewerMPL): +class ParticleProjectionPlot(PWViewerMPL, NormalPlot): r"""Creates a particle plot from a dataset - Given a ds object, an axis to slice along, and a field name + Given a ds object, a normal to project along, and a field name string, this will return a PWViewerMPL object containing the plot. @@ -93,9 +191,12 @@ class ParticleProjectionPlot(PWViewerMPL): ds : `Dataset` This is the dataset object corresponding to the simulation output to be plotted. - axis : int or one of 'x', 'y', 'z' - An int corresponding to the axis to slice along (0=x, 1=y, 2=z) - or the axis name itself + normal : int, str, or 3-element sequence of floats + This specifies the normal vector to the projection. + Valid int values are 0, 1 and 2. Corresponding str values depend on the + geometry of the dataset and are generally given by `ds.coordinates.axis_order`. + E.g. in cartesian they are 'x', 'y' and 'z'. + An arbitrary normal vector may be specified as a 3-element sequence of floats. fields : string, list or None If a string or list, the name of the particle field(s) to be used on the colorbar. The color shown will correspond to the sum of the @@ -142,8 +243,8 @@ class ParticleProjectionPlot(PWViewerMPL): For example, (10, 'kpc') requests a plot window that is 10 kiloparsecs wide in the x and y directions, ((10,'kpc'),(15,'kpc')) requests a - window that is 10 kiloparsecs wide along the x axis and 15 - kiloparsecs wide along the y axis. In the other two examples, code + window that is 10 kiloparsecs wide along the x-axis and 15 + kiloparsecs wide along the y-axis. In the other two examples, code units are assumed, for example (0.2, 0.3) requests a plot that has an x width of 0.2 and a y width of 0.3 in code units. If units are provided the resulting plot axis labels will use the supplied units. @@ -165,7 +266,7 @@ class ParticleProjectionPlot(PWViewerMPL): represented by '-' separated string or a tuple of strings. In the first index the y-location is given by 'lower', 'upper', or 'center'. The second index is the x-location, given as 'left', 'right', or - 'center'. Finally, the whether the origin is applied in 'domain' + 'center'. Finally, whether the origin is applied in 'domain' space, plot 'window' space or 'native' simulation coordinate system is given. For example, both 'upper-right-domain' and ['upper', 'right', 'domain'] both place the origin in the upper right hand @@ -197,8 +298,8 @@ class ParticleProjectionPlot(PWViewerMPL): including the margins but not the colorbar. aspect : float The aspect ratio of the plot. Set to None for 1. - data_source : YTSelectionContainer Object - Object to be used for data selection. Defaults to a region covering + data_source : YTSelectionContainer object + The object to be used for data selection. Defaults to a region covering the entire simulation. deposition : string Controls the order of the interpolation of the particles onto the @@ -208,11 +309,16 @@ class ParticleProjectionPlot(PWViewerMPL): If True, the quantity to be projected will be divided by the area of the cells, to make a projected density of the quantity. The plot name and units will also reflect this. Default: False + north_vector : a sequence of floats + A vector defining the 'up' direction in off-axis particle projection plots; + not used if the plot is on-axis. This option sets the orientation of the + projected plane. If not set, an arbitrary grid-aligned north-vector is + chosen. Examples -------- - This will save an image the the file + This will save an image to the file 'galaxy0030_Particle_z_particle_mass.png' >>> from yt import load @@ -227,7 +333,7 @@ class ParticleProjectionPlot(PWViewerMPL): def __init__( self, ds, - axis, + normal=None, fields=None, color="b", center="center", @@ -243,15 +349,25 @@ def __init__( data_source=None, deposition="ngp", density=False, + *, + north_vector=None, + axis=None, ): + if axis is not None: + issue_deprecation_warning( + "The 'axis' argument is a deprecated alias for the 'normal' argument. ", + since="4.2.0", + ) + normal = axis + if normal is None: + raise TypeError( + "ParticleProjectionPlot() missing 1 required positional argument: 'normal'" + ) # this will handle time series data and controllers ts = self._initialize_dataset(ds) self.ts = ts ds = self.ds = ts[0] - axis = fix_axis(axis, ds) - (bounds, center, display_center) = get_window_parameters( - axis, center, width, ds - ) + normal = self.sanitize_normal_vector(ds, normal) if field_parameters is None: field_parameters = {} @@ -270,28 +386,67 @@ def __init__( self._use_cbar = False splat_color = color - depth = ds.coordinates.sanitize_depth(depth) + if isinstance(normal, str): + axis = fix_axis(normal, ds) + (bounds, center, display_center) = get_window_parameters( + axis, center, width, ds + ) + x_coord = ds.coordinates.x_axis[axis] + y_coord = ds.coordinates.y_axis[axis] + + depth = ds.coordinates.sanitize_depth(depth) + + width = np.zeros_like(center) + width[x_coord] = bounds[1] - bounds[0] + width[y_coord] = bounds[3] - bounds[2] + width[axis] = depth[0].in_units(width[x_coord].units) + + ParticleSource = ParticleAxisAlignedDummyDataSource( + center, + ds, + axis, + width, + fields, + weight_field=weight_field, + field_parameters=field_parameters, + data_source=data_source, + deposition=deposition, + density=density, + ) + + oblique = False + plt_origin = origin + periodic = True - x_coord = ds.coordinates.x_axis[axis] - y_coord = ds.coordinates.y_axis[axis] - - width = np.zeros_like(center) - width[x_coord] = bounds[1] - bounds[0] - width[y_coord] = bounds[3] - bounds[2] - width[axis] = depth[0].in_units(width[x_coord].units) - - ParticleSource = ParticleAxisAlignedDummyDataSource( - center, - ds, - axis, - width, - fields, - weight_field, - field_parameters=field_parameters, - data_source=data_source, - deposition=deposition, - density=density, - ) + else: + (bounds, center_rot) = get_oblique_window_parameters( + normal, center, width, ds, depth=depth + ) + + width = ds.coordinates.sanitize_width(normal, width, depth) + + ParticleSource = ParticleOffAxisDummyDataSource( + center_rot, + ds, + normal, + width, + fields, + weight_field=weight_field, + field_parameters=field_parameters, + data_source=None, + deposition=deposition, + density=density, + north_vector=north_vector, + ) + + oblique = True + periodic = False + if origin != "center-window": + mylog.warning( + "The 'origin' keyword is ignored for off-axis " + "particle projections, it is always 'center-window'" + ) + plt_origin = "center-window" self.projected = weight_field is None @@ -299,13 +454,15 @@ def __init__( self, ParticleSource, bounds, - origin=origin, + origin=plt_origin, fontsize=fontsize, fields=fields, window_size=window_size, aspect=aspect, splat_color=splat_color, geometry=ds.geometry, + periodic=periodic, + oblique=oblique, ) self.set_axes_unit(axes_unit) diff --git a/yt/visualization/plot_window.py b/yt/visualization/plot_window.py index 5c33be9c107..9e1fee7bcd7 100644 --- a/yt/visualization/plot_window.py +++ b/yt/visualization/plot_window.py @@ -1618,7 +1618,7 @@ class ProjectionPlot(NormalPlot): This is the dataset object corresponding to the simulation output to be plotted. normal : int, str, or 3-element sequence of floats - This specifies the normal vector to the slice. + This specifies the normal vector to the projection. Valid int values are 0, 1 and 2. Corresponding str values depend on the geometry of the dataset and are generally given by `ds.coordinates.axis_order`. E.g. in cartesian they are 'x', 'y' and 'z'. diff --git a/yt/visualization/tests/test_fits_image.py b/yt/visualization/tests/test_fits_image.py index 8fe96f8d743..2ba606d0253 100644 --- a/yt/visualization/tests/test_fits_image.py +++ b/yt/visualization/tests/test_fits_image.py @@ -12,6 +12,7 @@ FITSImageData, FITSOffAxisProjection, FITSOffAxisSlice, + FITSParticleOffAxisProjection, FITSParticleProjection, FITSProjection, FITSSlice, @@ -221,6 +222,11 @@ def test_fits_image(): assert pfid["particle_mass"].header["BTYPE"] == "particle_mass" assert pfid["particle_mass"].units == "g" + pofid = FITSParticleOffAxisProjection(ds, [1, 1, 1], ("io", "particle_mass")) + assert pofid["particle_mass"].name == "particle_mass" + assert pofid["particle_mass"].header["BTYPE"] == "particle_mass" + assert pofid["particle_mass"].units == "g" + pdfid = FITSParticleProjection(ds, "x", ("io", "particle_mass"), density=True) assert pdfid["particle_mass"].name == "particle_mass" assert pdfid["particle_mass"].header["BTYPE"] == "particle_mass" @@ -282,6 +288,7 @@ def _vlsq(field, data): new_fid3, pfid, pdfid, + pofid, ): fid.close() diff --git a/yt/visualization/tests/test_particle_plot.py b/yt/visualization/tests/test_particle_plot.py index 0a627b4b0c3..0716b8bd4a6 100644 --- a/yt/visualization/tests/test_particle_plot.py +++ b/yt/visualization/tests/test_particle_plot.py @@ -119,6 +119,27 @@ def test_particle_projection_answers(): yield test +@requires_ds(g30, big_data=True) +def test_particle_offaxis_projection_answers(): + plot_field = ("all", "particle_mass") + decimals = 12 + ds = data_dir_load(g30) + attr_name = "set_cmap" + attr_args = ((("all", "particle_mass"), "RdBu"), {}) + L = [1, 1, 1] + test = PlotWindowAttributeTest( + ds, + plot_field, + L, + attr_name, + attr_args, + decimals, + "ParticleProjectionPlot", + ) + test_particle_offaxis_projection_answers.__name__ = test.description + yield test + + @requires_ds(g30, big_data=True) def test_particle_projection_filter(): """ @@ -410,6 +431,19 @@ def test_particle_plot_wf(self): ): pplot_wf.save() + def test_particle_plot_offaxis(self): + test_ds = fake_particle_ds() + Ls = [[1, 1, 1], [0, 1, -0.5]] + Ns = [None, [1, 1, 1]] + for L, N in zip(Ls, Ns): + pplot_off = ParticleProjectionPlot( + test_ds, L, ("all", "particle_mass"), north_vector=N + ) + with mock.patch( + "yt.visualization._mpl_imports.FigureCanvasAgg.print_figure" + ): + pplot_off.save() + def test_creation_with_width(self): test_ds = fake_particle_ds() for width, (xlim, ylim, pwidth, _aun) in WIDTH_SPECS.items():