[yt-svn] commit/yt: MatthewTurk: Merged in jzuhone/yt (pull request #1280)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Nov 6 07:33:43 PST 2014
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/e433a9791f5d/
Changeset: e433a9791f5d
Branch: yt
User: MatthewTurk
Date: 2014-11-06 15:33:31+00:00
Summary: Merged in jzuhone/yt (pull request #1280)
Improved FITS support, with a number of small additional refinements
Affected #: 20 files
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de 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 ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de 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 ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de 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 ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de 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 ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de 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 ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de 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 ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de 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,10 +115,19 @@
----------
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.data_containers.YTSelectionContainer, optional
@@ -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,10 +181,19 @@
----------
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.data_containers.YTSelectionContainer, optional
@@ -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 ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -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.
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/frontends/fits/data_structures.py
--- a/yt/frontends/fits/data_structures.py
+++ b/yt/frontends/fits/data_structures.py
@@ -42,11 +42,11 @@
unit_prefixes
from yt.units import dimensions
from yt.units.yt_array import YTQuantity
-from yt.utilities.on_demand_imports import _astropy
+from yt.utilities.on_demand_imports import _astropy, NotAModule
-lon_prefixes = ["X","RA","GLON"]
-lat_prefixes = ["Y","DEC","GLAT"]
-delimiters = ["*", "/", "-", "^"]
+lon_prefixes = ["X","RA","GLON","LINEAR"]
+lat_prefixes = ["Y","DEC","GLAT","LINEAR"]
+delimiters = ["*", "/", "-", "^", "(", ")"]
delimiters += [str(i) for i in xrange(10)]
regex_pattern = '|'.join(re.escape(_) for _ in delimiters)
@@ -552,7 +552,7 @@
x = 0
for p in lon_prefixes+lat_prefixes+spec_names.keys():
y = np_char.startswith(self.axis_names[:self.dimensionality], p)
- x += y.sum()
+ x += np.any(y)
if x == self.dimensionality: self._setup_spec_cube()
def _setup_spec_cube(self):
@@ -582,6 +582,12 @@
self.lon_axis = np.where(self.lon_axis)[0][0]
self.lon_name = ctypes[self.lon_axis].split("-")[0].lower()
+ if self.lat_axis == self.lon_axis and self.lat_name == self.lon_name:
+ self.lat_axis = 1
+ self.lon_axis = 0
+ self.lat_name = "Y"
+ self.lon_name = "X"
+
if self.wcs.naxis > 2:
self.spec_axis = np.zeros((end-1), dtype="bool")
@@ -658,6 +664,8 @@
ext = args[0].rsplit(".", 1)[0].rsplit(".", 1)[-1]
if ext.upper() not in ("FITS", "FTS"):
return False
+ elif isinstance(_astropy.pyfits, NotAModule):
+ raise RuntimeError("This appears to be a FITS file, but AstroPy is not installed.")
try:
with warnings.catch_warnings():
warnings.filterwarnings('ignore', category=UserWarning, append=True)
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/geometry/coordinates/coordinate_handler.py
--- a/yt/geometry/coordinates/coordinate_handler.py
+++ b/yt/geometry/coordinates/coordinate_handler.py
@@ -183,7 +183,15 @@
elif isinstance(center, YTArray):
return self.ds.arr(center), self.convert_to_cartesian(center)
elif iterable(center):
- if iterable(center[0]) and isinstance(center[1], basestring):
+ if isinstance(center[0], basestring) and isinstance(center[1], basestring):
+ if center[0].lower() == "min":
+ v, center = self.ds.find_min(center[1])
+ elif center[0].lower() == "max":
+ v, center = self.ds.find_max(center[1])
+ else:
+ raise RuntimeError("center keyword \"%s\" not recognized" % center)
+ center = self.ds.arr(center, 'code_length')
+ elif iterable(center[0]) and isinstance(center[1], basestring):
center = self.ds.arr(center[0], center[1])
else:
center = self.ds.arr(center, 'code_length')
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/geometry/coordinates/spec_cube_coordinates.py
--- a/yt/geometry/coordinates/spec_cube_coordinates.py
+++ b/yt/geometry/coordinates/spec_cube_coordinates.py
@@ -26,8 +26,19 @@
self.axis_name = {}
self.axis_id = {}
- for axis, axis_name in zip([ds.lon_axis, ds.lat_axis, ds.spec_axis],
- ["Image\ x", "Image\ y", ds.spec_name]):
+ self.default_unit_label = {}
+ if ds.lon_name == "X" and ds.lat_name == "Y":
+ names = ["x","y"]
+ else:
+ names = ["Image\ x", "Image\ y"]
+ self.default_unit_label[ds.lon_axis] = "pixel"
+ self.default_unit_label[ds.lat_axis] = "pixel"
+ names.append(ds.spec_name)
+ axes = [ds.lon_axis, ds.lat_axis, ds.spec_axis]
+ self.default_unit_label[ds.spec_axis] = ds.spec_unit
+
+ for axis, axis_name in zip(axes, names):
+
lower_ax = "xyz"[axis]
upper_ax = lower_ax.upper()
@@ -40,11 +51,6 @@
self.axis_id[axis] = axis
self.axis_id[axis_name] = axis
- self.default_unit_label = {}
- self.default_unit_label[ds.lon_axis] = "pixel"
- self.default_unit_label[ds.lat_axis] = "pixel"
- self.default_unit_label[ds.spec_axis] = ds.spec_unit
-
def _spec_axis(ax, x, y):
p = (x,y)[ax]
return [self.ds.pixel2spec(pp).v for pp in p]
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/units/unit_lookup_table.py
--- a/yt/units/unit_lookup_table.py
+++ b/yt/units/unit_lookup_table.py
@@ -172,6 +172,10 @@
'y': 1e-24, # yocto
}
+latex_prefixes = {
+ "u" : "\\mu",
+ }
+
prefixable_units = (
"m",
"pc",
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/units/unit_object.py
--- a/yt/units/unit_object.py
+++ b/yt/units/unit_object.py
@@ -26,7 +26,8 @@
from yt.units.unit_lookup_table import \
latex_symbol_lut, unit_prefixes, \
prefixable_units, cgs_base_units, \
- mks_base_units, cgs_conversions
+ mks_base_units, latex_prefixes, \
+ cgs_conversions
from yt.units.unit_registry import UnitRegistry
import copy
@@ -544,9 +545,14 @@
prefix_value = unit_prefixes[possible_prefix]
if symbol_str not in latex_symbol_lut:
+ if possible_prefix in latex_prefixes:
+ sstr = symbol_str.replace(possible_prefix,
+ '{'+latex_prefixes[possible_prefix]+'}')
+ else:
+ sstr = symbol_str
latex_symbol_lut[symbol_str] = \
latex_symbol_lut[symbol_wo_prefix].replace(
- '{'+symbol_wo_prefix+'}', '{'+symbol_str+'}')
+ '{'+symbol_wo_prefix+'}', '{'+sstr+'}')
# don't forget to account for the prefix value!
return (unit_data[0] * prefix_value, unit_data[1])
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -15,7 +15,8 @@
from yt.visualization.fixed_resolution import FixedResolutionBuffer
from yt.data_objects.construction_data_containers import YTCoveringGridBase
from yt.utilities.on_demand_imports import _astropy
-from yt.units.yt_array import YTQuantity
+from yt.units.yt_array import YTQuantity, YTArray
+import re
pyfits = _astropy.pyfits
pywcs = _astropy.pywcs
@@ -27,8 +28,7 @@
class FITSImageBuffer(HDUList):
- def __init__(self, data, fields=None, units="cm",
- center=None, scale=None, wcs=None):
+ def __init__(self, data, fields=None, units="cm", wcs=None):
r""" Initialize a FITSImageBuffer object.
FITSImageBuffer contains a list of FITS ImageHDU instances, and
@@ -49,29 +49,32 @@
keys, it will use these for the fields. If *data* is just a
single array one field name must be specified.
units : string
- The units of the WCS coordinates, default "cm".
- center : array_like, optional
- The coordinates [xctr,yctr] of the images in units
- *units*. If *units* is not specified, defaults to the origin.
- scale : tuple of floats, optional
- Pixel scale in unit *units*. Will be ignored if *data* is
- a FixedResolutionBuffer or a YTCoveringGrid. Must be
- specified otherwise, or if *units* is "deg".
+ The units of the WCS coordinates. Only applies
+ to FixedResolutionBuffers or YTCoveringGrids. Defaults to "cm".
wcs : `astropy.wcs.WCS` instance, optional
- Supply an AstroPy WCS instance to override automatic WCS creation.
+ Supply an AstroPy WCS instance. Will override automatic WCS
+ creation from FixedResolutionBuffers and YTCoveringGrids.
Examples
--------
+ >>> # This example uses a FRB.
>>> ds = load("sloshing_nomag2_hdf5_plt_cnt_0150")
>>> prj = ds.proj(2, "kT", weight_field="density")
>>> frb = prj.to_frb((0.5, "Mpc"), 800)
>>> # This example just uses the FRB and puts the coords in kpc.
>>> f_kpc = FITSImageBuffer(frb, fields="kT", units="kpc")
- >>> # This example specifies sky coordinates.
- >>> scale = [1./3600.]*2 # One arcsec per pixel
- >>> f_deg = FITSImageBuffer(frb, fields="kT", units="deg",
- scale=scale, center=(30., 45.))
+ >>> # This example specifies a specific WCS.
+ >>> from astropy.wcs import WCS
+ >>> w = WCS(naxis=self.dimensionality)
+ >>> w.wcs.crval = [30., 45.] # RA, Dec in degrees
+ >>> w.wcs.cunit = ["deg"]*2
+ >>> nx, ny = 800, 800
+ >>> w.wcs.crpix = [0.5*(nx+1), 0.5*(ny+1)]
+ >>> w.wcs.ctype = ["RA---TAN","DEC--TAN"]
+ >>> scale = 1./3600. # One arcsec per pixel
+ >>> w.wcs.cdelt = [-scale, scale]
+ >>> f_deg = FITSImageBuffer(frb, fields="kT", wcs=w)
>>> f_deg.writeto("temp.fits")
"""
@@ -98,8 +101,14 @@
first = True
+ self.field_units = {}
+
for key in fields:
if key not in exclude_fields:
+ if hasattr(img_data[key], "units"):
+ self.field_units[key] = str(img_data[key].units)
+ else:
+ self.field_units[key] = "dimensionless"
mylog.info("Making a FITS image of field %s" % (key))
if first:
hdu = pyfits.PrimaryHDU(np.array(img_data[key]))
@@ -109,7 +118,7 @@
hdu.name = key
hdu.header["btype"] = key
if hasattr(img_data[key], "units"):
- hdu.header["bunit"] = str(img_data[key].units)
+ hdu.header["bunit"] = re.sub('()', '', str(img_data[key].units))
self.append(hdu)
self.dimensionality = len(self[0].data.shape)
@@ -121,25 +130,11 @@
has_coords = (isinstance(img_data, FixedResolutionBuffer) or
isinstance(img_data, YTCoveringGridBase))
-
- if center is None:
- if units == "deg":
- mylog.error("Please specify center=(RA, Dec) in degrees.")
- raise ValueError
- elif not has_coords:
- mylog.warning("Setting center to the origin.")
- center = [0.0]*self.dimensionality
-
- if scale is None:
- if units == "deg" or not has_coords and wcs is None:
- mylog.error("Please specify scale=(dx,dy[,dz]) in %s." % (units))
- raise ValueError
if wcs is None:
w = pywcs.WCS(header=self[0].header, naxis=self.dimensionality)
w.wcs.crpix = 0.5*(np.array(self.shape)+1)
- proj_type = ["linear"]*self.dimensionality
- if isinstance(img_data, FixedResolutionBuffer) and units != "deg":
+ if isinstance(img_data, FixedResolutionBuffer):
# FRBs are a special case where we have coordinate
# information, so we take advantage of this and
# construct the WCS object
@@ -151,28 +146,20 @@
elif isinstance(img_data, YTCoveringGridBase):
dx, dy, dz = img_data.dds.in_units(units)
center = 0.5*(img_data.left_edge+img_data.right_edge).in_units(units)
- elif units == "deg" and self.dimensionality == 2:
- dx = -scale[0]
- dy = scale[1]
- proj_type = ["RA---TAN","DEC--TAN"]
else:
- dx = scale[0]
- dy = scale[1]
- if self.dimensionality == 3: dz = scale[2]
-
+ # We default to pixel coordinates if nothing is provided
+ dx, dy, dz = 1.0, 1.0, 1.0
+ center = 0.5*(np.array(self.shape)+1)
w.wcs.crval = center
- w.wcs.cunit = [units]*self.dimensionality
- w.wcs.ctype = proj_type
-
+ if has_coords:
+ w.wcs.cunit = [units]*self.dimensionality
if self.dimensionality == 2:
w.wcs.cdelt = [dx,dy]
elif self.dimensionality == 3:
w.wcs.cdelt = [dx,dy,dz]
-
+ w.wcs.ctype = ["linear"]*self.dimensionality
self._set_wcs(w)
-
else:
-
self._set_wcs(wcs)
def _set_wcs(self, wcs):
@@ -194,7 +181,7 @@
for img in self: img.header[key] = value
def keys(self):
- return [f.name for f in self]
+ return [f.name.lower() for f in self]
def has_key(self, key):
return key in self.keys()
@@ -249,33 +236,71 @@
import aplpy
return aplpy.FITSFigure(self, **kwargs)
+ def get_data(self, field):
+ return YTArray(self[field].data, self.field_units[field])
+
+ def set_unit(self, field, units):
+ """
+ Set the units of *field* to *units*.
+ """
+ new_data = YTArray(self[field].data, self.field_units[field]).in_units(units)
+ self[field].data = new_data.v
+ self[field].header["bunit"] = units
+ self.field_units[field] = units
+
axis_wcs = [[1,2],[0,2],[0,1]]
-def construct_image(data_source, center=None):
+def sanitize_fits_unit(unit):
+ if unit == "Mpc":
+ mylog.info("Changing FITS file unit to kpc.")
+ unit = "kpc"
+ elif unit == "au":
+ unit = "AU"
+ return unit
+
+def construct_image(data_source, center=None, width=None, image_res=None):
ds = data_source.ds
axis = data_source.axis
+ if center is None or width is None:
+ center = ds.domain_center[axis_wcs[axis]]
+ if width is None:
+ width = ds.domain_width[axis_wcs[axis]]
+ mylog.info("Making an image of the entire domain, "+
+ "so setting the center to the domain center.")
+ else:
+ width = ds.coordinates.sanitize_width(axis, width, None)
+ if image_res is None:
+ dd = ds.all_data()
+ dx, dy = [dd.quantities.extrema("d%s" % "xyz"[idx])[0]
+ for idx in axis_wcs[axis]]
+ nx = int((width[0]/dx).in_units("dimensionless"))
+ ny = int((width[1]/dy).in_units("dimensionless"))
+ else:
+ if iterable(image_res):
+ nx, ny = image_res
+ else:
+ nx, ny = image_res, image_res
+ dx, dy = width[0]/nx, width[1]/ny
+ crpix = [0.5*(nx+1), 0.5*(ny+1)]
if hasattr(ds, "wcs"):
- # This is a FITS dataset
- nx, ny = ds.domain_dimensions[axis_wcs[axis]]
- crpix = [ds.wcs.wcs.crpix[idx] for idx in axis_wcs[axis]]
- cdelt = [ds.wcs.wcs.cdelt[idx] for idx in axis_wcs[axis]]
- crval = [ds.wcs.wcs.crval[idx] for idx in axis_wcs[axis]]
+ # This is a FITS dataset, so we use it to construct the WCS
cunit = [str(ds.wcs.wcs.cunit[idx]) for idx in axis_wcs[axis]]
ctype = [ds.wcs.wcs.ctype[idx] for idx in axis_wcs[axis]]
+ cdelt = [ds.wcs.wcs.cdelt[idx] for idx in axis_wcs[axis]]
+ ctr_pix = center.in_units("code_length")[:self.dimensionality].v
+ crval = ds.wcs.wcs_pix2world(ctr_pix.reshape(1,self.dimensionality))[0]
+ crval = [crval[idx] for idx in axis_wcs[axis]]
else:
- # This is some other kind of dataset
- unit = ds.get_smallest_appropriate_unit(ds.domain_width.max())
- if center is None:
- crval = [0.0,0.0]
- else:
- crval = [(ds.domain_center-center)[idx].in_units(unit) for idx in axis_wcs[axis]]
- dx = ds.index.get_smallest_dx()
- nx, ny = (ds.domain_width[axis_wcs[axis]]/dx).ndarray_view().astype("int")
- crpix = [0.5*(nx+1), 0.5*(ny+1)]
- cdelt = [dx.in_units(unit)]*2
+ # This is some other kind of dataset
+ unit = str(width[0].units)
+ if unit == "unitary":
+ unit = ds.get_smallest_appropriate_unit(ds.domain_width.max())
+ unit = sanitize_fits_unit(unit)
cunit = [unit]*2
ctype = ["LINEAR"]*2
- frb = data_source.to_frb((1.0,"unitary"), (nx,ny))
+ cdelt = [dx.in_units(unit)]*2
+ crval = [center[idx].in_units(unit) for idx in axis_wcs[axis]]
+ frb = data_source.to_frb(width[0], (nx,ny), center=center, height=width[1])
w = pywcs.WCS(naxis=2)
w.wcs.crpix = crpix
w.wcs.cdelt = cdelt
@@ -296,21 +321,47 @@
The axis of the slice. One of "x","y","z", or 0,1,2.
fields : string or list of strings
The fields to slice
- center : A sequence floats, a string, or a tuple.
- The coordinate of the origin of the image. If set to 'c', 'center' or
+ 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. Units can be specified by passing in center
+ ('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
+ 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 windows with 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') 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
+ 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.
+ image_res : an int or 2-tuple of ints
+ Specify the resolution of the resulting image. If not provided, it will be
+ determined based on the minimum cell size of the dataset.
"""
- def __init__(self, ds, axis, fields, center="c", **kwargs):
+ def __init__(self, ds, axis, fields, center="c", width=None, **kwargs):
fields = ensure_list(fields)
axis = fix_axis(axis, ds)
- center = ds.coordinates.sanitize_center(center, axis)
+ center, dcenter = ds.coordinates.sanitize_center(center, axis)
slc = ds.slice(axis, center[axis], **kwargs)
- w, frb = construct_image(slc, center=center)
+ w, frb = construct_image(slc, center=dcenter, width=width,
+ image_res=image_res)
super(FITSSlice, self).__init__(frb, fields=fields, wcs=w)
for i, field in enumerate(fields):
self[i].header["bunit"] = str(frb[field].units)
@@ -329,21 +380,48 @@
The fields to project
weight_field : string
The field used to weight the projection.
- center : A sequence floats, a string, or a tuple.
- The coordinate of the origin 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. 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.
+ 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 windows with 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') 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
+ 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.
+ image_res : an int or 2-tuple of ints
+ Specify the resolution of the resulting image. If not provided, it will be
+ determined based on the minimum cell size of the dataset.
"""
- def __init__(self, ds, axis, fields, center="c", weight_field=None, **kwargs):
+ def __init__(self, ds, axis, fields, center="c", width=None,
+ weight_field=None, image_res=None, **kwargs):
fields = ensure_list(fields)
axis = fix_axis(axis, ds)
- center = ds.coordinates.sanitize_center(center, axis)
+ center, dcenter = ds.coordinates.sanitize_center(center, axis)
prj = ds.proj(fields[0], axis, weight_field=weight_field, **kwargs)
- w, frb = construct_image(prj, center=center)
+ w, frb = construct_image(prj, center=dcenter, width=width,
+ image_res=image_res)
super(FITSProjection, self).__init__(frb, fields=fields, wcs=w)
for i, field in enumerate(fields):
self[i].header["bunit"] = str(frb[field].units)
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/utilities/on_demand_imports.py
--- a/yt/utilities/on_demand_imports.py
+++ b/yt/utilities/on_demand_imports.py
@@ -132,4 +132,15 @@
self._interpolate = interpolate
return self._interpolate
+ _special = None
+ @property
+ def special(self):
+ if self._special is None:
+ try:
+ import scipy.special as special
+ except ImportError:
+ special = NotAModule(self._name)
+ self._special = special
+ return self._special
+
_scipy = scipy_imports()
\ No newline at end of file
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/visualization/fixed_resolution.py
--- a/yt/visualization/fixed_resolution.py
+++ b/yt/visualization/fixed_resolution.py
@@ -418,9 +418,9 @@
width, dd.resolution, item,
weight=dd.weight_field, volume=dd.volume,
no_ghost=dd.no_ghost, interpolated=dd.interpolated,
- north_vector=dd.north_vector)
+ north_vector=dd.north_vector, method=dd.method)
units = Unit(dd.ds.field_info[item].units, registry=dd.ds.unit_registry)
- if dd.weight_field is None:
+ if dd.weight_field is None and dd.method == "integrate":
units *= Unit('cm', registry=dd.ds.unit_registry)
ia = ImageArray(buff.swapaxes(0,1), input_units=units, info=self._get_info(item))
self[item] = ia
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -50,6 +50,8 @@
Unit
from yt.units.unit_registry import \
UnitParseError
+from yt.units.unit_lookup_table import \
+ prefixable_units, latex_prefixes
from yt.units.yt_array import \
YTArray, YTQuantity
from yt.utilities.png_writer import \
@@ -771,9 +773,10 @@
if hasattr(self.ds.coordinates, "default_unit_label"):
axax = getattr(self.ds.coordinates,
"%s_axis" % ("xy"[i]))[axis_index]
- un = self.ds.coordinates.default_unit_label[axax]
- axes_unit_labels[i] = '\/\/\left('+un+'\right)'
- continue
+ unn = self.ds.coordinates.default_unit_label.get(axax, "")
+ if unn != "":
+ axes_unit_labels[i] = '\/\/\left('+unn+'\right)'
+ continue
# Use sympy to factor h out of the unit. In this context 'un'
# is a string, so we call the Unit constructor.
expr = Unit(un, registry=self.ds.unit_registry).expr
@@ -794,6 +797,11 @@
raise RuntimeError
if un in formatted_length_unit_names:
un = formatted_length_unit_names[un]
+ pp = un[0]
+ if pp in latex_prefixes:
+ symbol_wo_prefix = un[1:]
+ if symbol_wo_prefix in prefixable_units:
+ un = un.replace(pp, "{"+latex_prefixes[pp]+"}", 1)
if un not in ['1', 'u', 'unitary']:
if hinv:
un = un + '\,h^{-1}'
@@ -940,13 +948,15 @@
or the axis name itself
fields : string
The name of the field(s) to be plotted.
- center : A sequence floats, a string, or a tuple.
+ 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. Units can be specified by passing in center
+ ('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
+ 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 windows with variable
@@ -1061,13 +1071,15 @@
or the axis name itself
fields : string
The name of the field(s) to be plotted.
- center : A sequence floats, a string, or a tuple.
+ 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. Units can be specified by passing in center
+ ('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
+ 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 windows with variable
@@ -1144,7 +1156,9 @@
"sum" : This method is the same as integrate, except that it does not
multiply by a path length when performing the integration, and is
- just a straight summation of the field along the given axis.
+ just a straight summation of the field along the given axis. WARNING:
+ This should only be used for uniform resolution grid datasets, as other
+ datasets may result in unphysical images.
proj_style : string
The method of projection--same as method keyword. Deprecated as of
version 3.0.2. Please use method instead.
@@ -1218,13 +1232,15 @@
The vector normal to the slicing plane.
fields : string
The name of the field(s) to be plotted.
- center : A sequence floats, a string, or a tuple.
+ 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. Units can be specified by passing in center
+ ('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
+ 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 windows with variable
@@ -1284,12 +1300,11 @@
class OffAxisProjectionDummyDataSource(object):
_type_name = 'proj'
- method = 'integrate'
_key_fields = []
def __init__(self, center, ds, normal_vector, width, fields,
interpolated, resolution = (800,800), weight=None,
volume=None, no_ghost=False, le=None, re=None,
- north_vector=None):
+ north_vector=None, method="integrate"):
self.center = center
self.ds = ds
self.axis = 4 # always true for oblique data objects
@@ -1306,6 +1321,7 @@
self.le = le
self.re = re
self.north_vector = north_vector
+ self.method = method
def _determine_fields(self, *args):
return self.dd._determine_fields(*args)
@@ -1329,13 +1345,15 @@
The vector normal to the slicing plane.
fields : string
The name of the field(s) to be plotted.
- center : A sequence floats, a string, or a tuple.
+ 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. Units can be specified by passing in center
+ ('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
+ 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 windows with variable
@@ -1376,7 +1394,20 @@
set, an arbitrary grid-aligned north-vector is chosen.
fontsize : integer
The size of the fonts for the axis, colorbar, and tick labels.
+ method : string
+ The method of projection. Valid methods are:
+ "integrate" with no weight_field specified : integrate the requested
+ field along the line of sight.
+
+ "integrate" with a weight_field specified : weight the requested
+ field by the weighting field and integrate along the line of sight.
+
+ "sum" : This method is the same as integrate, except that it does not
+ multiply by a path length when performing the integration, and is
+ just a straight summation of the field along the given axis. WARNING:
+ This should only be used for uniform resolution grid datasets, as other
+ datasets may result in unphysical images.
"""
_plot_type = 'OffAxisProjection'
_frb_generator = OffAxisProjectionFixedResolutionBuffer
@@ -1384,7 +1415,7 @@
def __init__(self, ds, normal, fields, center='c', width=None,
depth=(1, '1'), axes_unit=None, weight_field=None,
max_level=None, north_vector=None, volume=None, no_ghost=False,
- le=None, re=None, interpolated=False, fontsize=18):
+ le=None, re=None, interpolated=False, fontsize=18, method="integrate"):
(bounds, center_rot) = \
get_oblique_window_parameters(normal,center,width,ds,depth=depth)
fields = ensure_list(fields)[:]
@@ -1394,7 +1425,7 @@
OffAxisProj = OffAxisProjectionDummyDataSource(
center_rot, ds, normal, oap_width, fields, interpolated,
weight=weight_field, volume=volume, no_ghost=no_ghost,
- le=le, re=re, north_vector=north_vector)
+ le=le, re=re, north_vector=north_vector, method=method)
# If a non-weighted, integral projection, assure field-label reflects that
if weight_field is None and OffAxisProj.method == "integrate":
self.projected = True
@@ -1736,9 +1767,11 @@
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. Units can be specified by passing in center
+ ('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
+ 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 windows with variable
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py
+++ b/yt/visualization/volume_rendering/camera.py
@@ -2171,7 +2171,8 @@
class ProjectionCamera(Camera):
def __init__(self, center, normal_vector, width, resolution,
field, weight=None, volume=None, no_ghost = False,
- north_vector=None, ds=None, interpolated=False):
+ north_vector=None, ds=None, interpolated=False,
+ method="integrate"):
if not interpolated:
volume = 1
@@ -2180,6 +2181,7 @@
self.field = field
self.weight = weight
self.resolution = resolution
+ self.method = method
fields = [field]
if self.weight is not None:
@@ -2250,13 +2252,14 @@
dd = ds.all_data()
field = dd._determine_fields([self.field])[0]
finfo = ds._get_field_info(*field)
- if self.weight is None:
- dl = self.width[2]
- image *= dl
- else:
- image[:,:,0] /= image[:,:,1]
+ dl = 1.0
+ if self.method == "integrate":
+ if self.weight is None:
+ dl = self.width[2].in_units("cm")
+ else:
+ image[:,:,0] /= image[:,:,1]
- return image[:,:,0]
+ return ImageArray(image[:,:,0], finfo.units)*dl
def _render(self, double_check, num_threads, image, sampler):
@@ -2342,7 +2345,7 @@
def off_axis_projection(ds, center, normal_vector, width, resolution,
field, weight = None,
volume = None, no_ghost = False, interpolated = False,
- north_vector = None):
+ north_vector = None, method = "integrate"):
r"""Project through a dataset, off-axis, and return the image plane.
This function will accept the necessary items to integrate through a volume
@@ -2386,6 +2389,20 @@
If True, the data is first interpolated to vertex-centered data,
then tri-linearly interpolated along the ray. Not suggested for
quantitative studies.
+ method : string
+ The method of projection. Valid methods are:
+
+ "integrate" with no weight_field specified : integrate the requested
+ field along the line of sight.
+
+ "integrate" with a weight_field specified : weight the requested
+ field by the weighting field and integrate along the line of sight.
+
+ "sum" : This method is the same as integrate, except that it does not
+ multiply by a path length when performing the integration, and is
+ just a straight summation of the field along the given axis. WARNING:
+ This should only be used for uniform resolution grid datasets, as other
+ datasets may result in unphysical images.
Returns
-------
@@ -2403,7 +2420,7 @@
projcam = ProjectionCamera(center, normal_vector, width, resolution,
field, weight=weight, ds=ds, volume=volume,
no_ghost=no_ghost, interpolated=interpolated,
- north_vector=north_vector)
+ north_vector=north_vector, method=method)
image = projcam.snapshot()
return image[:,:]
diff -r ab0c5ed83cc0eff0c8f72e3faca3cecfe9cc93e6 -r e433a9791f5dadcd5de80c5d9f801526335455de yt/visualization/volume_rendering/image_handling.py
--- a/yt/visualization/volume_rendering/image_handling.py
+++ b/yt/visualization/volume_rendering/image_handling.py
@@ -40,9 +40,7 @@
data["b"] = image[:,:,2]
data["a"] = image[:,:,3]
nx, ny = data["r"].shape
- fib = FITSImageBuffer(data, units="pixel",
- center=[0.5*(nx+1), 0.5*(ny+1)],
- scale=[1.]*2)
+ fib = FITSImageBuffer(data)
fib.writeto('%s.fits'%fn,clobber=True)
def import_rgba(name, h5=True):
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