[yt-svn] commit/yt: 3 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Nov 6 08:44:48 PST 2014
3 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/c3b5065b5fba/
Changeset: c3b5065b5fba
Branch: yt
User: jzuhone
Date: 2014-11-03 18:46:05+00:00
Summary: Adding has_equivalent method.
Affected #: 2 files
diff -r 7953b98ce611c9f58896f56f85440bd98e13aec5 -r c3b5065b5fbaad46920e867edce8b3ab4a839ec8 yt/units/equivalencies.py
--- a/yt/units/equivalencies.py
+++ b/yt/units/equivalencies.py
@@ -11,7 +11,6 @@
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
-import yt.utilities.physical_constants as pc
from yt.units.dimensions import temperature, mass, energy, length, rate, \
velocity, dimensionless, density, number_density, flux
from yt.extern.six import add_metaclass
@@ -29,15 +28,19 @@
class Equivalence(object):
_skip_add = False
+ def __init__(self):
+ import yt.utilities.physical_constants as pc
+ self.pc = pc
+
class NumberDensityEquivalence(Equivalence):
_type_name = "number_density"
dims = (density,number_density,)
def convert(self, x, new_dims, mu=0.6):
if new_dims == number_density:
- return x/(mu*pc.mh)
+ return x/(mu*self.pc.mh)
elif new_dims == density:
- return x*mu*pc.mh
+ return x*mu*self.pc.mh
def __str__(self):
return "number density: density <-> number density"
@@ -48,9 +51,9 @@
def convert(self, x, new_dims):
if new_dims == energy:
- return pc.kboltz*x
+ return self.pc.kboltz*x
elif new_dims == temperature:
- return x/pc.kboltz
+ return x/self.pc.kboltz
def __str__(self):
return "thermal: temperature <-> energy"
@@ -61,9 +64,9 @@
def convert(self, x, new_dims):
if new_dims == energy:
- return x*pc.clight*pc.clight
+ return x*self.pc.clight*self.pc.clight
elif new_dims == mass:
- return x/(pc.clight*pc.clight)
+ return x/(self.pc.clight*self.pc.clight)
def __str__(self):
return "mass_energy: mass <-> energy"
@@ -75,20 +78,20 @@
def convert(self, x, new_dims):
if new_dims == energy:
if x.units.dimensions == length:
- nu = pc.clight/x
+ nu = self.pc.clight/x
elif x.units.dimensions == rate:
nu = x
- return pc.hcgs*nu
+ return self.pc.hcgs*nu
elif new_dims == length:
if x.units.dimensions == rate:
- return pc.clight/x
+ return self.pc.clight/x
elif x.units.dimensions == energy:
- return pc.hcgs*pc.clight/x
+ return self.pc.hcgs*self.pc.clight/x
elif new_dims == rate:
if x.units.dimensions == length:
- return pc.clight/x
+ return self.pc.clight/x
elif x.units.dimensions == energy:
- return x/pc.hcgs
+ return x/self.pc.hcgs
def __str__(self):
return "spectral: length <-> rate <-> energy"
@@ -100,14 +103,14 @@
def convert(self, x, new_dims, mu=0.6, gamma=5./3.):
if new_dims == velocity:
if x.units.dimensions == temperature:
- kT = pc.kboltz*x
+ kT = self.pc.kboltz*x
elif x.units.dimensions == energy:
kT = x
- return np.sqrt(gamma*kT/(mu*pc.mh))
+ return np.sqrt(gamma*kT/(mu*self.pc.mh))
else:
- kT = x*x*mu*pc.mh/gamma
+ kT = x*x*mu*self.pc.mh/gamma
if new_dims == temperature:
- return kT/pc.kboltz
+ return kT/self.pc.kboltz
else:
return kT
@@ -120,10 +123,10 @@
def convert(self, x, new_dims):
if new_dims == dimensionless:
- beta = x.in_cgs()/pc.clight
+ beta = x.in_cgs()/self.pc.clight
return 1./np.sqrt(1.-beta**2)
elif new_dims == velocity:
- return pc.clight*np.sqrt(1.-1./(x*x))
+ return self.pc.clight*np.sqrt(1.-1./(x*x))
def __str__(self):
return "lorentz: velocity <-> dimensionless"
@@ -134,9 +137,9 @@
def convert(self, x, new_dims):
if new_dims == length:
- return 2.*pc.G*x/(pc.clight*pc.clight)
+ return 2.*self.pc.G*x/(self.pc.clight*self.pc.clight)
elif new_dims == mass:
- return 0.5*x*pc.clight*pc.clight/pc.G
+ return 0.5*x*self.pc.clight*self.pc.clight/self.pc.G
def __str__(self):
return "schwarzschild: mass <-> length"
@@ -146,7 +149,7 @@
dims = (mass,length,)
def convert(self, x, new_dims):
- return pc.hcgs/(x*pc.clight)
+ return self.pc.hcgs/(x*self.pc.clight)
def __str__(self):
return "compton: mass <-> length"
@@ -157,9 +160,9 @@
def convert(self, x, new_dims):
if new_dims == flux:
- return pc.stefan_boltzmann_constant_cgs*x**4
+ return self.pc.stefan_boltzmann_constant_cgs*x**4
elif new_dims == temperature:
- return (x/pc.stefan_boltzmann_constant_cgs)**0.25
+ return (x/self.pc.stefan_boltzmann_constant_cgs)**0.25
def __str__(self):
return "effective_temperature: flux <-> temperature"
diff -r 7953b98ce611c9f58896f56f85440bd98e13aec5 -r c3b5065b5fbaad46920e867edce8b3ab4a839ec8 yt/units/yt_array.py
--- a/yt/units/yt_array.py
+++ b/yt/units/yt_array.py
@@ -39,6 +39,7 @@
from yt.utilities.on_demand_imports import _astropy
from sympy import Rational
from yt.units.unit_lookup_table import unit_prefixes, prefixable_units
+from yt.units.equivalencies import equivalence_registry
NULL_UNIT = Unit()
@@ -479,12 +480,10 @@
>>> a = yt.YTArray(1.0e7,"K")
>>> a.to_equivalent("keV", "thermal")
"""
- from equivalencies import equivalence_registry
- this_equiv = equivalence_registry[equiv]()
- old_dims = self.units.dimensions
- new_dims = YTQuantity(1.0, unit, registry=self.units.registry).units.dimensions
- if old_dims in this_equiv.dims and new_dims in this_equiv.dims:
- return this_equiv.convert(self, new_dims, **kwargs).in_units(unit)
+ unit_quan = YTQuantity(1.0, unit, registry=self.units.registry)
+ if self.has_equivalent(equiv) and unit_quan.has_equivalent(equiv):
+ this_equiv = equivalence_registry[equiv]()
+ return this_equiv.convert(self, unit_quan.units.dimensions, **kwargs).in_units(unit)
else:
raise YTInvalidUnitEquivalence(equiv, self.units, unit)
@@ -493,11 +492,22 @@
Lists the possible equivalencies associated with this YTArray or
YTQuantity.
"""
- from equivalencies import equivalence_registry
for k,v in equivalence_registry.items():
- if self.units.dimensions in v.dims:
+ if self.has_equivalent(k):
print v()
+ def has_equivalent(self, equiv):
+ """
+ Check to see if this YTArray or YTQuantity has an equivalent unit in
+ *equiv*.
+ """
+ try:
+ this_equiv = equivalence_registry[equiv]()
+ except KeyError:
+ raise KeyError("No such equivalence \"%s\"." % equiv)
+ old_dims = self.units.dimensions
+ return old_dims in this_equiv.dims
+
def ndarray_view(self):
"""
Returns a view into the array, but as an ndarray rather than ytarray.
https://bitbucket.org/yt_analysis/yt/commits/e9d1135f9e91/
Changeset: e9d1135f9e91
Branch: yt
User: jzuhone
Date: 2014-11-04 16:27:25+00:00
Summary: Adding has_equivalent to docs
Affected #: 1 file
diff -r c3b5065b5fbaad46920e867edce8b3ab4a839ec8 -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 doc/source/analyzing/units/6)_Unit_Equivalencies.ipynb
--- a/doc/source/analyzing/units/6)_Unit_Equivalencies.ipynb
+++ b/doc/source/analyzing/units/6)_Unit_Equivalencies.ipynb
@@ -1,7 +1,7 @@
{
"metadata": {
"name": "",
- "signature": "sha256:cee652d703dd3369d81ebc670882d3734f73d0274aab98823a784d8039355480"
+ "signature": "sha256:b62d83c168828afa81bcf0603bb37d3183f2a810258f25963254ffb24a0acd82"
},
"nbformat": 3,
"nbformat_minor": 0,
@@ -159,6 +159,24 @@
"cell_type": "markdown",
"metadata": {},
"source": [
+ "You can check if a `YTArray` has a given equivalence with `has_equivalent`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "print mp.has_equivalent(\"compton\")\n",
+ "print mp.has_equivalent(\"thermal\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
"To list the equivalencies available for a given `YTArray` or `YTQuantity`, use the `list_equivalencies` method:"
]
},
https://bitbucket.org/yt_analysis/yt/commits/4c67e2a36f2c/
Changeset: 4c67e2a36f2c
Branch: yt
User: jzuhone
Date: 2014-11-06 15:45:03+00:00
Summary: Merge
Affected #: 36 files
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 doc/source/analyzing/analysis_modules/PPVCube.ipynb
--- a/doc/source/analyzing/analysis_modules/PPVCube.ipynb
+++ b/doc/source/analyzing/analysis_modules/PPVCube.ipynb
@@ -1,7 +1,7 @@
{
"metadata": {
"name": "",
- "signature": "sha256:56a8d72735e3cc428ff04b241d4b2ce6f653019818c6fc7a4148840d99030c85"
+ "signature": "sha256:b83e125278c2e58da4d99ac9d2ca2a136d01f1094e1b83497925e0f9b9b056c2"
},
"nbformat": 3,
"nbformat_minor": 0,
@@ -32,7 +32,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "To demonstrate this functionality, we'll create a simple unigrid dataset from scratch of a rotating disk galaxy. We create a thin disk in the x-y midplane of the domain of three cells in height in either direction, and a radius of 10 kpc. The density and azimuthal velocity profiles of the disk as a function of radius will be given by the following functions:"
+ "To demonstrate this functionality, we'll create a simple unigrid dataset from scratch of a rotating disk. We create a thin disk in the x-y midplane of the domain of three cells in height in either direction, and a radius of 10 kpc. The density and azimuthal velocity profiles of the disk as a function of radius will be given by the following functions:"
]
},
{
@@ -84,7 +84,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "Second, we'll construct the data arrays for the density and the velocity of the disk. Since we have the tangential velocity profile, we have to use the polar coordinates we derived earlier to compute `velx` and `vely`. Everywhere outside the disk, all fields are set to zero. "
+ "Second, we'll construct the data arrays for the density, temperature, and velocity of the disk. Since we have the tangential velocity profile, we have to use the polar coordinates we derived earlier to compute `velx` and `vely`. Everywhere outside the disk, all fields are set to zero. "
]
},
{
@@ -93,12 +93,15 @@
"input": [
"dens = np.zeros((nx,ny,nz))\n",
"dens[:,:,nz/2-3:nz/2+3] = (r**alpha).reshape(nx,ny,1) # the density profile of the disk\n",
- "vel_theta = r/(1.+(r/r_0)**beta) # the azimuthal velocity profile of the disk\n",
+ "temp = np.zeros((nx,ny,nz))\n",
+ "temp[:,:,nz/2-3:nz/2+3] = 1.0e5 # Isothermal\n",
+ "vel_theta = 100.*r/(1.+(r/r_0)**beta) # the azimuthal velocity profile of the disk\n",
"velx = np.zeros((nx,ny,nz))\n",
"vely = np.zeros((nx,ny,nz))\n",
"velx[:,:,nz/2-3:nz/2+3] = (-vel_theta*np.sin(theta)).reshape(nx,ny,1) # convert polar to cartesian\n",
"vely[:,:,nz/2-3:nz/2+3] = (vel_theta*np.cos(theta)).reshape(nx,ny,1) # convert polar to cartesian\n",
"dens[r > R] = 0.0\n",
+ "temp[r > R] = 0.0\n",
"velx[r > R] = 0.0\n",
"vely[r > R] = 0.0"
],
@@ -119,6 +122,7 @@
"input": [
"data = {}\n",
"data[\"density\"] = (dens,\"g/cm**3\")\n",
+ "data[\"temperature\"] = (temp, \"K\")\n",
"data[\"velocity_x\"] = (velx, \"km/s\")\n",
"data[\"velocity_y\"] = (vely, \"km/s\")\n",
"data[\"velocity_z\"] = (np.zeros((nx,ny,nz)), \"km/s\") # zero velocity in the z-direction\n",
@@ -189,7 +193,7 @@
"cell_type": "code",
"collapsed": false,
"input": [
- "cube = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-1.5,1.5,\"km/s\"))"
+ "cube = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-150.,150.,\"km/s\"))"
],
"language": "python",
"metadata": {},
@@ -199,14 +203,33 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "Following this, we can now write this cube to a FITS file:"
+ "Following this, we can now write this cube to a FITS file. The x and y axes of the file can be in length units, which can be optionally specified by `length_unit`:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
- "cube.write_fits(\"cube.fits\", clobber=True, length_unit=(5.0,\"deg\"))"
+ "cube.write_fits(\"cube.fits\", clobber=True, length_unit=\"kpc\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Or one can use the `sky_scale` and `sky_center` keywords to set up the coordinates in RA and Dec:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "sky_scale = (1.0, \"arcsec/kpc\")\n",
+ "sky_center = (30., 45.) # RA, Dec in degrees\n",
+ "cube.write_fits(\"cube_sky.fits\", clobber=True, sky_scale=sky_scale, sky_center=sky_center)"
],
"language": "python",
"metadata": {},
@@ -223,7 +246,7 @@
"cell_type": "code",
"collapsed": false,
"input": [
- "ds = yt.load(\"cube.fits\")"
+ "ds_cube = yt.load(\"cube.fits\")"
],
"language": "python",
"metadata": {},
@@ -234,7 +257,7 @@
"collapsed": false,
"input": [
"# Specifying no center gives us the center slice\n",
- "slc = yt.SlicePlot(ds, \"z\", [\"density\"])\n",
+ "slc = yt.SlicePlot(ds_cube, \"z\", [\"density\"])\n",
"slc.show()"
],
"language": "python",
@@ -247,9 +270,9 @@
"input": [
"import yt.units as u\n",
"# Picking different velocities for the slices\n",
- "new_center = ds.domain_center\n",
- "new_center[2] = ds.spec2pixel(-1.0*u.km/u.s)\n",
- "slc = yt.SlicePlot(ds, \"z\", [\"density\"], center=new_center)\n",
+ "new_center = ds_cube.domain_center\n",
+ "new_center[2] = ds_cube.spec2pixel(-100.*u.km/u.s)\n",
+ "slc = yt.SlicePlot(ds_cube, \"z\", [\"density\"], center=new_center)\n",
"slc.show()"
],
"language": "python",
@@ -260,8 +283,8 @@
"cell_type": "code",
"collapsed": false,
"input": [
- "new_center[2] = ds.spec2pixel(0.7*u.km/u.s)\n",
- "slc = yt.SlicePlot(ds, \"z\", [\"density\"], center=new_center)\n",
+ "new_center[2] = ds_cube.spec2pixel(70.0*u.km/u.s)\n",
+ "slc = yt.SlicePlot(ds_cube, \"z\", [\"density\"], center=new_center)\n",
"slc.show()"
],
"language": "python",
@@ -272,8 +295,8 @@
"cell_type": "code",
"collapsed": false,
"input": [
- "new_center[2] = ds.spec2pixel(-0.3*u.km/u.s)\n",
- "slc = yt.SlicePlot(ds, \"z\", [\"density\"], center=new_center)\n",
+ "new_center[2] = ds_cube.spec2pixel(-30.0*u.km/u.s)\n",
+ "slc = yt.SlicePlot(ds_cube, \"z\", [\"density\"], center=new_center)\n",
"slc.show()"
],
"language": "python",
@@ -291,7 +314,7 @@
"cell_type": "code",
"collapsed": false,
"input": [
- "prj = yt.ProjectionPlot(ds, \"z\", [\"density\"], method=\"sum\")\n",
+ "prj = yt.ProjectionPlot(ds_cube, \"z\", [\"density\"], method=\"sum\")\n",
"prj.set_log(\"density\", True)\n",
"prj.set_zlim(\"density\", 1.0e-3, 0.2)\n",
"prj.show()"
@@ -299,9 +322,100 @@
"language": "python",
"metadata": {},
"outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `thermal_broad` keyword allows one to simulate thermal line broadening based on the temperature, and the `atomic_weight` argument is used to specify the atomic weight of the particle that is doing the emitting."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "cube2 = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-150,150,\"km/s\"), thermal_broad=True, \n",
+ " atomic_weight=12.0)\n",
+ "cube2.write_fits(\"cube2.fits\", clobber=True, length_unit=\"kpc\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Taking a slice of this cube shows:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "ds_cube2 = yt.load(\"cube2.fits\")\n",
+ "new_center = ds_cube2.domain_center\n",
+ "new_center[2] = ds_cube2.spec2pixel(70.0*u.km/u.s)\n",
+ "slc = yt.SlicePlot(ds_cube2, \"z\", [\"density\"], center=new_center)\n",
+ "slc.show()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "new_center[2] = ds_cube2.spec2pixel(-100.*u.km/u.s)\n",
+ "slc = yt.SlicePlot(ds_cube2, \"z\", [\"density\"], center=new_center)\n",
+ "slc.show()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where we can see the emission has been smeared into this velocity slice from neighboring slices due to the thermal broadening. \n",
+ "\n",
+ "Finally, the \"velocity\" or \"spectral\" axis of the cube can be changed to a different unit, such as wavelength, frequency, or energy: "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "print cube2.vbins[0], cube2.vbins[-1]\n",
+ "cube2.transform_spectral_axis(400.0,\"nm\")\n",
+ "print cube2.vbins[0], cube2.vbins[-1]"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "If a FITS file is now written from the cube, the spectral axis will be in the new units. To reset the spectral axis back to the original velocity units:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "cube2.reset_spectral_axis()\n",
+ "print cube2.vbins[0], cube2.vbins[-1]"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
}
],
"metadata": {}
}
]
-}
+}
\ No newline at end of file
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 doc/source/analyzing/objects.rst
--- a/doc/source/analyzing/objects.rst
+++ b/doc/source/analyzing/objects.rst
@@ -96,7 +96,7 @@
**Point**
| Class :class:`~yt.data_objects.selection_data_containers.YTPointBase`
- | Usage: ``point(coord, ds=None, field_parameters=None)``
+ | Usage: ``point(coord, ds=None, field_parameters=None, data_source=None)``
| A point defined by a single cell at specified coordinates.
1D Objects
@@ -104,14 +104,14 @@
**Ray (Axis-Aligned)**
| Class :class:`~yt.data_objects.selection_data_containers.YTOrthoRayBase`
- | Usage: ``ortho_ray(axis, coord, ds=None, field_parameters=None)``
+ | Usage: ``ortho_ray(axis, coord, ds=None, field_parameters=None, data_source=None)``
| A line (of data cells) stretching through the full domain
aligned with one of the x,y,z axes. Defined by an axis and a point
to be intersected.
**Ray (Arbitrarily-Aligned)**
| Class :class:`~yt.data_objects.selection_data_containers.YTRayBase`
- | Usage: ``ray(start_coord, end_coord, ds=None, field_parameters=None)``
+ | Usage: ``ray(start_coord, end_coord, ds=None, field_parameters=None, data_source=None)``
| A line (of data cells) defined by arbitrary start and end coordinates.
2D Objects
@@ -119,13 +119,13 @@
**Slice (Axis-Aligned)**
| Class :class:`~yt.data_objects.selection_data_containers.YTSliceBase`
- | Usage: ``slice(axis, coord, center=None, ds=None, field_parameters=None)``
+ | Usage: ``slice(axis, coord, center=None, ds=None, field_parameters=None, data_source=None)``
| A plane normal to one of the axes and intersecting a particular
coordinate.
**Slice (Arbitrarily-Aligned)**
| Class :class:`~yt.data_objects.selection_data_containers.YTCuttingPlaneBase`
- | Usage: ``cutting(normal, coord, north_vector=None, ds=None, field_parameters=None)``
+ | Usage: ``cutting(normal, coord, north_vector=None, ds=None, field_parameters=None, data_source=None)``
| A plane normal to a specified vector and intersecting a particular
coordinate.
@@ -141,8 +141,8 @@
**Box Region**
| Class :class:`~yt.data_objects.selection_data_containers.YTRegionBase`
- | Usage: ``region(center, left_edge, right_edge, fields=None, ds=None, field_parameters=None)``
- | Alternatively: ``box(left_edge, right_edge, fields=None, ds=None, field_parameters=None)``
+ | Usage: ``region(center, left_edge, right_edge, fields=None, ds=None, field_parameters=None, data_source=None)``
+ | Alternatively: ``box(left_edge, right_edge, fields=None, ds=None, field_parameters=None, data_source=None)``
| A box-like region aligned with the grid axis orientation. It is
defined by a left_edge, a right_edge, and a center. The left_edge
and right_edge are the minimum and maximum bounds in the three axes
@@ -152,14 +152,14 @@
**Disk/Cylinder**
| Class: :class:`~yt.data_objects.selection_data_containers.YTDiskBase`
- | Usage: ``disk(center, normal, radius, height, fields=None, ds=None, field_parameters=None)``
+ | Usage: ``disk(center, normal, radius, height, fields=None, ds=None, field_parameters=None, data_source=None)``
| A cylinder defined by a point at the center of one of the circular bases,
a normal vector to it defining the orientation of the length of the
cylinder, and radius and height values for the cylinder's dimensions.
**Ellipsoid**
| Class :class:`~yt.data_objects.selection_data_containers.YTEllipsoidBase`
- | Usage: ``ellipsoid(center, semi_major_axis_length, semi_medium_axis_length, semi_minor_axis_length, semi_major_vector, tilt, fields=None, ds=None, field_parameters=None)``
+ | Usage: ``ellipsoid(center, semi_major_axis_length, semi_medium_axis_length, semi_minor_axis_length, semi_major_vector, tilt, fields=None, ds=None, field_parameters=None, data_source=None)``
| An ellipsoid with axis magnitudes set by semi_major_axis_length,
semi_medium_axis_length, and semi_minor_axis_length. semi_major_vector
sets the direction of the semi_major_axis. tilt defines the orientation
@@ -167,7 +167,7 @@
**Sphere**
| Class :class:`~yt.data_objects.selection_data_containers.YTSphereBase`
- | Usage: ``sphere(center, radius, ds=None, field_parameters=None)``
+ | Usage: ``sphere(center, radius, ds=None, field_parameters=None, data_source=None)``
| A sphere defined by a central coordinate and a radius.
@@ -176,6 +176,12 @@
See also the section on :ref:`filtering-data`.
+**Intersecting Regions**
+ | Most Region objects provide a data_source parameter, which allows you to subselect
+ | one region from another (in the coordinate system of the DataSet). Note, this can
+ | easily lead to empty data for non-intersecting regions.
+ | Usage: ``slice(axis, coord, ds, data_source=sph)``
+
**Boolean Regions**
| **Note: not yet implemented in yt 3.0**
| Usage: ``boolean()``
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 doc/source/visualizing/plots.rst
--- a/doc/source/visualizing/plots.rst
+++ b/doc/source/visualizing/plots.rst
@@ -124,11 +124,16 @@
slice. To instead use the coordinates as defined in the dataset, use
the optional argument: ``origin="native"``
-If supplied
-without units, the center is assumed by in code units. Optionally, you can
-supply 'c' or 'm' for the center. These two choices will center the plot on the
-center of the simulation box and the coordinate of the maximum density cell,
-respectively.
+If supplied without units, the center is assumed by in code units. There are also
+the following alternative options for the `center` keyword:
+
+* `"center"`, "c"`: the domain center
+* `"max"`, `"m"`: the position of the maximum density
+* `("min", field)`: the position of the minimum of `field`
+* `("max", field)`: the position of the maximum of `field`
+
+where for the last two objects any spatial field, such as `"density"`, `"velocity_z"`,
+etc., may be used, e.g. `center=("min","temperature")`
Here is an example that combines all of the options we just discussed.
@@ -278,7 +283,8 @@
straight summation of the field along the given axis. The units of the
projected field will be the same as those of the unprojected field. This
method is typically only useful for datasets such as 3D FITS cubes where
- the third axis of the dataset is something like velocity or frequency.
+ the third axis of the dataset is something like velocity or frequency, and
+ should _only_ be used with fixed-resolution grid-based datasets.
.. _off-axis-projections:
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/analysis_modules/halo_finding/halo_objects.py
--- a/yt/analysis_modules/halo_finding/halo_objects.py
+++ b/yt/analysis_modules/halo_finding/halo_objects.py
@@ -598,7 +598,7 @@
# 1. From saved_fields, e.g. we've already got it.
# 2. From the halo binary files off disk.
# 3. Use the unique particle indexes of the halo to select a missing
- # field from an AMR Sphere.
+ # field from a Sphere.
if key in self._saved_fields:
# We've already got it.
return self._saved_fields[key]
@@ -675,7 +675,7 @@
Returns
-------
- ellipsoid : `yt.data_objects.api.AMREllipsoidBase`
+ ellipsoid : `yt.data_objects.data_containers.YTEllipsoidBase`
The ellipsoidal data object.
Examples
@@ -754,7 +754,7 @@
# 1. From saved_fields, e.g. we've already got it.
# 2. From the halo h5 files off disk.
# 3. Use the unique particle indexes of the halo to select a missing
- # field from an AMR Sphere.
+ # field from a Sphere.
if key in self._saved_fields:
# We've already got it.
return self._saved_fields[key]
@@ -868,7 +868,7 @@
Returns
-------
- ellipsoid : `yt.data_objects.api.AMREllipsoidBase`
+ ellipsoid : `yt.data_objects.data_containers.YTEllipsoidBase`
The ellipsoidal data object.
Examples
@@ -1670,7 +1670,7 @@
----------
ds : `Dataset`
The dataset on which halo finding will be conducted.
- subvolume : `yt.data_objects.api.AMRData`, optional
+ subvolume : `yt.data_objects.data_containers.YTSelectionContainer`, optional
A region over which HOP will be run, which can be used to run HOP
on a subvolume of the full volume. Default = None, which defaults
to the full volume automatically.
@@ -1772,7 +1772,7 @@
----------
ds : `Dataset`
The dataset on which halo finding will be conducted.
- subvolume : `yt.data_objects.api.AMRData`, optional
+ subvolume : `yt.data_objects.data_containers.YTSelectionContainer`, optional
A region over which HOP will be run, which can be used to run HOP
on a subvolume of the full volume. Default = None, which defaults
to the full volume automatically.
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/analysis_modules/photon_simulator/photon_simulator.py
--- a/yt/analysis_modules/photon_simulator/photon_simulator.py
+++ b/yt/analysis_modules/photon_simulator/photon_simulator.py
@@ -150,7 +150,7 @@
Parameters
----------
- data_source : `yt.data_objects.api.AMRData`
+ data_source : `yt.data_objects.data_containers.YTSelectionContainer`
The data source from which the photons will be generated.
redshift : float
The cosmological redshift for the photons.
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/analysis_modules/ppv_cube/ppv_cube.py
--- a/yt/analysis_modules/ppv_cube/ppv_cube.py
+++ b/yt/analysis_modules/ppv_cube/ppv_cube.py
@@ -13,40 +13,81 @@
import numpy as np
from yt.utilities.on_demand_imports import _astropy
from yt.utilities.orientation import Orientation
-from yt.utilities.fits_image import FITSImageBuffer
+from yt.utilities.fits_image import FITSImageBuffer, sanitize_fits_unit
from yt.visualization.volume_rendering.camera import off_axis_projection
from yt.funcs import get_pbar
+from yt.utilities.physical_constants import clight, mh
+import yt.units.dimensions as ytdims
+from yt.units.yt_array import YTQuantity
+from yt.funcs import iterable
+from yt.utilities.parallel_tools.parallel_analysis_interface import \
+ parallel_root_only, parallel_objects
+import re
+import ppv_utils
+from yt.funcs import is_root
-def create_vlos(z_hat):
- def _v_los(field, data):
- vz = data["velocity_x"]*z_hat[0] + \
- data["velocity_y"]*z_hat[1] + \
- data["velocity_z"]*z_hat[2]
- return -vz
+def create_vlos(normal):
+ if isinstance(normal, basestring):
+ def _v_los(field, data):
+ return -data["velocity_%s" % normal]
+ else:
+ orient = Orientation(normal)
+ los_vec = orient.unit_vectors[2]
+ def _v_los(field, data):
+ vz = data["velocity_x"]*los_vec[0] + \
+ data["velocity_y"]*los_vec[1] + \
+ data["velocity_z"]*los_vec[2]
+ return -vz
return _v_los
+fits_info = {"velocity":("m/s","VELOCITY","v"),
+ "frequency":("Hz","FREQUENCY","f"),
+ "energy":("eV","ENERGY","E"),
+ "wavelength":("angstrom","WAVELENG","lambda")}
+
class PPVCube(object):
- def __init__(self, ds, normal, field, width=(1.0,"unitary"),
- dims=(100,100,100), velocity_bounds=None):
+ def __init__(self, ds, normal, field, center="c", width=(1.0,"unitary"),
+ dims=(100,100,100), velocity_bounds=None, thermal_broad=False,
+ atomic_weight=56., method="integrate"):
r""" Initialize a PPVCube object.
Parameters
----------
ds : dataset
The dataset.
- normal : array_like
- The normal vector along with to make the projections.
+ normal : array_like or string
+ The normal vector along with to make the projections. If an array, it
+ will be normalized. If a string, it will be assumed to be along one of the
+ principal axes of the domain ("x","y", or "z").
field : string
The field to project.
- width : float or tuple, optional
- The width of the projection in length units. Specify a float
- for code_length units or a tuple (value, units).
+ 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 : float, tuple, or YTQuantity.
+ The width of the projection. A float will assume the width is in code units.
+ A (value, unit) tuple or YTQuantity allows for the units of the width to be
+ specified.
dims : tuple, optional
A 3-tuple of dimensions (nx,ny,nv) for the cube.
velocity_bounds : tuple, optional
A 3-tuple of (vmin, vmax, units) for the velocity bounds to
integrate over. If None, the largest velocity of the
dataset will be used, e.g. velocity_bounds = (-v.max(), v.max())
+ atomic_weight : float, optional
+ Set this value to the atomic weight of the particle that is emitting the line
+ if *thermal_broad* is True. Defaults to 56 (Fe).
+ method : string, optional
+ Set the projection method to be used.
+ "integrate" : line of sight integration over the line element.
+ "sum" : straight summation over the line of sight.
Examples
--------
@@ -55,21 +96,22 @@
>>> cube = PPVCube(ds, L, "density", width=(10.,"kpc"),
... velocity_bounds=(-5.,4.,"km/s"))
"""
+
self.ds = ds
self.field = field
self.width = width
+ self.particle_mass = atomic_weight*mh
+ self.thermal_broad = thermal_broad
+
+ self.center = ds.coordinates.sanitize_center(center, normal)[0]
self.nx = dims[0]
self.ny = dims[1]
self.nv = dims[2]
- normal = np.array(normal)
- normal /= np.sqrt(np.dot(normal, normal))
- vecs = np.identity(3)
- t = np.cross(normal, vecs).sum(axis=1)
- ax = t.argmax()
- north = np.cross(normal, vecs[ax,:]).ravel()
- orient = Orientation(normal, north_vector=north)
+ if method not in ["integrate","sum"]:
+ raise RuntimeError("Only the 'integrate' and 'sum' projection +"
+ "methods are supported in PPVCube.")
dd = ds.all_data()
@@ -85,27 +127,104 @@
ds.quan(velocity_bounds[1], velocity_bounds[2]))
self.vbins = np.linspace(self.v_bnd[0], self.v_bnd[1], num=self.nv+1)
+ self._vbins = self.vbins.copy()
self.vmid = 0.5*(self.vbins[1:]+self.vbins[:-1])
- self.dv = (self.v_bnd[1]-self.v_bnd[0])/self.nv
+ self.vmid_cgs = self.vmid.in_cgs().v
+ self.dv = self.vbins[1]-self.vbins[0]
+ self.dv_cgs = self.dv.in_cgs().v
- _vlos = create_vlos(orient.unit_vectors[2])
- ds.field_info.add_field(("gas","v_los"), function=_vlos, units="cm/s")
+ self.current_v = 0.0
- self.data = ds.arr(np.zeros((self.nx,self.ny,self.nv)), self.field_units)
+ _vlos = create_vlos(normal)
+ self.ds.add_field(("gas","v_los"), function=_vlos, units="cm/s")
+
+ _intensity = self.create_intensity()
+ self.ds.add_field(("gas","intensity"), function=_intensity, units=self.field_units)
+
+ if method == "integrate":
+ self.proj_units = str(ds.quan(1.0, self.field_units+"*cm").units)
+ elif method == "sum":
+ self.proj_units = self.field_units
+
+ self.data = ds.arr(np.zeros((self.nx,self.ny,self.nv)), self.proj_units)
+ storage = {}
pbar = get_pbar("Generating cube.", self.nv)
- for i in xrange(self.nv):
- _intensity = self._create_intensity(i)
- ds.add_field(("gas","intensity"), function=_intensity, units=self.field_units)
- prj = off_axis_projection(ds, ds.domain_center, normal, width,
- (self.nx, self.ny), "intensity")
- self.data[:,:,i] = prj[:,:]
- ds.field_info.pop(("gas","intensity"))
+ for sto, i in parallel_objects(xrange(self.nv), storage=storage):
+ self.current_v = self.vmid_cgs[i]
+ if isinstance(normal, basestring):
+ prj = ds.proj("intensity", ds.coordinates.axis_id[normal], method=method)
+ buf = prj.to_frb(width, self.nx, center=self.center)["intensity"]
+ else:
+ buf = off_axis_projection(ds, self.center, normal, width,
+ (self.nx, self.ny), "intensity",
+ no_ghost=True, method=method)[::-1]
+ sto.result_id = i
+ sto.result = buf
pbar.update(i)
-
pbar.finish()
- def write_fits(self, filename, clobber=True, length_unit=(10.0, "kpc"),
- sky_center=(30.,45.)):
+ self.data = ds.arr(np.zeros((self.nx,self.ny,self.nv)), self.proj_units)
+ if is_root():
+ for i, buf in sorted(storage.items()):
+ self.data[:,:,i] = buf[:,:]
+
+ self.axis_type = "velocity"
+
+ # Now fix the width
+ if iterable(self.width):
+ self.width = ds.quan(self.width[0], self.width[1])
+ elif isinstance(self.width, YTQuantity):
+ self.width = width
+ else:
+ self.width = ds.quan(self.width, "code_length")
+
+ def create_intensity(self):
+ def _intensity(field, data):
+ v = self.current_v-data["v_los"].v
+ T = data["temperature"].v
+ w = ppv_utils.compute_weight(self.thermal_broad, self.dv_cgs,
+ self.particle_mass, v.flatten(), T.flatten())
+ w[np.isnan(w)] = 0.0
+ return data[self.field]*w.reshape(v.shape)
+ return _intensity
+
+ def transform_spectral_axis(self, rest_value, units):
+ """
+ Change the units of the spectral axis to some equivalent unit, such
+ as energy, wavelength, or frequency, by providing a *rest_value* and the
+ *units* of the new spectral axis. This corresponds to the Doppler-shifting
+ of lines due to gas motions and thermal broadening.
+ """
+ if self.axis_type != "velocity":
+ self.reset_spectral_axis()
+ x0 = self.ds.quan(rest_value, units)
+ if x0.units.dimensions == ytdims.rate or x0.units.dimensions == ytdims.energy:
+ self.vbins = x0*(1.-self.vbins.in_cgs()/clight)
+ elif x0.units.dimensions == ytdims.length:
+ self.vbins = x0/(1.-self.vbins.in_cgs()/clight)
+ self.vmid = 0.5*(self.vbins[1:]+self.vbins[:-1])
+ self.dv = self.vbins[1]-self.vbins[0]
+ dims = self.dv.units.dimensions
+ if dims == ytdims.rate:
+ self.axis_type = "frequency"
+ elif dims == ytdims.length:
+ self.axis_type = "wavelength"
+ elif dims == ytdims.energy:
+ self.axis_type = "energy"
+ elif dims == ytdims.velocity:
+ self.axis_type = "velocity"
+
+ def reset_spectral_axis(self):
+ """
+ Reset the spectral axis to the original velocity range and units.
+ """
+ self.vbins = self._vbins.copy()
+ self.vmid = 0.5*(self.vbins[1:]+self.vbins[:-1])
+ self.dv = self.vbins[1]-self.vbins[0]
+
+ @parallel_root_only
+ def write_fits(self, filename, clobber=True, length_unit=None,
+ sky_scale=None, sky_center=None):
r""" Write the PPVCube to a FITS file.
Parameters
@@ -114,51 +233,71 @@
The name of the file to write.
clobber : boolean
Whether or not to clobber an existing file with the same name.
- length_unit : tuple, optional
- The length that corresponds to the width of the projection in
- (value, unit) form. Accepts a length unit or 'deg'.
+ length_unit : string
+ The units to convert the coordinates to in the file.
+ sky_scale : tuple or YTQuantity
+ Conversion between an angle unit and a length unit, if sky
+ coordinates are desired.
+ Examples: (1.0, "arcsec/kpc"), YTQuantity(0.001, "deg/kpc")
sky_center : tuple, optional
The (RA, Dec) coordinate in degrees of the central pixel if
- *length_unit* is 'deg'.
+ *sky_scale* has been specified.
Examples
--------
- >>> cube.write_fits("my_cube.fits", clobber=False, length_unit=(5,"deg"))
+ >>> cube.write_fits("my_cube.fits", clobber=False, sky_scale=(1.0,"arcsec/kpc"))
"""
- if length_unit[1] == "deg":
- center = sky_center
- types = ["RA---SIN","DEC--SIN"]
+ if sky_scale is None:
+ center = (0.0,0.0)
+ types = ["LINEAR","LINEAR"]
else:
- center = [0.0,0.0]
- types = ["LINEAR","LINEAR"]
+ if iterable(sky_scale):
+ sky_scale = self.ds.quan(sky_scale[0], sky_scale[1])
+ if sky_center is None:
+ center = (30.,45.)
+ else:
+ center = sky_center
+ types = ["RA---TAN","DEC--TAN"]
- v_center = 0.5*(self.v_bnd[0]+self.v_bnd[1]).in_units("m/s").value
+ vunit = fits_info[self.axis_type][0]
+ vtype = fits_info[self.axis_type][1]
- dx = length_unit[0]/self.nx
- dy = length_unit[0]/self.ny
- dv = self.dv.in_units("m/s").value
+ v_center = 0.5*(self.vbins[0]+self.vbins[-1]).in_units(vunit).value
- if length_unit[1] == "deg":
+ if sky_scale:
+ dx = (self.width*sky_scale).in_units("deg").v/self.nx
+ units = "deg"
+ else:
+ if length_unit is None:
+ units = str(self.ds.get_smallest_appropriate_unit(self.width))
+ else:
+ units = length_unit
+ units = sanitize_fits_unit(units)
+ dx = self.width.in_units(units).v/self.nx
+ dy = dx
+ dv = self.dv.in_units(vunit).v
+
+ if sky_scale:
dx *= -1.
w = _astropy.pywcs.WCS(naxis=3)
w.wcs.crpix = [0.5*(self.nx+1), 0.5*(self.ny+1), 0.5*(self.nv+1)]
w.wcs.cdelt = [dx,dy,dv]
- w.wcs.crval = [center[0], center[1], v_center]
- w.wcs.cunit = [length_unit[1],length_unit[1],"m/s"]
- w.wcs.ctype = [types[0],types[1],"VELO-LSR"]
+ w.wcs.crval = [center[0],center[1],v_center]
+ w.wcs.cunit = [units,units,vunit]
+ w.wcs.ctype = [types[0],types[1],vtype]
fib = FITSImageBuffer(self.data.transpose(), fields=self.field, wcs=w)
- fib[0].header["bunit"] = self.field_units
+ fib[0].header["bunit"] = re.sub('()', '', str(self.proj_units))
fib[0].header["btype"] = self.field
fib.writeto(filename, clobber=clobber)
- def _create_intensity(self, i):
- def _intensity(field, data):
- vlos = data["v_los"]
- w = np.abs(vlos-self.vmid[i])/self.dv.in_units(vlos.units)
- w = 1.-w
- w[w < 0.0] = 0.0
- return data[self.field]*w
- return _intensity
+ def __repr__(self):
+ return "PPVCube [%d %d %d] (%s < %s < %s)" % (self.nx, self.ny, self.nv,
+ self.vbins[0],
+ fits_info[self.axis_type][2],
+ self.vbins[-1])
+
+ def __getitem__(self, item):
+ return self.data[item]
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/analysis_modules/ppv_cube/ppv_utils.pyx
--- /dev/null
+++ b/yt/analysis_modules/ppv_cube/ppv_utils.pyx
@@ -0,0 +1,40 @@
+import numpy as np
+cimport numpy as np
+cimport cython
+from yt.utilities.physical_constants import kboltz
+
+cdef extern from "math.h":
+ double exp(double x) nogil
+ double fabs(double x) nogil
+ double sqrt(double x) nogil
+
+cdef double kb = kboltz.v
+cdef double pi = np.pi
+
+ at cython.cdivision(True)
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+def compute_weight(np.uint8_t thermal_broad,
+ double dv,
+ double m_part,
+ np.ndarray[np.float64_t, ndim=1] v,
+ np.ndarray[np.float64_t, ndim=1] T):
+
+ cdef int i, n
+ cdef double v2_th, x
+ cdef np.ndarray[np.float64_t, ndim=1] w
+
+ n = v.shape[0]
+ w = np.zeros(n)
+
+ for i in range(n):
+ if thermal_broad:
+ if T[i] > 0.0:
+ v2_th = 2.*kb*T[i]/m_part
+ w[i] = dv*exp(-v[i]*v[i]/v2_th)/sqrt(v2_th*pi)
+ else:
+ x = 1.-fabs(v[i])/dv
+ if x > 0.0:
+ w[i] = x
+
+ return w
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/analysis_modules/ppv_cube/setup.py
--- a/yt/analysis_modules/ppv_cube/setup.py
+++ b/yt/analysis_modules/ppv_cube/setup.py
@@ -8,7 +8,10 @@
def configuration(parent_package='', top_path=None):
from numpy.distutils.misc_util import Configuration
config = Configuration('ppv_cube', parent_package, top_path)
- #config.add_subpackage("tests")
+ config.add_extension("ppv_utils",
+ ["yt/analysis_modules/ppv_cube/ppv_utils.pyx"],
+ libraries=["m"])
+ config.add_subpackage("tests")
config.make_config_py() # installs __config__.py
#config.make_svn_version_py()
return config
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/analysis_modules/ppv_cube/tests/test_ppv.py
--- /dev/null
+++ b/yt/analysis_modules/ppv_cube/tests/test_ppv.py
@@ -0,0 +1,62 @@
+"""
+Unit test the sunyaev_zeldovich analysis module.
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+from yt.frontends.stream.api import load_uniform_grid
+from yt.analysis_modules.ppv_cube.api import PPVCube
+import yt.units as u
+from yt.utilities.physical_constants import kboltz, mh, clight
+import numpy as np
+from yt.testing import *
+
+def setup():
+ """Test specific setup."""
+ from yt.config import ytcfg
+ ytcfg["yt", "__withintesting"] = "True"
+
+def test_ppv():
+
+ np.random.seed(seed=0x4d3d3d3)
+
+ dims = (8,8,1024)
+ v_shift = 1.0e7*u.cm/u.s
+ sigma_v = 2.0e7*u.cm/u.s
+ T_0 = 1.0e8*u.Kelvin
+ data = {"density":(np.ones(dims),"g/cm**3"),
+ "temperature":(T_0.v*np.ones(dims), "K"),
+ "velocity_x":(np.zeros(dims),"cm/s"),
+ "velocity_y":(np.zeros(dims),"cm/s"),
+ "velocity_z":(np.random.normal(loc=v_shift.v,scale=sigma_v.v,size=dims), "cm/s")}
+
+ ds = load_uniform_grid(data, dims)
+
+ cube = PPVCube(ds, "z", "density", dims=dims,
+ velocity_bounds=(-300., 300., "km/s"),
+ thermal_broad=True)
+
+ dv = cube.dv
+ v_th = np.sqrt(2.*kboltz*T_0/(56.*mh) + 2.*sigma_v**2).in_units("km/s")
+ a = cube.data.mean(axis=(0,1)).v
+ b = dv*np.exp(-((cube.vmid+v_shift)/v_th)**2)/(np.sqrt(np.pi)*v_th)
+
+ yield assert_allclose, a, b, 1.0e-2
+
+ E_0 = 6.8*u.keV
+
+ cube.transform_spectral_axis(E_0.v, str(E_0.units))
+
+ dE = -cube.dv
+ delta_E = E_0*v_th.in_cgs()/clight
+ E_shift = E_0*(1.+v_shift/clight)
+
+ c = dE*np.exp(-((cube.vmid-E_shift)/delta_E)**2)/(np.sqrt(np.pi)*delta_E)
+
+ yield assert_allclose, a, c, 1.0e-2
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/analysis_modules/sunyaev_zeldovich/projection.py
--- a/yt/analysis_modules/sunyaev_zeldovich/projection.py
+++ b/yt/analysis_modules/sunyaev_zeldovich/projection.py
@@ -25,6 +25,7 @@
from yt.utilities.parallel_tools.parallel_analysis_interface import \
communication_system, parallel_root_only
from yt import units
+from yt.utilities.on_demand_imports import _astropy
import numpy as np
@@ -114,13 +115,22 @@
----------
axis : integer or string
The axis of the simulation domain along which to make the SZprojection.
- center : array_like or string, optional
- The center of the projection.
- width : float or tuple
- The width of the projection.
+ 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 : float, tuple, or YTQuantity.
+ The width of the projection. A float will assume the width is in code units.
+ A (value, unit) tuple or YTQuantity allows for the units of the width to be specified.
nx : integer, optional
The dimensions on a side of the projection image.
- source : yt.data_objects.api.AMRData, optional
+ source : yt.data_objects.data_containers.YTSelectionContainer, optional
If specified, this will be the data source used for selecting regions to project.
Examples
@@ -128,13 +138,8 @@
>>> szprj.on_axis("y", center="max", width=(1.0, "Mpc"), source=my_sphere)
"""
axis = fix_axis(axis, self.ds)
-
- if center == "c":
- ctr = self.ds.domain_center
- elif center == "max":
- v, ctr = self.ds.find_max("density")
- else:
- ctr = center
+ ctr, dctr = self.ds.coordinates.sanitize_center(center, axis)
+ width = self.ds.coordinates.sanitize_width(axis, width, None)
L = np.zeros((3))
L[axis] = 1.0
@@ -143,7 +148,7 @@
self.ds.add_field(("gas","beta_par"), function=beta_par, units="g/cm**3")
setup_sunyaev_zeldovich_fields(self.ds)
proj = self.ds.proj("density", axis, center=ctr, data_source=source)
- frb = proj.to_frb(width, nx)
+ frb = proj.to_frb(width[0], nx, height=width[1])
dens = frb["density"]
Te = frb["t_sz"]/dens
bpar = frb["beta_par"]/dens
@@ -176,13 +181,22 @@
----------
L : array_like
The normal vector of the projection.
- center : array_like or string, optional
- The center of the projection.
- width : float or tuple
- The width of the projection.
+ 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 : float, tuple, or YTQuantity.
+ The width of the projection. A float will assume the width is in code units.
+ A (value, unit) tuple or YTQuantity allows for the units of the width to be specified.
nx : integer, optional
The dimensions on a side of the projection image.
- source : yt.data_objects.api.AMRData, optional
+ source : yt.data_objects.data_containers.YTSelectionContainer, optional
If specified, this will be the data source used for selecting regions to project.
Currently unsupported in yt 2.x.
@@ -197,12 +211,7 @@
w = width.in_units("code_length").value
else:
w = width
- if center == "c":
- ctr = self.ds.domain_center
- elif center == "max":
- v, ctr = self.ds.find_max("density")
- else:
- ctr = center
+ ctr, dctr = self.ds.coordinates.sanitize_center(center, L)
if source is not None:
mylog.error("Source argument is not currently supported for off-axis S-Z projections.")
@@ -279,8 +288,7 @@
self.data["TeSZ"] = self.ds.arr(Te, "keV")
@parallel_root_only
- def write_fits(self, filename, units="kpc", sky_center=None, sky_scale=None,
- time_units="Gyr", clobber=True):
+ def write_fits(self, filename, sky_scale=None, sky_center=None, clobber=True):
r""" Export images to a FITS file. Writes the SZ distortion in all
specified frequencies as well as the mass-weighted temperature and the
optical depth. Distance units are in kpc, unless *sky_center*
@@ -290,12 +298,13 @@
----------
filename : string
The name of the FITS file to be written.
- sky_center : tuple of floats, optional
- The center of the observation in (RA, Dec) in degrees. Only used if
- converting to sky coordinates.
- sky_scale : float, optional
- Scale between degrees and kpc. Only used if
- converting to sky coordinates.
+ sky_scale : tuple or YTQuantity
+ Conversion between an angle unit and a length unit, if sky
+ coordinates are desired.
+ Examples: (1.0, "arcsec/kpc"), YTQuantity(0.001, "deg/kpc")
+ sky_center : tuple, optional
+ The (RA, Dec) coordinate in degrees of the central pixel if
+ *sky_scale* has been specified.
clobber : boolean, optional
If the file already exists, do we overwrite?
@@ -304,27 +313,43 @@
>>> # This example just writes out a FITS file with kpc coords
>>> szprj.write_fits("SZbullet.fits", clobber=False)
>>> # This example uses sky coords
- >>> sky_scale = 1./3600. # One arcsec per kpc
+ >>> sky_scale = (1., "arcsec/kpc") # One arcsec per kpc
>>> sky_center = (30., 45.) # In degrees
>>> szprj.write_fits("SZbullet.fits", sky_center=sky_center, sky_scale=sky_scale)
"""
+ from yt.utilities.fits_image import FITSImageBuffer, sanitize_fits_unit
- deltas = np.array([self.dx.in_units(units),
- self.dy.in_units(units)])
+ if sky_scale is None:
+ center = (0.0,0.0)
+ types = ["LINEAR","LINEAR"]
+ else:
+ if iterable(sky_scale):
+ sky_scale = self.ds.quan(sky_scale[0], sky_scale[1])
+ if sky_center is None:
+ center = (30.,45.)
+ else:
+ center = sky_center
+ types = ["RA---TAN","DEC--TAN"]
- if sky_center is None:
- center = [0.0]*2
- else:
- center = sky_center
+ units = self.ds.get_smallest_appropriate_unit(self.width)
+ units = sanitize_fits_unit(units)
+ # Hack because FITS is stupid and doesn't understand case
+ dx = self.dx.in_units(units)
+ if sky_scale:
+ dx = (dx*sky_scale).in_units("deg")
units = "deg"
- deltas *= sky_scale
- deltas[0] *= -1.
+ dy = dx
+ if sky_scale:
+ dx *= -1.
- from yt.utilities.fits_image import FITSImageBuffer
- fib = FITSImageBuffer(self.data, fields=self.data.keys(),
- center=center, units=units,
- scale=deltas)
- fib.update_all_headers("Time", float(self.ds.current_time.in_units(time_units).value))
+ w = _astropy.pywcs.WCS(naxis=2)
+ w.wcs.crpix = [0.5*(self.nx+1)]*2
+ w.wcs.cdelt = [dx.v,dy.v]
+ w.wcs.crval = center
+ w.wcs.cunit = [units]*2
+ w.wcs.ctype = types
+
+ fib = FITSImageBuffer(self.data, fields=self.data.keys(), wcs=w)
fib.writeto(filename, clobber=clobber)
@parallel_root_only
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -42,7 +42,7 @@
from yt.utilities.minimal_representation import \
MinimalProjectionData
from yt.utilities.parallel_tools.parallel_analysis_interface import \
- parallel_objects, parallel_root_only, ParallelAnalysisInterface
+ parallel_objects, parallel_root_only
from yt.units.unit_object import Unit
import yt.geometry.particle_deposit as particle_deposit
from yt.utilities.grid_data_format.writer import write_to_gdf
@@ -183,7 +183,7 @@
center : array_like, optional
The 'center' supplied to fields that use it. Note that this does
not have to have `coord` as one value. Strictly optional.
- data_source : `yt.data_objects.api.AMRData`, optional
+ data_source : `yt.data_objects.data_containers.YTSelectionContainer`, optional
If specified, this will be the data source used for selecting
regions to project.
method : string, optional
@@ -191,6 +191,8 @@
"integrate" : integration along the axis
"mip" : maximum intensity projection
"sum" : same as "integrate", except that we don't multiply by the path length
+ WARNING: The "sum" option should only be used for uniform resolution grid
+ datasets, as other datasets may result in unphysical images.
style : string, optional
The same as the method keyword. Deprecated as of version 3.0.2.
Please use method keyword instead.
@@ -833,7 +835,7 @@
new_fields.append(output_field)
level_state.fields = new_fields
-class YTSurfaceBase(YTSelectionContainer3D, ParallelAnalysisInterface):
+class YTSurfaceBase(YTSelectionContainer3D):
r"""This surface object identifies isocontours on a cell-by-cell basis,
with no consideration of global connectedness, and returns the vertices
of the Triangles in that isocontour.
@@ -850,7 +852,7 @@
Parameters
----------
- data_source : AMR3DDataObject
+ data_source : YTSelectionContainer
This is the object which will used as a source
surface_field : string
Any field that can be obtained in a data object. This is the field
@@ -886,7 +888,6 @@
("index", "z"))
vertices = None
def __init__(self, data_source, surface_field, field_value):
- ParallelAnalysisInterface.__init__(self)
self.data_source = data_source
self.surface_field = surface_field
self.field_value = field_value
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -41,6 +41,8 @@
from yt.fields.derived_field import \
ValidateSpatial
import yt.geometry.selection_routines
+from yt.geometry.selection_routines import \
+ compose_selector
from yt.extern.six import add_metaclass
def force_array(item, shape):
@@ -101,8 +103,15 @@
sets its initial set of fields, and the remainder of the arguments
are passed as field_parameters.
"""
- if ds != None:
+ # ds is typically set in the new object type created in Dataset._add_object_class
+ # but it can also be passed as a parameter to the constructor, in which case it will
+ # override the default. This code ensures it is never not set.
+ if ds is not None:
self.ds = ds
+ else:
+ if not hasattr(self, "ds"):
+ raise RuntimeError("Error: ds must be set either through class type or parameter to the constructor")
+
self._current_particle_type = "all"
self._current_fluid_type = self.ds.default_fluid_type
self.ds.objects.append(weakref.proxy(self))
@@ -542,10 +551,22 @@
_sort_by = None
_selector = None
_current_chunk = None
+ _data_source = None
+ _dimensionality = None
- def __init__(self, *args, **kwargs):
- super(YTSelectionContainer, self).__init__(*args, **kwargs)
-
+ def __init__(self, ds, field_parameters, data_source=None):
+ ParallelAnalysisInterface.__init__(self)
+ super(YTSelectionContainer, self).__init__(ds, field_parameters)
+ self._data_source = data_source
+ if data_source is not None:
+ if data_source.ds is not self.ds:
+ raise RuntimeError("Attempted to construct a DataContainer with a data_source from a different DataSet", ds, data_source.ds)
+ else:
+ print "DataSets: ", self.ds, data_source.ds
+ if data_source._dimensionality < self._dimensionality:
+ raise RuntimeError("Attempted to construct a DataContainer with a data_source of lower dimensionality (%u vs %u)" %
+ (data_source._dimensionality, self._dimensionality))
+
@property
def selector(self):
if self._selector is not None: return self._selector
@@ -555,7 +576,11 @@
"%s_selector" % self._type_name, None)
if sclass is None:
raise YTDataSelectorNotImplemented(self._type_name)
- self._selector = sclass(self)
+
+ if self._data_source is not None:
+ self._selector = compose_selector(self, self._data_source.selector, sclass(self))
+ else:
+ self._selector = sclass(self)
return self._selector
def chunks(self, fields, chunking_style, **kwargs):
@@ -765,30 +790,32 @@
class YTSelectionContainer0D(YTSelectionContainer):
_spatial = False
- def __init__(self, ds, field_parameters):
+ _dimensionality = 0
+ def __init__(self, ds, field_parameters = None, data_source = None):
super(YTSelectionContainer0D, self).__init__(
- ds, field_parameters)
+ ds, field_parameters, data_source)
class YTSelectionContainer1D(YTSelectionContainer):
_spatial = False
- def __init__(self, ds, field_parameters):
+ _dimensionality = 1
+ def __init__(self, ds, field_parameters = None, data_source = None):
super(YTSelectionContainer1D, self).__init__(
- ds, field_parameters)
+ ds, field_parameters, data_source)
self._grids = None
self._sortkey = None
self._sorted = {}
class YTSelectionContainer2D(YTSelectionContainer):
_key_fields = ['px','py','pdx','pdy']
+ _dimensionality = 2
"""
Prepares the YTSelectionContainer2D, normal to *axis*. If *axis* is 4, we are not
aligned with any axis.
"""
_spatial = False
- def __init__(self, axis, ds, field_parameters):
- ParallelAnalysisInterface.__init__(self)
+ def __init__(self, axis, ds, field_parameters = None, data_source = None):
super(YTSelectionContainer2D, self).__init__(
- ds, field_parameters)
+ ds, field_parameters, data_source)
# We need the ds, which will exist by now, for fix_axis.
self.axis = fix_axis(axis, self.ds)
self.set_field_parameter("axis", axis)
@@ -910,9 +937,9 @@
_key_fields = ['x','y','z','dx','dy','dz']
_spatial = False
_num_ghost_zones = 0
- def __init__(self, center, ds = None, field_parameters = None):
- ParallelAnalysisInterface.__init__(self)
- super(YTSelectionContainer3D, self).__init__(ds, field_parameters)
+ _dimensionality = 3
+ def __init__(self, center, ds, field_parameters = None, data_source = None):
+ super(YTSelectionContainer3D, self).__init__(ds, field_parameters, data_source)
self._set_center(center)
self.coords = None
self._grids = None
@@ -1273,9 +1300,9 @@
"""
_type_name = "boolean"
_con_args = ("regions",)
- def __init__(self, regions, fields = None, ds = None, **kwargs):
+ def __init__(self, regions, fields = None, ds = None, field_parameters = None, data_source = None):
# Center is meaningless, but we'll define it all the same.
- YTSelectionContainer3D.__init__(self, [0.5]*3, fields, ds, **kwargs)
+ YTSelectionContainer3D.__init__(self, [0.5]*3, fields, ds, field_parameters, data_source)
self.regions = regions
self._all_regions = []
self._some_overlap = []
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/data_objects/profiles.py
--- a/yt/data_objects/profiles.py
+++ b/yt/data_objects/profiles.py
@@ -761,6 +761,7 @@
self.field_data = YTFieldData()
if weight_field is not None:
self.variance = YTFieldData()
+ weight_field = self.data_source._determine_fields(weight_field)[0]
self.weight_field = weight_field
self.field_units = {}
ParallelAnalysisInterface.__init__(self, comm=data_source.comm)
@@ -774,7 +775,7 @@
A list of fields to create profile histograms for
"""
- fields = ensure_list(fields)
+ fields = self.data_source._determine_fields(fields)
temp_storage = ProfileFieldAccumulator(len(fields), self.size)
cfields = fields + list(self.bin_fields)
citer = self.data_source.chunks(cfields, "io")
@@ -907,9 +908,11 @@
if not np.any(filter): return None
arr = np.zeros((bin_fields[0].size, len(fields)), dtype="float64")
for i, field in enumerate(fields):
- arr[:,i] = chunk[field][filter]
+ units = chunk.ds.field_info[field].units
+ arr[:,i] = chunk[field][filter].in_units(units)
if self.weight_field is not None:
- weight_data = chunk[self.weight_field]
+ units = chunk.ds.field_info[self.weight_field].units
+ weight_data = chunk[self.weight_field].in_units(units)
else:
weight_data = np.ones(filter.size, dtype="float64")
weight_data = weight_data[filter]
@@ -1230,6 +1233,16 @@
self.z_bins.convert_to_units(new_unit)
self.z = 0.5*(self.z_bins[1:]+self.z_bins[:-1])
+
+def sanitize_field_tuple_keys(input_dict, data_source):
+ if input_dict is not None:
+ dummy = {}
+ for item in input_dict:
+ dummy[data_source._determine_fields(item)[0]] = input_dict[item]
+ return dummy
+ else:
+ return input_dict
+
def create_profile(data_source, bin_fields, fields, n_bins=64,
extrema=None, logs=None, units=None,
weight_field="cell_mass",
@@ -1242,7 +1255,7 @@
Parameters
----------
- data_source : AMR3DData Object
+ data_source : YTSelectionContainer Object
The data object to be profiled.
bin_fields : list of strings
List of the binning fields for profiling.
@@ -1293,7 +1306,7 @@
>>> print profile["gas", "temperature"]
"""
- bin_fields = ensure_list(bin_fields)
+ bin_fields = data_source._determine_fields(bin_fields)
fields = ensure_list(fields)
if len(bin_fields) == 1:
cls = Profile1D
@@ -1305,16 +1318,9 @@
raise NotImplementedError
bin_fields = data_source._determine_fields(bin_fields)
fields = data_source._determine_fields(fields)
- if units is not None:
- dummy = {}
- for item in units:
- dummy[data_source._determine_fields(item)[0]] = units[item]
- units.update(dummy)
- if extrema is not None:
- dummy = {}
- for item in extrema:
- dummy[data_source._determine_fields(item)[0]] = extrema[item]
- extrema.update(dummy)
+ units = sanitize_field_tuple_keys(units, data_source)
+ extrema = sanitize_field_tuple_keys(extrema, data_source)
+ logs = sanitize_field_tuple_keys(logs, data_source)
if weight_field is not None:
weight_field, = data_source._determine_fields([weight_field])
if not iterable(n_bins):
@@ -1322,18 +1328,21 @@
if not iterable(accumulation):
accumulation = [accumulation] * len(bin_fields)
if logs is None:
- logs = [data_source.ds._get_field_info(f[0],f[1]).take_log
- for f in bin_fields]
- else:
- logs = [logs[bin_field[-1]] for bin_field in bin_fields]
+ logs = {}
+ logs_list = []
+ for bin_field in bin_fields:
+ if bin_field in logs:
+ logs_list.append(logs[bin_field])
+ else:
+ logs_list.append(data_source.ds.field_info[bin_field].take_log)
+ logs = logs_list
if extrema is None:
ex = [data_source.quantities["Extrema"](f, non_zero=l)
for f, l in zip(bin_fields, logs)]
else:
ex = []
for bin_field in bin_fields:
- bf_units = data_source.ds._get_field_info(
- bin_field[0], bin_field[1]).units
+ bf_units = data_source.ds.field_info[bin_field].units
try:
field_ex = list(extrema[bin_field[-1]])
except KeyError:
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -51,8 +51,11 @@
ds: Dataset, optional
An optional dataset to use rather than self.ds
field_parameters : dictionary
- A dictionary of field parameters than can be accessed by derived
- fields.
+ A dictionary of field parameters than can be accessed by derived
+ fields.
+ data_source: optional
+ Draw the selection from the provided data source rather than
+ all data associated with the data_set
Examples
--------
@@ -64,8 +67,8 @@
"""
_type_name = "point"
_con_args = ('p',)
- def __init__(self, p, ds = None, field_parameters = None):
- super(YTPointBase, self).__init__(ds, field_parameters)
+ def __init__(self, p, ds=None, field_parameters=None, data_source=None):
+ super(YTPointBase, self).__init__(ds, field_parameters, data_source)
self.p = p
class YTOrthoRayBase(YTSelectionContainer1D):
@@ -92,6 +95,9 @@
field_parameters : dictionary
A dictionary of field parameters than can be accessed by derived
fields.
+ data_source: optional
+ Draw the selection from the provided data source rather than
+ all data associated with the data_set
Examples
--------
@@ -104,8 +110,9 @@
_key_fields = ['x','y','z','dx','dy','dz']
_type_name = "ortho_ray"
_con_args = ('axis', 'coords')
- def __init__(self, axis, coords, ds=None, field_parameters=None):
- super(YTOrthoRayBase, self).__init__(ds, field_parameters)
+ def __init__(self, axis, coords, ds=None,
+ field_parameters=None, data_source=None):
+ super(YTOrthoRayBase, self).__init__(ds, field_parameters, data_source)
self.axis = axis
xax = self.ds.coordinates.x_axis[self.axis]
yax = self.ds.coordinates.y_axis[self.axis]
@@ -144,6 +151,9 @@
field_parameters : dictionary
A dictionary of field parameters than can be accessed by derived
fields.
+ data_source: optional
+ Draw the selection from the provided data source rather than
+ all data associated with the data_set
Examples
--------
@@ -156,8 +166,9 @@
_type_name = "ray"
_con_args = ('start_point', 'end_point')
_container_fields = ("t", "dts")
- def __init__(self, start_point, end_point, ds=None, field_parameters=None):
- super(YTRayBase, self).__init__(ds, field_parameters)
+ def __init__(self, start_point, end_point, ds=None,
+ field_parameters=None, data_source=None):
+ super(YTRayBase, self).__init__(ds, field_parameters, data_source)
self.start_point = self.ds.arr(start_point,
'code_length', dtype='float64')
self.end_point = self.ds.arr(end_point,
@@ -204,6 +215,9 @@
field_parameters : dictionary
A dictionary of field parameters than can be accessed by derived
fields.
+ data_source: optional
+ Draw the selection from the provided data source rather than
+ all data associated with the data_set
Examples
--------
@@ -217,10 +231,10 @@
_type_name = "slice"
_con_args = ('axis', 'coord')
_container_fields = ("px", "py", "pdx", "pdy")
-
def __init__(self, axis, coord, center=None, ds=None,
- field_parameters = None):
- YTSelectionContainer2D.__init__(self, axis, ds, field_parameters)
+ field_parameters=None, data_source=None):
+ YTSelectionContainer2D.__init__(self, axis, ds,
+ field_parameters, data_source)
self._set_center(center)
self.coord = coord
@@ -285,6 +299,9 @@
field_parameters : dictionary
A dictionary of field parameters than can be accessed by derived
fields.
+ data_source: optional
+ Draw the selection from the provided data source rather than
+ all data associated with the data_set
Notes
-----
@@ -308,10 +325,10 @@
_type_name = "cutting"
_con_args = ('normal', 'center')
_container_fields = ("px", "py", "pz", "pdx", "pdy", "pdz")
-
- def __init__(self, normal, center, north_vector = None,
- ds = None, field_parameters = None):
- YTSelectionContainer2D.__init__(self, 4, ds, field_parameters)
+ def __init__(self, normal, center, north_vector=None,
+ ds=None, field_parameters=None, data_source=None):
+ YTSelectionContainer2D.__init__(self, 4, ds,
+ field_parameters, data_source)
self._set_center(center)
self.set_field_parameter('center',center)
# Let's set up our plane equation
@@ -465,7 +482,7 @@
Parameters
----------
- center : array_like
+ center : array_like
coordinate to which the normal, radius, and height all reference
normal : array_like
the normal vector defining the direction of lengthwise part of the
@@ -482,6 +499,9 @@
field_parameters : dictionary
A dictionary of field parameters than can be accessed by derived
fields.
+ data_source: optional
+ Draw the selection from the provided data source rather than
+ all data associated with the data_set
Examples
--------
@@ -494,8 +514,9 @@
_type_name = "disk"
_con_args = ('center', '_norm_vec', 'radius', 'height')
def __init__(self, center, normal, radius, height, fields=None,
- ds=None, **kwargs):
- YTSelectionContainer3D.__init__(self, center, fields, ds, **kwargs)
+ ds=None, field_parameters=None, data_source=None):
+ YTSelectionContainer3D.__init__(self, center, ds,
+ field_parameters, data_source)
self._norm_vec = np.array(normal)/np.sqrt(np.dot(normal,normal))
self.set_field_parameter("normal", self._norm_vec)
self.set_field_parameter("center", self.center)
@@ -523,9 +544,10 @@
"""
_type_name = "region"
_con_args = ('center', 'left_edge', 'right_edge')
- def __init__(self, center, left_edge, right_edge, fields = None,
- ds = None, **kwargs):
- YTSelectionContainer3D.__init__(self, center, ds, **kwargs)
+ def __init__(self, center, left_edge, right_edge, fields=None,
+ ds=None, field_parameters=None, data_source=None):
+ YTSelectionContainer3D.__init__(self, center, ds,
+ field_parameters, data_source)
if not isinstance(left_edge, YTArray):
self.left_edge = self.ds.arr(left_edge, 'code_length')
else:
@@ -542,8 +564,10 @@
"""
_type_name = "data_collection"
_con_args = ("_obj_list",)
- def __init__(self, obj_list, ds=None, field_parameters=None, center=None):
- YTSelectionContainer3D.__init__(self, center, ds, field_parameters)
+ def __init__(self, obj_list, ds=None, field_parameters=None,
+ data_source=None, center=None):
+ YTSelectionContainer3D.__init__(self, center, ds,
+ field_parameters, data_source)
self._obj_ids = np.array([o.id - o._id_offset for o in obj_list],
dtype="int64")
self._obj_list = obj_list
@@ -569,8 +593,10 @@
"""
_type_name = "sphere"
_con_args = ('center', 'radius')
- def __init__(self, center, radius, ds = None, field_parameters = None):
- super(YTSphereBase, self).__init__(center, ds, field_parameters)
+ def __init__(self, center, radius, ds=None,
+ field_parameters=None, data_source=None):
+ super(YTSphereBase, self).__init__(center, ds,
+ field_parameters, data_source)
# Unpack the radius, if necessary
radius = fix_length(radius, self.ds)
if radius < self.index.get_smallest_dx():
@@ -615,8 +641,9 @@
_type_name = "ellipsoid"
_con_args = ('center', '_A', '_B', '_C', '_e0', '_tilt')
def __init__(self, center, A, B, C, e0, tilt, fields=None,
- ds=None, field_parameters = None):
- YTSelectionContainer3D.__init__(self, center, ds, field_parameters)
+ ds=None, field_parameters=None, data_source=None):
+ YTSelectionContainer3D.__init__(self, center, ds,
+ field_parameters, data_source)
# make sure the magnitudes of semi-major axes are in order
if A<B or B<C:
raise YTEllipsoidOrdering(ds, A, B, C)
@@ -625,10 +652,10 @@
self._B = self.ds.quan(B, 'code_length')
self._C = self.ds.quan(C, 'code_length')
if self._C < self.index.get_smallest_dx():
- raise YTSphereTooSmall(ds, self._C, self.index.get_smallest_dx())
+ raise YTSphereTooSmall(self.ds, self._C, self.index.get_smallest_dx())
self._e0 = e0 = e0 / (e0**2.0).sum()**0.5
self._tilt = tilt
-
+
# find the t1 angle needed to rotate about z axis to align e0 to x
t1 = np.arctan(e0[1] / e0[0])
# rotate e0 by -t1
@@ -684,9 +711,10 @@
"""
_type_name = "cut_region"
_con_args = ("base_object", "conditionals")
- def __init__(self, base_object, conditionals, ds = None,
- field_parameters = None):
- super(YTCutRegionBase, self).__init__(base_object.center, ds, field_parameters)
+ def __init__(self, base_object, conditionals, ds=None,
+ field_parameters=None, data_source=None):
+ super(YTCutRegionBase, self).__init__(base_object.center, ds,
+ field_parameters, data_source)
self.conditionals = ensure_list(conditionals)
self.base_object = base_object
self._selector = None
@@ -762,4 +790,3 @@
@property
def fwidth(self):
return self.base_object.fwidth[self._cond_ind,:]
-
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/data_objects/tests/test_compose.py
--- /dev/null
+++ b/yt/data_objects/tests/test_compose.py
@@ -0,0 +1,146 @@
+from yt.testing import *
+from yt.fields.local_fields import add_field
+from yt.units.yt_array import YTArray, uintersect1d
+
+def setup():
+ from yt.config import ytcfg
+ ytcfg["yt","__withintesting"] = "True"
+
+# Copied from test_boolean for computing a unique identifier for
+# each cell from cell positions
+def _IDFIELD(field, data):
+ width = data.ds.domain_right_edge - data.ds.domain_left_edge
+ min_dx = YTArray(1.0/8192, input_units='code_length',
+ registry=data.ds.unit_registry)
+ delta = width / min_dx
+ x = data['x'] - min_dx / 2.
+ y = data['y'] - min_dx / 2.
+ z = data['z'] - min_dx / 2.
+ xi = x / min_dx
+ yi = y / min_dx
+ zi = z / min_dx
+ index = xi + delta[0] * (yi + delta[1] * zi)
+ index = index.astype('int64')
+ return index
+
+def test_compose_no_overlap():
+ r"""Test to make sure that composed data objects that don't
+ overlap behave the way we expect (return empty collections)
+ """
+ empty = np.array([])
+ for n in [1, 2, 4, 8]:
+ ds = fake_random_ds(64, nprocs=n)
+ ds.add_field("ID", function=_IDFIELD)
+
+ # position parameters for initial region
+ center = [0.25]*3
+ left_edge = [0.1]*3
+ right_edge = [0.4]*3
+ normal = [1, 0, 0]
+ radius = height = 0.15
+
+ # initial 3D regions
+ sources = [ds.sphere(center, radius),
+ ds.region(center, left_edge, right_edge),
+ ds.disk(center, normal, radius, height)]
+
+ # position parameters for non-overlapping regions
+ center = [0.75]*3
+ left_edge = [0.6]*3
+ right_edge = [0.9]*3
+
+ # subselect non-overlapping 0, 1, 2, 3D regions
+ for data1 in sources:
+ data2 = ds.sphere(center, radius, data_source=data1)
+ yield assert_array_equal, data2['ID'], empty
+
+ data2 = ds.region(center, left_edge, right_edge, data_source=data1)
+ yield assert_array_equal, data2['ID'], empty
+
+ data2 = ds.disk(center, normal, radius, height, data_source=data1)
+ yield assert_array_equal, data2['ID'], empty
+
+ for d in range(3):
+ data2 = ds.slice(d, center[d], data_source=data1)
+ yield assert_array_equal, data2['ID'], empty
+
+ for d in range(3):
+ data2 = ds.ortho_ray(d, center[0:d] + center[d+1:], data_source=data1)
+ yield assert_array_equal, data2['ID'], empty
+
+ data2 = ds.point(center, data_source=data1)
+ yield assert_array_equal, data2['ID'], empty
+
+def test_compose_overlap():
+ r"""Test to make sure that composed data objects that do
+ overlap behave the way we expect
+ """
+ empty = np.array([])
+ for n in [1, 2, 4, 8]:
+ ds = fake_random_ds(64, nprocs=n)
+ ds.add_field("ID", function=_IDFIELD)
+
+ # position parameters for initial region
+ center = [0.4, 0.5, 0.5]
+ left_edge = [0.1]*3
+ right_edge = [0.7]*3
+ normal = [1, 0, 0]
+ radius = height = 0.15
+
+ # initial 3D regions
+ sources = [ds.sphere(center, radius),
+ ds.region(center, left_edge, right_edge),
+ ds.disk(center, normal, radius, height)]
+
+ # position parameters for overlapping regions
+ center = [0.6, 0.5, 0.5]
+ left_edge = [0.3]*3
+ right_edge = [0.9]*3
+
+ # subselect non-overlapping 0, 1, 2, 3D regions
+ for data1 in sources:
+ id1 = data1['ID']
+
+ data2 = ds.sphere(center, radius)
+ data3 = ds.sphere(center, radius, data_source=data1)
+ id2 = data2['ID']
+ id3 = data3['ID']
+ id3.sort()
+ yield assert_array_equal, uintersect1d(id1, id2), id3
+
+ data2 = ds.region(center, left_edge, right_edge)
+ data3 = ds.region(center, left_edge, right_edge, data_source=data1)
+ id2 = data2['ID']
+ id3 = data3['ID']
+ id3.sort()
+ yield assert_array_equal, uintersect1d(id1, id2), id3
+
+ data2 = ds.disk(center, normal, radius, height)
+ data3 = ds.disk(center, normal, radius, height, data_source=data1)
+ id2 = data2['ID']
+ id3 = data3['ID']
+ id3.sort()
+ yield assert_array_equal, uintersect1d(id1, id2), id3
+
+ for d in range(3):
+ data2 = ds.slice(d, center[d])
+ data3 = ds.slice(d, center[d], data_source=data1)
+ id2 = data2['ID']
+ id3 = data3['ID']
+ id3.sort()
+ yield assert_array_equal, uintersect1d(id1, id2), id3
+
+ for d in range(3):
+ data2 = ds.ortho_ray(d, center[0:d] + center[d+1:])
+ data3 = ds.ortho_ray(d, center[0:d] + center[d+1:], data_source=data1)
+ id2 = data2['ID']
+ id3 = data3['ID']
+ id3.sort()
+ yield assert_array_equal, uintersect1d(id1, id2), id3
+
+ data2 = ds.point(center)
+ data3 = ds.point(center, data_source=data1)
+ id2 = data2['ID']
+ id3 = data3['ID']
+ id3.sort()
+ yield assert_array_equal, uintersect1d(id1, id2), id3
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/data_objects/tests/test_profiles.py
--- a/yt/data_objects/tests/test_profiles.py
+++ b/yt/data_objects/tests/test_profiles.py
@@ -1,7 +1,7 @@
from yt.testing import *
from yt.data_objects.profiles import \
BinnedProfile1D, BinnedProfile2D, BinnedProfile3D, \
- Profile1D, Profile2D, Profile3D
+ Profile1D, Profile2D, Profile3D, create_profile
_fields = ("density", "temperature", "dinosaurs", "tribbles")
_units = ("g/cm**3", "K", "dyne", "erg")
@@ -87,13 +87,26 @@
for nb in [8, 16, 32, 64]:
# We log all the fields or don't log 'em all. No need to do them
# individually.
- for lf in [True, False]:
- p1d = Profile1D(dd,
- "density", nb, rmi*e1, rma*e2, lf,
- weight_field = None)
- p1d.add_fields(["ones", "temperature"])
- yield assert_equal, p1d["ones"].sum(), nv
- yield assert_rel_equal, tt, p1d["temperature"].sum(), 7
+ for lf in [True, False]:
+ direct_profile = Profile1D(
+ dd, "density", nb, rmi*e1, rma*e2, lf, weight_field = None)
+ direct_profile.add_fields(["ones", "temperature"])
+
+ indirect_profile_s = create_profile(
+ dd, "density", ["ones", "temperature"], n_bins=nb,
+ extrema={'density': (rmi*e1, rma*e2)}, logs={'density': lf},
+ weight_field=None)
+
+ indirect_profile_t = create_profile(
+ dd, ("gas", "density"),
+ [("index", "ones"), ("gas", "temperature")], n_bins=nb,
+ extrema={'density': (rmi*e1, rma*e2)}, logs={'density': lf},
+ weight_field=None)
+
+ for p1d in [direct_profile, indirect_profile_s,
+ indirect_profile_t]:
+ yield assert_equal, p1d["index", "ones"].sum(), nv
+ yield assert_rel_equal, tt, p1d["gas", "temperature"].sum(), 7
p2d = Profile2D(dd,
"density", nb, rmi*e1, rma*e2, lf,
@@ -154,6 +167,12 @@
p3d.add_fields(["ones"])
yield assert_equal, p3d["ones"], np.ones((nb,nb,nb))
+extrema_s = {'particle_position_x': (0, 1)}
+logs_s = {'particle_position_x': False}
+
+extrema_t = {('all', 'particle_position_x'): (0, 1)}
+logs_t = {('all', 'particle_position_x'): False}
+
def test_particle_profiles():
for nproc in [1, 2, 4, 8]:
ds = fake_random_ds(32, nprocs=nproc, particles = 32**3)
@@ -164,6 +183,18 @@
p1d.add_fields(["particle_ones"])
yield assert_equal, p1d["particle_ones"].sum(), 32**3
+ p1d = create_profile(dd, ["particle_position_x"], ["particle_ones"],
+ weight_field=None, n_bins=128, extrema=extrema_s,
+ logs=logs_s)
+ yield assert_equal, p1d["particle_ones"].sum(), 32**3
+
+ p1d = create_profile(dd,
+ [("all", "particle_position_x")],
+ [("all", "particle_ones")],
+ weight_field=None, n_bins=128, extrema=extrema_t,
+ logs=logs_t)
+ yield assert_equal, p1d["particle_ones"].sum(), 32**3
+
p2d = Profile2D(dd, "particle_position_x", 128, 0.0, 1.0, False,
"particle_position_y", 128, 0.0, 1.0, False,
weight_field = None)
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/frontends/boxlib/fields.py
--- a/yt/frontends/boxlib/fields.py
+++ b/yt/frontends/boxlib/fields.py
@@ -97,6 +97,8 @@
# Now, let's figure out what fields are included.
if any(f[1] == "xmom" for f in self.field_list):
self.setup_momentum_to_velocity()
+ elif any(f[1] == "xvel" for f in self.field_list):
+ self.setup_velocity_to_momentum()
self.add_field(("gas", "thermal_energy"),
function=_thermal_energy,
units="erg/g")
@@ -112,11 +114,22 @@
def _get_vel(axis):
def velocity(field, data):
return data["%smom" % axis]/data["density"]
+ return velocity
for ax in 'xyz':
self.add_field(("gas", "velocity_%s" % ax),
function=_get_vel(ax),
units="cm/s")
+ def setup_velocity_to_momentum(self):
+ def _get_mom(axis):
+ def momentum(field, data):
+ return data["%svel" % axis]*data["density"]
+ return momentum
+ for ax in 'xyz':
+ self.add_field(("gas", "momentum_%s" % ax),
+ function=_get_mom(ax),
+ units=mom_units)
+
class CastroFieldInfo(FieldInfoContainer):
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/frontends/boxlib/tests/test_orion.py
--- a/yt/frontends/boxlib/tests/test_orion.py
+++ b/yt/frontends/boxlib/tests/test_orion.py
@@ -43,6 +43,14 @@
test_radtube.__name__ = test.description
yield test
+star = "StarParticles/plrd01000"
+ at requires_ds(star)
+def test_star():
+ ds = data_dir_load(star)
+ yield assert_equal, str(ds), "plrd01000"
+ for test in small_patch_amr(star, _fields):
+ test_star.__name__ = test.description
+ yield test
@requires_file(rt)
def test_OrionDataset():
diff -r e9d1135f9e91ad18f52fe9745789d68fa97d8225 -r 4c67e2a36f2c6f008a8360b7f7a3a4dd1dc04e28 yt/frontends/chombo/tests/test_outputs.py
--- a/yt/frontends/chombo/tests/test_outputs.py
+++ b/yt/frontends/chombo/tests/test_outputs.py
@@ -47,6 +47,15 @@
test_tb.__name__ = test.description
yield test
+iso = "IsothermalSphere/data.0000.3d.hdf5"
+ at requires_ds(iso)
+def test_iso():
+ ds = data_dir_load(iso)
+ yield assert_equal, str(ds), "data.0000.3d.hdf5"
+ for test in small_patch_amr(iso, _fields):
+ test_iso.__name__ = test.description
+ yield test
+
_zp_fields = ("rhs", "phi", "gravitational_field_x",
"gravitational_field_y")
zp = "ZeldovichPancake/plt32.2d.hdf5"
This diff is so big that we needed to truncate the remainder.
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list