[yt-svn] commit/yt: 9 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Dec 3 09:43:30 PST 2014
9 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/47b9e7c00c67/
Changeset: 47b9e7c00c67
Branch: yt
User: jzuhone
Date: 2014-11-13 16:38:20+00:00
Summary: Setting up a function to create a celestial coordinate WCS from a cartesian one. Adding some keyword arguments back to FITSImageBuffer.
Affected #: 3 files
diff -r f7c75759e0395861b52d16921d8ce3ad6e36f89f -r 47b9e7c00c674a3d6aeddb4cf7d6f53c1434808a 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,7 +13,8 @@
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, sanitize_fits_unit
+from yt.utilities.fits_image import FITSImageBuffer, sanitize_fits_unit, \
+ create_sky_wcs
from yt.visualization.volume_rendering.camera import off_axis_projection
from yt.funcs import get_pbar
from yt.utilities.physical_constants import clight, mh
@@ -235,57 +236,40 @@
Whether or not to clobber an existing file with the same name.
length_unit : string
The units to convert the coordinates to in the file.
- sky_scale : tuple or YTQuantity
+ sky_scale : tuple, optional
Conversion between an angle unit and a length unit, if sky
- coordinates are desired.
- Examples: (1.0, "arcsec/kpc"), YTQuantity(0.001, "deg/kpc")
+ coordinates are desired, e.g. (1.0, "arcsec/kpc")
sky_center : tuple, optional
- The (RA, Dec) coordinate in degrees of the central pixel if
- *sky_scale* has been specified.
+ The (RA, Dec) coordinate in degrees of the central pixel. Must
+ be specified with *sky_scale*.
Examples
--------
>>> cube.write_fits("my_cube.fits", clobber=False, sky_scale=(1.0,"arcsec/kpc"))
"""
- 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"]
-
vunit = fits_info[self.axis_type][0]
vtype = fits_info[self.axis_type][1]
v_center = 0.5*(self.vbins[0]+self.vbins[-1]).in_units(vunit).value
- if sky_scale:
- dx = (self.width*sky_scale).in_units("deg").v/self.nx
- units = "deg"
+ if length_unit is None:
+ units = str(self.ds.get_smallest_appropriate_unit(self.width))
else:
- if length_unit is None:
- units = str(self.ds.get_smallest_appropriate_unit(self.width))
- else:
- units = length_unit
+ units = length_unit
units = sanitize_fits_unit(units)
dx = self.width.in_units(units).v/self.nx
- dy = dx
+ dy = self.width.in_units(units).v/self.ny
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.crval = [0.0,0.0,v_center]
w.wcs.cunit = [units,units,vunit]
- w.wcs.ctype = [types[0],types[1],vtype]
+ w.wcs.ctype = ["LINEAR","LINEAR",vtype]
+
+ if sky_scale is not None and sky_center is not None:
+ w = create_sky_wcs(w, sky_center, sky_scale)
fib = FITSImageBuffer(self.data.transpose(), fields=self.field, wcs=w)
fib[0].header["bunit"] = re.sub('()', '', str(self.proj_units))
diff -r f7c75759e0395861b52d16921d8ce3ad6e36f89f -r 47b9e7c00c674a3d6aeddb4cf7d6f53c1434808a yt/analysis_modules/sunyaev_zeldovich/projection.py
--- a/yt/analysis_modules/sunyaev_zeldovich/projection.py
+++ b/yt/analysis_modules/sunyaev_zeldovich/projection.py
@@ -298,13 +298,12 @@
----------
filename : string
The name of the FITS file to be written.
- sky_scale : tuple or YTQuantity
+ sky_scale : tuple
Conversion between an angle unit and a length unit, if sky
- coordinates are desired.
- Examples: (1.0, "arcsec/kpc"), YTQuantity(0.001, "deg/kpc")
+ coordinates are desired, e.g. (1.0, "arcsec/kpc")
sky_center : tuple, optional
- The (RA, Dec) coordinate in degrees of the central pixel if
- *sky_scale* has been specified.
+ The (RA, Dec) coordinate in degrees of the central pixel. Must
+ be specified with *sky_scale*.
clobber : boolean, optional
If the file already exists, do we overwrite?
@@ -314,40 +313,25 @@
>>> szprj.write_fits("SZbullet.fits", clobber=False)
>>> # This example uses sky coords
>>> sky_scale = (1., "arcsec/kpc") # One arcsec per kpc
- >>> sky_center = (30., 45.) # In degrees
+ >>> sky_center = (30., 45., "deg")
>>> szprj.write_fits("SZbullet.fits", sky_center=sky_center, sky_scale=sky_scale)
"""
- from yt.utilities.fits_image import FITSImageBuffer, sanitize_fits_unit
-
- 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"]
+ from yt.utilities.fits_image import FITSImageBuffer, sanitize_fits_unit, create_sky_wcs
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"
dy = dx
- if sky_scale:
- dx *= -1.
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.crval = [0.0,0.0]
w.wcs.cunit = [units]*2
- w.wcs.ctype = types
+ w.wcs.ctype = ["LINEAR"]*2
+
+ if sky_scale is not None and sky_center is not None:
+ w = create_sky_wcs(w, sky_center, sky_scale)
fib = FITSImageBuffer(self.data, fields=self.data.keys(), wcs=w)
fib.writeto(filename, clobber=clobber)
diff -r f7c75759e0395861b52d16921d8ce3ad6e36f89f -r 47b9e7c00c674a3d6aeddb4cf7d6f53c1434808a yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -28,7 +28,7 @@
class FITSImageBuffer(HDUList):
- def __init__(self, data, fields=None, units="cm", wcs=None):
+ def __init__(self, data, fields=None, units=None, pixel_scale=None, wcs=None):
r""" Initialize a FITSImageBuffer object.
FITSImageBuffer contains a list of FITS ImageHDU instances, and
@@ -49,8 +49,11 @@
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. Only applies
- to FixedResolutionBuffers or YTCoveringGrids. Defaults to "cm".
+ The units of the WCS coordinates. Defaults to "cm".
+ pixel_scale : float
+ The scale of the pixel, in *units*. Either a single float or
+ iterable of floats. Only used if this information is not already
+ provided by *data*.
wcs : `astropy.wcs.WCS` instance, optional
Supply an AstroPy WCS instance. Will override automatic WCS
creation from FixedResolutionBuffers and YTCoveringGrids.
@@ -77,7 +80,10 @@
>>> f_deg = FITSImageBuffer(frb, fields="kT", wcs=w)
>>> f_deg.writeto("temp.fits")
"""
-
+
+ if units is None: units = "cm"
+ if pixel_scale is None: pixel_scale = 1.0
+
super(FITSImageBuffer, self).__init__()
if isinstance(fields, basestring): fields = [fields]
@@ -91,13 +97,13 @@
img_data = {}
if fields is None:
mylog.error("Please specify a field name for this array.")
- raise KeyError
+ raise KeyError("Please specify a field name for this array.")
img_data[fields[0]] = data
if fields is None: fields = img_data.keys()
if len(fields) == 0:
mylog.error("Please specify one or more fields to write.")
- raise KeyError
+ raise KeyError("Please specify one or more fields to write.")
first = True
@@ -128,12 +134,8 @@
elif self.dimensionality == 3:
self.nx, self.ny, self.nz = self[0].data.shape
- has_coords = (isinstance(img_data, FixedResolutionBuffer) or
- isinstance(img_data, YTCoveringGridBase))
-
if wcs is None:
w = pywcs.WCS(header=self[0].header, naxis=self.dimensionality)
- w.wcs.crpix = 0.5*(np.array(self.shape)+1)
if isinstance(img_data, FixedResolutionBuffer):
# FRBs are a special case where we have coordinate
# information, so we take advantage of this and
@@ -143,26 +145,28 @@
xctr = 0.5*(img_data.bounds[1]+img_data.bounds[0]).in_units(units)
yctr = 0.5*(img_data.bounds[3]+img_data.bounds[2]).in_units(units)
center = [xctr, yctr]
+ cdelt = [dx,dy]
elif isinstance(img_data, YTCoveringGridBase):
- dx, dy, dz = img_data.dds.in_units(units)
+ cdelt = img_data.dds.in_units(units).v
center = 0.5*(img_data.left_edge+img_data.right_edge).in_units(units)
else:
- # 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)
+ # If img_data is just an array, we assume the center is the origin
+ # and use *pixel_scale* to determine the cell widths
+ if iterable(pixel_scale):
+ cdelt = pixel_scale
+ else:
+ cdelt = [pixel_scale]*self.dimensionality
+ center = [0.0]*self.dimensionality
+ w.wcs.crpix = 0.5*(np.array(self.shape)+1)
w.wcs.crval = center
- 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.cdelt = cdelt
w.wcs.ctype = ["linear"]*self.dimensionality
- self._set_wcs(w)
+ w.wcs.cunit = [units]*self.dimensionality
+ self.set_wcs(w)
else:
- self._set_wcs(wcs)
+ self.set_wcs(wcs)
- def _set_wcs(self, wcs):
+ def set_wcs(self, wcs):
"""
Set the WCS coordinate information for all images
with a WCS object *wcs*.
@@ -250,6 +254,50 @@
axis_wcs = [[1,2],[0,2],[0,1]]
+def create_sky_wcs(old_wcs, sky_center, sky_scale,
+ ctype=["RA---TAN","DEC--TAN"], crota=None):
+ """
+ Takes an astropy.wcs.WCS instance created in yt from a
+ simulation that has a Cartesian coordinate system and
+ converts it to one in a celestial coordinate system.
+
+ Parameters
+ ----------
+ old_wcs : astropy.wcs.WCS
+ The original WCS to be converted.
+ sky_center : tuple
+ Reference coordinates of the WCS in degrees.
+ sky_scale : tuple
+ Conversion between an angle unit and a length unit,
+ e.g. (3.0, "arcsec/kpc")
+ ctype : list of strings, optional
+ The type of the coordinate system to create.
+ crota : list of floats, optional
+ Rotation angles between cartesian coordinates and
+ the celestial coordinates.
+ """
+ naxis = old_wcs.naxis
+ crval = [sky_center[0], sky_center[1]]
+ scaleq = YTQuantity(sky_scale[0],sky_scale[1])
+ dx, dy = old_wcs.wcs.cdelt
+ xunit, yunit = [str(unit) for unit in old_wcs.wcs.cunit]
+ new_dx = (YTQuantity(-dx, xunit)*scaleq).in_units("deg")
+ new_dy = (YTQuantity(dy, yunit)*scaleq).in_units("deg")
+ new_wcs = pywcs.WCS(naxis=naxis)
+ cdelt = [new_dx.v, new_dy.v]
+ if naxis == 3:
+ crval.append(old_wcs.wcs.crval[2])
+ cdelt.append(old_wcs.wcs.cdelt[2])
+ ctype.append(old_wcs.wcs.ctype[2])
+ new_wcs.wcs.crpix = old_wcs.wcs.crpix
+ new_wcs.wcs.cdelt = cdelt
+ new_wcs.wcs.crval = crval
+ new_wcs.wcs.cunit = ["deg"]*naxis
+ new_wcs.wcs.ctype = ctype
+ if crota is not None:
+ new_wcs.wcs.crota = crota
+ return new_wcs
+
def sanitize_fits_unit(unit):
if unit == "Mpc":
mylog.info("Changing FITS file unit to kpc.")
https://bitbucket.org/yt_analysis/yt/commits/2f28f708378c/
Changeset: 2f28f708378c
Branch: yt
User: jzuhone
Date: 2014-12-02 20:53:17+00:00
Summary: Making sure create_sky_wcs works for 3D FITS cubes. Also fixing up the PPVCube docs.
Affected #: 2 files
diff -r 47b9e7c00c674a3d6aeddb4cf7d6f53c1434808a -r 2f28f708378c953c424db2b7ebf4a8d7667a195f doc/source/analyzing/analysis_modules/PPVCube.ipynb
--- a/doc/source/analyzing/analysis_modules/PPVCube.ipynb
+++ b/doc/source/analyzing/analysis_modules/PPVCube.ipynb
@@ -1,421 +1,451 @@
{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Detailed spectra of astrophysical objects sometimes allow for determinations of how much of the gas is moving with a certain velocity along the line of sight, thanks to Doppler shifting of spectral lines. This enables \"data cubes\" to be created in RA, Dec, and line-of-sight velocity space. In yt, we can use the `PPVCube` analysis module to project fields along a given line of sight traveling at different line-of-sight velocities, to \"mock-up\" what would be seen in observations."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "import yt\n",
+ "import numpy as np\n",
+ "\n",
+ "from yt.analysis_modules.ppv_cube.api import PPVCube"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "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:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Density: $\\rho(r) \\propto r^{\\alpha}$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Velocity: $v_{\\theta}(r) \\propto \\frac{r}{1+(r/r_0)^{\\beta}}$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where for simplicity we won't worry about the normalizations of these profiles. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First, we'll set up the grid and the parameters of the profiles:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "nx,ny,nz = (256,256,256) # domain dimensions\n",
+ "R = 10. # outer radius of disk, kpc\n",
+ "r_0 = 3. # scale radius, kpc\n",
+ "beta = 1.4 # for the tangential velocity profile\n",
+ "alpha = -1. # for the radial density profile\n",
+ "x, y = np.mgrid[-R:R:nx*1j,-R:R:ny*1j] # cartesian coordinates of x-y plane of disk\n",
+ "r = np.sqrt(x*x+y*y) # polar coordinates\n",
+ "theta = np.arctan2(y, x) # polar coordinates"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "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. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "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",
+ "temp = np.zeros((nx,ny,nz))\n",
+ "temp[:,:,nz/2-3:nz/2+3] = 1.0e4 # 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"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Finally, we'll package these data arrays up into a dictionary, which will then be shipped off to `load_uniform_grid`. We'll define the width of the grid to be `2*R` kpc, which will be equal to 1 `code_length`. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "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",
+ "bbox = np.array([[-0.5,0.5],[-0.5,0.5],[-0.5,0.5]]) # bbox of width 1 on a side with center (0,0,0)\n",
+ "ds = yt.load_uniform_grid(data, (nx,ny,nz), length_unit=(2*R,\"kpc\"), nprocs=1, bbox=bbox)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To get a sense of what the data looks like, we'll take a slice through the middle of the disk:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "slc = yt.SlicePlot(ds, \"z\", [\"density\",\"velocity_x\",\"velocity_y\",\"velocity_magnitude\"])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "slc.set_log(\"velocity_x\", False)\n",
+ "slc.set_log(\"velocity_y\", False)\n",
+ "slc.set_log(\"velocity_magnitude\", False)\n",
+ "slc.set_unit(\"velocity_magnitude\", \"km/s\")\n",
+ "slc.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Which shows a rotating disk with a specific density and velocity profile. Now, suppose we wanted to look at this disk galaxy from a certain orientation angle, and simulate a 3D FITS data cube where we can see the gas that is emitting at different velocities along the line of sight. We can do this using the `PPVCube` class. First, let's assume we rotate our viewing angle 60 degrees from face-on, from along the z-axis into the y-axis. We'll create a normal vector:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "i = 60.*np.pi/180.\n",
+ "L = [0.0,np.sin(i),np.sin(i)]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next, we need to specify a field that will serve as the \"intensity\" of the emission that we see. For simplicity, we'll simply choose the gas density as this field, though it could be any field (including derived fields) in principle. We also need to specify the dimensions of the data cube, and optionally we may choose the bounds in line-of-sight velocity that the data will be binned into. Otherwise, the bounds will simply be set to the negative and positive of the largest speed in the dataset."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "cube = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-150.,150.,\"km/s\"), method=\"sum\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "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",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "cube.write_fits(\"cube.fits\", clobber=True, length_unit=\"kpc\")"
+ ]
+ },
+ {
+ "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",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "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)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now, we'll look at the FITS dataset in yt and look at different slices along the velocity axis, which is the \"z\" axis:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "ds_cube = yt.load(\"cube.fits\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Specifying no center gives us the center slice\n",
+ "slc = yt.SlicePlot(ds_cube, \"z\", [\"density\"])\n",
+ "slc.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "import yt.units as u\n",
+ "# Picking different velocities for the slices\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()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "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()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "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()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "If we project all the emission at all the different velocities along the z-axis, we recover the entire disk:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "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()"
+ ]
+ },
+ {
+ "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",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "cube2 = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-150,150,\"km/s\"), thermal_broad=True, \n",
+ " atomic_weight=12.0, method=\"sum\")\n",
+ "cube2.write_fits(\"cube2.fits\", clobber=True, length_unit=\"kpc\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Taking a slice of this cube shows:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "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()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "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()"
+ ]
+ },
+ {
+ "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",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "print cube2.vbins[0], cube2.vbins[-1]\n",
+ "cube2.transform_spectral_axis(400.0,\"nm\")\n",
+ "print cube2.vbins[0], cube2.vbins[-1]"
+ ]
+ },
+ {
+ "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",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "cube2.reset_spectral_axis()\n",
+ "print cube2.vbins[0], cube2.vbins[-1]"
+ ]
+ }
+ ],
"metadata": {
- "name": "",
- "signature": "sha256:b83e125278c2e58da4d99ac9d2ca2a136d01f1094e1b83497925e0f9b9b056c2"
+ "kernelspec": {
+ "display_name": "IPython (Python 2)",
+ "name": "python2"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 2
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython2"
+ },
+ "signature": "sha256:81516d7b1f41271d88a9c89eb3fed25537418d40db93a50f34f3f58ac7679dc5"
},
- "nbformat": 3,
- "nbformat_minor": 0,
- "worksheets": [
- {
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Detailed spectra of astrophysical objects sometimes allow for determinations of how much of the gas is moving with a certain velocity along the line of sight, thanks to Doppler shifting of spectral lines. This enables \"data cubes\" to be created in RA, Dec, and line-of-sight velocity space. In yt, we can use the `PPVCube` analysis module to project fields along a given line of sight traveling at different line-of-sight velocities, to \"mock-up\" what would be seen in observations."
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "import yt\n",
- "import numpy as np\n",
- "\n",
- "from yt.analysis_modules.ppv_cube.api import PPVCube"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "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:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Density: $\\rho(r) \\propto r^{\\alpha}$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Velocity: $v_{\\theta}(r) \\propto \\frac{r}{1+(r/r_0)^{\\beta}}$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "where for simplicity we won't worry about the normalizations of these profiles. "
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "First, we'll set up the grid and the parameters of the profiles:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "nx,ny,nz = (256,256,256) # domain dimensions\n",
- "R = 10. # outer radius of disk, kpc\n",
- "r_0 = 3. # scale radius, kpc\n",
- "beta = 1.4 # for the tangential velocity profile\n",
- "alpha = -1. # for the radial density profile\n",
- "x, y = np.mgrid[-R:R:nx*1j,-R:R:ny*1j] # cartesian coordinates of x-y plane of disk\n",
- "r = np.sqrt(x*x+y*y) # polar coordinates\n",
- "theta = np.arctan2(y, x) # polar coordinates"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "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. "
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "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",
- "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"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Finally, we'll package these data arrays up into a dictionary, which will then be shipped off to `load_uniform_grid`. We'll define the width of the grid to be `2*R` kpc, which will be equal to 1 `code_length`. "
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "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",
- "bbox = np.array([[-0.5,0.5],[-0.5,0.5],[-0.5,0.5]]) # bbox of width 1 on a side with center (0,0,0)\n",
- "ds = yt.load_uniform_grid(data, (nx,ny,nz), length_unit=(2*R,\"kpc\"), nprocs=1, bbox=bbox)"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "To get a sense of what the data looks like, we'll take a slice through the middle of the disk:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "slc = yt.SlicePlot(ds, \"z\", [\"density\",\"velocity_x\",\"velocity_y\",\"velocity_magnitude\"])"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "slc.set_log(\"velocity_x\", False)\n",
- "slc.set_log(\"velocity_y\", False)\n",
- "slc.set_log(\"velocity_magnitude\", False)\n",
- "slc.set_unit(\"velocity_magnitude\", \"km/s\")\n",
- "slc.show()"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Which shows a rotating disk with a specific density and velocity profile. Now, suppose we wanted to look at this disk galaxy from a certain orientation angle, and simulate a 3D FITS data cube where we can see the gas that is emitting at different velocities along the line of sight. We can do this using the `PPVCube` class. First, let's assume we rotate our viewing angle 60 degrees from face-on, from along the z-axis into the y-axis. We'll create a normal vector:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "i = 60.*np.pi/180.\n",
- "L = [0.0,np.sin(i),np.sin(i)]"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Next, we need to specify a field that will serve as the \"intensity\" of the emission that we see. For simplicity, we'll simply choose the gas density as this field, though it could be any field (including derived fields) in principle. We also need to specify the dimensions of the data cube, and optionally we may choose the bounds in line-of-sight velocity that the data will be binned into. Otherwise, the bounds will simply be set to the negative and positive of the largest speed in the dataset."
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "cube = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-150.,150.,\"km/s\"))"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "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=\"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": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Now, we'll look at the FITS dataset in yt and look at different slices along the velocity axis, which is the \"z\" axis:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "ds_cube = yt.load(\"cube.fits\")"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "# Specifying no center gives us the center slice\n",
- "slc = yt.SlicePlot(ds_cube, \"z\", [\"density\"])\n",
- "slc.show()"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "import yt.units as u\n",
- "# Picking different velocities for the slices\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",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "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",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "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",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "If we project all the emission at all the different velocities along the z-axis, we recover the entire disk:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "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()"
- ],
- "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": {}
- }
- ]
+ "nbformat": 4,
+ "nbformat_minor": 0
}
\ No newline at end of file
diff -r 47b9e7c00c674a3d6aeddb4cf7d6f53c1434808a -r 2f28f708378c953c424db2b7ebf4a8d7667a195f yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -279,10 +279,10 @@
naxis = old_wcs.naxis
crval = [sky_center[0], sky_center[1]]
scaleq = YTQuantity(sky_scale[0],sky_scale[1])
- dx, dy = old_wcs.wcs.cdelt
- xunit, yunit = [str(unit) for unit in old_wcs.wcs.cunit]
- new_dx = (YTQuantity(-dx, xunit)*scaleq).in_units("deg")
- new_dy = (YTQuantity(dy, yunit)*scaleq).in_units("deg")
+ deltas = old_wcs.wcs.cdelt
+ units = [str(unit) for unit in old_wcs.wcs.cunit]
+ new_dx = (YTQuantity(-deltas[0], units[0])*scaleq).in_units("deg")
+ new_dy = (YTQuantity(deltas[1], units[1])*scaleq).in_units("deg")
new_wcs = pywcs.WCS(naxis=naxis)
cdelt = [new_dx.v, new_dy.v]
if naxis == 3:
https://bitbucket.org/yt_analysis/yt/commits/0e8cedc03e2f/
Changeset: 0e8cedc03e2f
Branch: yt
User: jzuhone
Date: 2014-12-02 21:18:51+00:00
Summary: For simplicity, change this back to kpc
Affected #: 1 file
diff -r 2f28f708378c953c424db2b7ebf4a8d7667a195f -r 0e8cedc03e2f67c6d5f27b9612687bed952acbf7 yt/analysis_modules/sunyaev_zeldovich/projection.py
--- a/yt/analysis_modules/sunyaev_zeldovich/projection.py
+++ b/yt/analysis_modules/sunyaev_zeldovich/projection.py
@@ -316,18 +316,16 @@
>>> sky_center = (30., 45., "deg")
>>> szprj.write_fits("SZbullet.fits", sky_center=sky_center, sky_scale=sky_scale)
"""
- from yt.utilities.fits_image import FITSImageBuffer, sanitize_fits_unit, create_sky_wcs
+ from yt.utilities.fits_image import FITSImageBuffer, create_sky_wcs
- units = self.ds.get_smallest_appropriate_unit(self.width)
- units = sanitize_fits_unit(units)
- dx = self.dx.in_units(units)
+ dx = self.dx.in_units("kpc")
dy = dx
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 = [0.0,0.0]
- w.wcs.cunit = [units]*2
+ w.wcs.cunit = ["kpc"]*2
w.wcs.ctype = ["LINEAR"]*2
if sky_scale is not None and sky_center is not None:
https://bitbucket.org/yt_analysis/yt/commits/5c0367b70806/
Changeset: 5c0367b70806
Branch: yt
User: jzuhone
Date: 2014-12-02 21:21:12+00:00
Summary: Change temp back to where it should go
Affected #: 1 file
diff -r 0e8cedc03e2f67c6d5f27b9612687bed952acbf7 -r 5c0367b708068edc612f2e5694c6b8e54d3f2c37 doc/source/analyzing/analysis_modules/PPVCube.ipynb
--- a/doc/source/analyzing/analysis_modules/PPVCube.ipynb
+++ b/doc/source/analyzing/analysis_modules/PPVCube.ipynb
@@ -92,7 +92,7 @@
"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",
"temp = np.zeros((nx,ny,nz))\n",
- "temp[:,:,nz/2-3:nz/2+3] = 1.0e4 # Isothermal\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",
@@ -444,7 +444,7 @@
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2"
},
- "signature": "sha256:81516d7b1f41271d88a9c89eb3fed25537418d40db93a50f34f3f58ac7679dc5"
+ "signature": "sha256:c0d10b9e420fd80ed69d1ad4dcfb9eb852999c7bf8524868f2df19e940e1fae8"
},
"nbformat": 4,
"nbformat_minor": 0
https://bitbucket.org/yt_analysis/yt/commits/506e06f0427a/
Changeset: 506e06f0427a
Branch: yt
User: jzuhone
Date: 2014-12-02 21:22:05+00:00
Summary: Get rid of these fields so they don't screw us up later
Affected #: 1 file
diff -r 5c0367b708068edc612f2e5694c6b8e54d3f2c37 -r 506e06f0427a6cacdc4679813fbf10cade833927 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
@@ -179,6 +179,9 @@
else:
self.width = ds.quan(self.width, "code_length")
+ self.ds.field_info.pop(("gas","intensity"))
+ self.ds.field_info.pop(("gas","v_los"))
+
def create_intensity(self):
def _intensity(field, data):
v = self.current_v-data["v_los"].v
https://bitbucket.org/yt_analysis/yt/commits/afb117e019cc/
Changeset: afb117e019cc
Branch: yt
User: jzuhone
Date: 2014-12-02 21:36:33+00:00
Summary: Fixing my mistake with the wrong notebook format
Affected #: 1 file
diff -r 506e06f0427a6cacdc4679813fbf10cade833927 -r afb117e019cc3e1a6b525c7d216f8a5ce312055a doc/source/analyzing/analysis_modules/PPVCube.ipynb
--- a/doc/source/analyzing/analysis_modules/PPVCube.ipynb
+++ b/doc/source/analyzing/analysis_modules/PPVCube.ipynb
@@ -1,451 +1,421 @@
{
- "cells": [
+ "metadata": {
+ "name": "",
+ "signature": "sha256:7ba601e381023e5b7c00a585ccaa55996316d2dfe7f8c96284b4539e1ee5b846"
+ },
+ "nbformat": 3,
+ "nbformat_minor": 0,
+ "worksheets": [
{
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Detailed spectra of astrophysical objects sometimes allow for determinations of how much of the gas is moving with a certain velocity along the line of sight, thanks to Doppler shifting of spectral lines. This enables \"data cubes\" to be created in RA, Dec, and line-of-sight velocity space. In yt, we can use the `PPVCube` analysis module to project fields along a given line of sight traveling at different line-of-sight velocities, to \"mock-up\" what would be seen in observations."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "import yt\n",
- "import numpy as np\n",
- "\n",
- "from yt.analysis_modules.ppv_cube.api import PPVCube"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "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:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Density: $\\rho(r) \\propto r^{\\alpha}$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Velocity: $v_{\\theta}(r) \\propto \\frac{r}{1+(r/r_0)^{\\beta}}$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "where for simplicity we won't worry about the normalizations of these profiles. "
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "First, we'll set up the grid and the parameters of the profiles:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "nx,ny,nz = (256,256,256) # domain dimensions\n",
- "R = 10. # outer radius of disk, kpc\n",
- "r_0 = 3. # scale radius, kpc\n",
- "beta = 1.4 # for the tangential velocity profile\n",
- "alpha = -1. # for the radial density profile\n",
- "x, y = np.mgrid[-R:R:nx*1j,-R:R:ny*1j] # cartesian coordinates of x-y plane of disk\n",
- "r = np.sqrt(x*x+y*y) # polar coordinates\n",
- "theta = np.arctan2(y, x) # polar coordinates"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "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. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "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",
- "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"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Finally, we'll package these data arrays up into a dictionary, which will then be shipped off to `load_uniform_grid`. We'll define the width of the grid to be `2*R` kpc, which will be equal to 1 `code_length`. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "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",
- "bbox = np.array([[-0.5,0.5],[-0.5,0.5],[-0.5,0.5]]) # bbox of width 1 on a side with center (0,0,0)\n",
- "ds = yt.load_uniform_grid(data, (nx,ny,nz), length_unit=(2*R,\"kpc\"), nprocs=1, bbox=bbox)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "To get a sense of what the data looks like, we'll take a slice through the middle of the disk:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "slc = yt.SlicePlot(ds, \"z\", [\"density\",\"velocity_x\",\"velocity_y\",\"velocity_magnitude\"])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "slc.set_log(\"velocity_x\", False)\n",
- "slc.set_log(\"velocity_y\", False)\n",
- "slc.set_log(\"velocity_magnitude\", False)\n",
- "slc.set_unit(\"velocity_magnitude\", \"km/s\")\n",
- "slc.show()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Which shows a rotating disk with a specific density and velocity profile. Now, suppose we wanted to look at this disk galaxy from a certain orientation angle, and simulate a 3D FITS data cube where we can see the gas that is emitting at different velocities along the line of sight. We can do this using the `PPVCube` class. First, let's assume we rotate our viewing angle 60 degrees from face-on, from along the z-axis into the y-axis. We'll create a normal vector:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "i = 60.*np.pi/180.\n",
- "L = [0.0,np.sin(i),np.sin(i)]"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Next, we need to specify a field that will serve as the \"intensity\" of the emission that we see. For simplicity, we'll simply choose the gas density as this field, though it could be any field (including derived fields) in principle. We also need to specify the dimensions of the data cube, and optionally we may choose the bounds in line-of-sight velocity that the data will be binned into. Otherwise, the bounds will simply be set to the negative and positive of the largest speed in the dataset."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "cube = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-150.,150.,\"km/s\"), method=\"sum\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "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",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "cube.write_fits(\"cube.fits\", clobber=True, length_unit=\"kpc\")"
- ]
- },
- {
- "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",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "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)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Now, we'll look at the FITS dataset in yt and look at different slices along the velocity axis, which is the \"z\" axis:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "ds_cube = yt.load(\"cube.fits\")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "# Specifying no center gives us the center slice\n",
- "slc = yt.SlicePlot(ds_cube, \"z\", [\"density\"])\n",
- "slc.show()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "import yt.units as u\n",
- "# Picking different velocities for the slices\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()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "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()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "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()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "If we project all the emission at all the different velocities along the z-axis, we recover the entire disk:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "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()"
- ]
- },
- {
- "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",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "cube2 = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-150,150,\"km/s\"), thermal_broad=True, \n",
- " atomic_weight=12.0, method=\"sum\")\n",
- "cube2.write_fits(\"cube2.fits\", clobber=True, length_unit=\"kpc\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Taking a slice of this cube shows:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "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()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "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()"
- ]
- },
- {
- "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",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "print cube2.vbins[0], cube2.vbins[-1]\n",
- "cube2.transform_spectral_axis(400.0,\"nm\")\n",
- "print cube2.vbins[0], cube2.vbins[-1]"
- ]
- },
- {
- "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",
- "execution_count": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "cube2.reset_spectral_axis()\n",
- "print cube2.vbins[0], cube2.vbins[-1]"
- ]
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Detailed spectra of astrophysical objects sometimes allow for determinations of how much of the gas is moving with a certain velocity along the line of sight, thanks to Doppler shifting of spectral lines. This enables \"data cubes\" to be created in RA, Dec, and line-of-sight velocity space. In yt, we can use the `PPVCube` analysis module to project fields along a given line of sight traveling at different line-of-sight velocities, to \"mock-up\" what would be seen in observations."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "import yt\n",
+ "import numpy as np\n",
+ "\n",
+ "from yt.analysis_modules.ppv_cube.api import PPVCube"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "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:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Density: $\\rho(r) \\propto r^{\\alpha}$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Velocity: $v_{\\theta}(r) \\propto \\frac{r}{1+(r/r_0)^{\\beta}}$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where for simplicity we won't worry about the normalizations of these profiles. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First, we'll set up the grid and the parameters of the profiles:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "nx,ny,nz = (256,256,256) # domain dimensions\n",
+ "R = 10. # outer radius of disk, kpc\n",
+ "r_0 = 3. # scale radius, kpc\n",
+ "beta = 1.4 # for the tangential velocity profile\n",
+ "alpha = -1. # for the radial density profile\n",
+ "x, y = np.mgrid[-R:R:nx*1j,-R:R:ny*1j] # cartesian coordinates of x-y plane of disk\n",
+ "r = np.sqrt(x*x+y*y) # polar coordinates\n",
+ "theta = np.arctan2(y, x) # polar coordinates"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "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. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "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",
+ "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"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Finally, we'll package these data arrays up into a dictionary, which will then be shipped off to `load_uniform_grid`. We'll define the width of the grid to be `2*R` kpc, which will be equal to 1 `code_length`. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "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",
+ "bbox = np.array([[-0.5,0.5],[-0.5,0.5],[-0.5,0.5]]) # bbox of width 1 on a side with center (0,0,0)\n",
+ "ds = yt.load_uniform_grid(data, (nx,ny,nz), length_unit=(2*R,\"kpc\"), nprocs=1, bbox=bbox)"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To get a sense of what the data looks like, we'll take a slice through the middle of the disk:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "slc = yt.SlicePlot(ds, \"z\", [\"density\",\"velocity_x\",\"velocity_y\",\"velocity_magnitude\"])"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "slc.set_log(\"velocity_x\", False)\n",
+ "slc.set_log(\"velocity_y\", False)\n",
+ "slc.set_log(\"velocity_magnitude\", False)\n",
+ "slc.set_unit(\"velocity_magnitude\", \"km/s\")\n",
+ "slc.show()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Which shows a rotating disk with a specific density and velocity profile. Now, suppose we wanted to look at this disk galaxy from a certain orientation angle, and simulate a 3D FITS data cube where we can see the gas that is emitting at different velocities along the line of sight. We can do this using the `PPVCube` class. First, let's assume we rotate our viewing angle 60 degrees from face-on, from along the z-axis into the y-axis. We'll create a normal vector:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "i = 60.*np.pi/180.\n",
+ "L = [0.0,np.sin(i),np.sin(i)]"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next, we need to specify a field that will serve as the \"intensity\" of the emission that we see. For simplicity, we'll simply choose the gas density as this field, though it could be any field (including derived fields) in principle. We also need to specify the dimensions of the data cube, and optionally we may choose the bounds in line-of-sight velocity that the data will be binned into. Otherwise, the bounds will simply be set to the negative and positive of the largest speed in the dataset."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "cube = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-150.,150.,\"km/s\"), method=\"sum\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "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=\"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": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now, we'll look at the FITS dataset in yt and look at different slices along the velocity axis, which is the \"z\" axis:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "ds_cube = yt.load(\"cube.fits\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "# Specifying no center gives us the center slice\n",
+ "slc = yt.SlicePlot(ds_cube, \"z\", [\"density\"])\n",
+ "slc.show()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "import yt.units as u\n",
+ "# Picking different velocities for the slices\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",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "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",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "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",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "If we project all the emission at all the different velocities along the z-axis, we recover the entire disk:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "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()"
+ ],
+ "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, method=\"sum\")\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": {}
}
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "IPython (Python 2)",
- "name": "python2"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 2
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython2"
- },
- "signature": "sha256:c0d10b9e420fd80ed69d1ad4dcfb9eb852999c7bf8524868f2df19e940e1fae8"
- },
- "nbformat": 4,
- "nbformat_minor": 0
+ ]
}
\ No newline at end of file
https://bitbucket.org/yt_analysis/yt/commits/6e0ee61ccebb/
Changeset: 6e0ee61ccebb
Branch: yt
User: jzuhone
Date: 2014-12-02 22:16:45+00:00
Summary: Tuples for field names. Making the profile plots look nicer.
Affected #: 1 file
diff -r afb117e019cc3e1a6b525c7d216f8a5ce312055a -r 6e0ee61ccebb9186c2f572e76fca9db2e1e10423 doc/source/cookbook/fits_xray_images.ipynb
--- a/doc/source/cookbook/fits_xray_images.ipynb
+++ b/doc/source/cookbook/fits_xray_images.ipynb
@@ -1,7 +1,7 @@
{
"metadata": {
"name": "",
- "signature": "sha256:650e3fc7f66951a5fcdb18332bbc625f6f6e449fc919acd01da01e1fbbf92ee1"
+ "signature": "sha256:73eb48af6c5628619b5cab88840f6ac2b3eb278ae057c34d567f4a127181b8fb"
},
"nbformat": 3,
"nbformat_minor": 0,
@@ -74,15 +74,15 @@
"def _counts(field, data):\n",
" exposure_time = data.get_field_parameter(\"exposure_time\")\n",
" return data[\"flux\"]*data[\"pixel\"]*exposure_time\n",
- "ds.add_field(name=\"counts\", function=_counts, units=\"counts\", take_log=False)\n",
+ "ds.add_field((\"gas\",\"counts\"), function=_counts, units=\"counts\", take_log=False)\n",
"\n",
"def _pp(field, data):\n",
" return np.sqrt(data[\"counts\"])*data[\"projected_temperature\"]\n",
- "ds.add_field(name=\"pseudo_pressure\", function=_pp, units=\"sqrt(counts)*keV\", take_log=False)\n",
+ "ds.add_field((\"gas\",\"pseudo_pressure\"), function=_pp, units=\"sqrt(counts)*keV\", take_log=False)\n",
"\n",
"def _pe(field, data):\n",
" return data[\"projected_temperature\"]*data[\"counts\"]**(-1./3.)\n",
- "ds.add_field(name=\"pseudo_entropy\", function=_pe, units=\"keV*(counts)**(-1/3)\", take_log=False)"
+ "ds.add_field((\"gas\",\"pseudo_entropy\"), function=_pe, units=\"keV*(counts)**(-1/3)\", take_log=False)"
],
"language": "python",
"metadata": {},
@@ -180,10 +180,11 @@
"input": [
"radial_profile = yt.ProfilePlot(my_sphere, \"radius\", \n",
" [\"counts\",\"pseudo_pressure\",\"pseudo_entropy\"], \n",
- " n_bins=50, weight_field=\"ones\")\n",
+ " n_bins=30, weight_field=\"ones\")\n",
"radial_profile.set_log(\"counts\", True)\n",
"radial_profile.set_log(\"pseudo_pressure\", True)\n",
"radial_profile.set_log(\"pseudo_entropy\", True)\n",
+ "radial_profile.set_xlim(3,100.)\n",
"radial_profile.show()"
],
"language": "python",
@@ -401,4 +402,4 @@
"metadata": {}
}
]
-}
+}
\ No newline at end of file
https://bitbucket.org/yt_analysis/yt/commits/9b5859482dcd/
Changeset: 9b5859482dcd
Branch: yt
User: jzuhone
Date: 2014-12-02 22:17:11+00:00
Summary: These should have been raw strings
Affected #: 1 file
diff -r 6e0ee61ccebb9186c2f572e76fca9db2e1e10423 -r 9b5859482dcd2c9334bd56caac6ba89af84aa3c0 yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -775,7 +775,7 @@
"%s_axis" % ("xy"[i]))[axis_index]
unn = self.ds.coordinates.default_unit_label.get(axax, "")
if unn != "":
- axes_unit_labels[i] = '\/\/\left('+unn+'\right)'
+ axes_unit_labels[i] = r'\/\/\left('+unn+r'\right)'
continue
# Use sympy to factor h out of the unit. In this context 'un'
# is a string, so we call the Unit constructor.
https://bitbucket.org/yt_analysis/yt/commits/ee4276f9be71/
Changeset: ee4276f9be71
Branch: yt
User: xarthisius
Date: 2014-12-03 17:43:19+00:00
Summary: Merged in jzuhone/yt (pull request #1326)
[bugfix] Fixing problems in docs build, adding some keywords back to FITSImageBuffer
Affected #: 6 files
diff -r 02ef23fcc82e96c58e4786027e42534cc6d09be1 -r ee4276f9be712f57aa6dc829d7990fe6f8511f83 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:b83e125278c2e58da4d99ac9d2ca2a136d01f1094e1b83497925e0f9b9b056c2"
+ "signature": "sha256:7ba601e381023e5b7c00a585ccaa55996316d2dfe7f8c96284b4539e1ee5b846"
},
"nbformat": 3,
"nbformat_minor": 0,
@@ -193,7 +193,7 @@
"cell_type": "code",
"collapsed": false,
"input": [
- "cube = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-150.,150.,\"km/s\"))"
+ "cube = PPVCube(ds, L, \"density\", dims=(200,100,50), velocity_bounds=(-150.,150.,\"km/s\"), method=\"sum\")"
],
"language": "python",
"metadata": {},
@@ -335,7 +335,7 @@
"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",
+ " atomic_weight=12.0, method=\"sum\")\n",
"cube2.write_fits(\"cube2.fits\", clobber=True, length_unit=\"kpc\")"
],
"language": "python",
diff -r 02ef23fcc82e96c58e4786027e42534cc6d09be1 -r ee4276f9be712f57aa6dc829d7990fe6f8511f83 doc/source/cookbook/fits_xray_images.ipynb
--- a/doc/source/cookbook/fits_xray_images.ipynb
+++ b/doc/source/cookbook/fits_xray_images.ipynb
@@ -1,7 +1,7 @@
{
"metadata": {
"name": "",
- "signature": "sha256:650e3fc7f66951a5fcdb18332bbc625f6f6e449fc919acd01da01e1fbbf92ee1"
+ "signature": "sha256:73eb48af6c5628619b5cab88840f6ac2b3eb278ae057c34d567f4a127181b8fb"
},
"nbformat": 3,
"nbformat_minor": 0,
@@ -74,15 +74,15 @@
"def _counts(field, data):\n",
" exposure_time = data.get_field_parameter(\"exposure_time\")\n",
" return data[\"flux\"]*data[\"pixel\"]*exposure_time\n",
- "ds.add_field(name=\"counts\", function=_counts, units=\"counts\", take_log=False)\n",
+ "ds.add_field((\"gas\",\"counts\"), function=_counts, units=\"counts\", take_log=False)\n",
"\n",
"def _pp(field, data):\n",
" return np.sqrt(data[\"counts\"])*data[\"projected_temperature\"]\n",
- "ds.add_field(name=\"pseudo_pressure\", function=_pp, units=\"sqrt(counts)*keV\", take_log=False)\n",
+ "ds.add_field((\"gas\",\"pseudo_pressure\"), function=_pp, units=\"sqrt(counts)*keV\", take_log=False)\n",
"\n",
"def _pe(field, data):\n",
" return data[\"projected_temperature\"]*data[\"counts\"]**(-1./3.)\n",
- "ds.add_field(name=\"pseudo_entropy\", function=_pe, units=\"keV*(counts)**(-1/3)\", take_log=False)"
+ "ds.add_field((\"gas\",\"pseudo_entropy\"), function=_pe, units=\"keV*(counts)**(-1/3)\", take_log=False)"
],
"language": "python",
"metadata": {},
@@ -180,10 +180,11 @@
"input": [
"radial_profile = yt.ProfilePlot(my_sphere, \"radius\", \n",
" [\"counts\",\"pseudo_pressure\",\"pseudo_entropy\"], \n",
- " n_bins=50, weight_field=\"ones\")\n",
+ " n_bins=30, weight_field=\"ones\")\n",
"radial_profile.set_log(\"counts\", True)\n",
"radial_profile.set_log(\"pseudo_pressure\", True)\n",
"radial_profile.set_log(\"pseudo_entropy\", True)\n",
+ "radial_profile.set_xlim(3,100.)\n",
"radial_profile.show()"
],
"language": "python",
@@ -401,4 +402,4 @@
"metadata": {}
}
]
-}
+}
\ No newline at end of file
diff -r 02ef23fcc82e96c58e4786027e42534cc6d09be1 -r ee4276f9be712f57aa6dc829d7990fe6f8511f83 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,7 +13,8 @@
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, sanitize_fits_unit
+from yt.utilities.fits_image import FITSImageBuffer, sanitize_fits_unit, \
+ create_sky_wcs
from yt.visualization.volume_rendering.camera import off_axis_projection
from yt.funcs import get_pbar
from yt.utilities.physical_constants import clight, mh
@@ -188,6 +189,9 @@
else:
self.width = ds.quan(self.width, "code_length")
+ self.ds.field_info.pop(("gas","intensity"))
+ self.ds.field_info.pop(("gas","v_los"))
+
def create_intensity(self):
def _intensity(field, data):
v = self.current_v-data["v_los"].v
@@ -245,57 +249,40 @@
Whether or not to clobber an existing file with the same name.
length_unit : string
The units to convert the coordinates to in the file.
- sky_scale : tuple or YTQuantity
+ sky_scale : tuple, optional
Conversion between an angle unit and a length unit, if sky
- coordinates are desired.
- Examples: (1.0, "arcsec/kpc"), YTQuantity(0.001, "deg/kpc")
+ coordinates are desired, e.g. (1.0, "arcsec/kpc")
sky_center : tuple, optional
- The (RA, Dec) coordinate in degrees of the central pixel if
- *sky_scale* has been specified.
+ The (RA, Dec) coordinate in degrees of the central pixel. Must
+ be specified with *sky_scale*.
Examples
--------
>>> cube.write_fits("my_cube.fits", clobber=False, sky_scale=(1.0,"arcsec/kpc"))
"""
- 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"]
-
vunit = fits_info[self.axis_type][0]
vtype = fits_info[self.axis_type][1]
v_center = 0.5*(self.vbins[0]+self.vbins[-1]).in_units(vunit).value
- if sky_scale:
- dx = (self.width*sky_scale).in_units("deg").v/self.nx
- units = "deg"
+ if length_unit is None:
+ units = str(self.ds.get_smallest_appropriate_unit(self.width))
else:
- if length_unit is None:
- units = str(self.ds.get_smallest_appropriate_unit(self.width))
- else:
- units = length_unit
+ units = length_unit
units = sanitize_fits_unit(units)
dx = self.width.in_units(units).v/self.nx
- dy = dx
+ dy = self.width.in_units(units).v/self.ny
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.crval = [0.0,0.0,v_center]
w.wcs.cunit = [units,units,vunit]
- w.wcs.ctype = [types[0],types[1],vtype]
+ w.wcs.ctype = ["LINEAR","LINEAR",vtype]
+
+ if sky_scale is not None and sky_center is not None:
+ w = create_sky_wcs(w, sky_center, sky_scale)
fib = FITSImageBuffer(self.data.transpose(2,0,1), fields=self.field, wcs=w)
fib[0].header["bunit"] = re.sub('()', '', str(self.proj_units))
diff -r 02ef23fcc82e96c58e4786027e42534cc6d09be1 -r ee4276f9be712f57aa6dc829d7990fe6f8511f83 yt/analysis_modules/sunyaev_zeldovich/projection.py
--- a/yt/analysis_modules/sunyaev_zeldovich/projection.py
+++ b/yt/analysis_modules/sunyaev_zeldovich/projection.py
@@ -298,13 +298,12 @@
----------
filename : string
The name of the FITS file to be written.
- sky_scale : tuple or YTQuantity
+ sky_scale : tuple
Conversion between an angle unit and a length unit, if sky
- coordinates are desired.
- Examples: (1.0, "arcsec/kpc"), YTQuantity(0.001, "deg/kpc")
+ coordinates are desired, e.g. (1.0, "arcsec/kpc")
sky_center : tuple, optional
- The (RA, Dec) coordinate in degrees of the central pixel if
- *sky_scale* has been specified.
+ The (RA, Dec) coordinate in degrees of the central pixel. Must
+ be specified with *sky_scale*.
clobber : boolean, optional
If the file already exists, do we overwrite?
@@ -314,40 +313,23 @@
>>> szprj.write_fits("SZbullet.fits", clobber=False)
>>> # This example uses sky coords
>>> sky_scale = (1., "arcsec/kpc") # One arcsec per kpc
- >>> sky_center = (30., 45.) # In degrees
+ >>> sky_center = (30., 45., "deg")
>>> szprj.write_fits("SZbullet.fits", sky_center=sky_center, sky_scale=sky_scale)
"""
- from yt.utilities.fits_image import FITSImageBuffer, sanitize_fits_unit
+ from yt.utilities.fits_image import FITSImageBuffer, create_sky_wcs
- 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"]
-
- 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"
+ dx = self.dx.in_units("kpc")
dy = dx
- if sky_scale:
- dx *= -1.
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
+ w.wcs.crval = [0.0,0.0]
+ w.wcs.cunit = ["kpc"]*2
+ w.wcs.ctype = ["LINEAR"]*2
+
+ if sky_scale is not None and sky_center is not None:
+ w = create_sky_wcs(w, sky_center, sky_scale)
fib = FITSImageBuffer(self.data, fields=self.data.keys(), wcs=w)
fib.writeto(filename, clobber=clobber)
diff -r 02ef23fcc82e96c58e4786027e42534cc6d09be1 -r ee4276f9be712f57aa6dc829d7990fe6f8511f83 yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -28,7 +28,7 @@
class FITSImageBuffer(HDUList):
- def __init__(self, data, fields=None, units="cm", wcs=None):
+ def __init__(self, data, fields=None, units=None, pixel_scale=None, wcs=None):
r""" Initialize a FITSImageBuffer object.
FITSImageBuffer contains a list of FITS ImageHDU instances, and
@@ -49,8 +49,11 @@
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. Only applies
- to FixedResolutionBuffers or YTCoveringGrids. Defaults to "cm".
+ The units of the WCS coordinates. Defaults to "cm".
+ pixel_scale : float
+ The scale of the pixel, in *units*. Either a single float or
+ iterable of floats. Only used if this information is not already
+ provided by *data*.
wcs : `astropy.wcs.WCS` instance, optional
Supply an AstroPy WCS instance. Will override automatic WCS
creation from FixedResolutionBuffers and YTCoveringGrids.
@@ -77,7 +80,10 @@
>>> f_deg = FITSImageBuffer(frb, fields="kT", wcs=w)
>>> f_deg.writeto("temp.fits")
"""
-
+
+ if units is None: units = "cm"
+ if pixel_scale is None: pixel_scale = 1.0
+
super(FITSImageBuffer, self).__init__()
if isinstance(fields, basestring): fields = [fields]
@@ -91,13 +97,13 @@
img_data = {}
if fields is None:
mylog.error("Please specify a field name for this array.")
- raise KeyError
+ raise KeyError("Please specify a field name for this array.")
img_data[fields[0]] = data
if fields is None: fields = img_data.keys()
if len(fields) == 0:
mylog.error("Please specify one or more fields to write.")
- raise KeyError
+ raise KeyError("Please specify one or more fields to write.")
first = True
@@ -128,12 +134,8 @@
elif self.dimensionality == 3:
self.nx, self.ny, self.nz = self[0].data.shape
- has_coords = (isinstance(img_data, FixedResolutionBuffer) or
- isinstance(img_data, YTCoveringGridBase))
-
if wcs is None:
w = pywcs.WCS(header=self[0].header, naxis=self.dimensionality)
- w.wcs.crpix = 0.5*(np.array(self.shape)+1)
if isinstance(img_data, FixedResolutionBuffer):
# FRBs are a special case where we have coordinate
# information, so we take advantage of this and
@@ -143,26 +145,28 @@
xctr = 0.5*(img_data.bounds[1]+img_data.bounds[0]).in_units(units)
yctr = 0.5*(img_data.bounds[3]+img_data.bounds[2]).in_units(units)
center = [xctr, yctr]
+ cdelt = [dx,dy]
elif isinstance(img_data, YTCoveringGridBase):
- dx, dy, dz = img_data.dds.in_units(units)
+ cdelt = img_data.dds.in_units(units).v
center = 0.5*(img_data.left_edge+img_data.right_edge).in_units(units)
else:
- # 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)
+ # If img_data is just an array, we assume the center is the origin
+ # and use *pixel_scale* to determine the cell widths
+ if iterable(pixel_scale):
+ cdelt = pixel_scale
+ else:
+ cdelt = [pixel_scale]*self.dimensionality
+ center = [0.0]*self.dimensionality
+ w.wcs.crpix = 0.5*(np.array(self.shape)+1)
w.wcs.crval = center
- 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.cdelt = cdelt
w.wcs.ctype = ["linear"]*self.dimensionality
- self._set_wcs(w)
+ w.wcs.cunit = [units]*self.dimensionality
+ self.set_wcs(w)
else:
- self._set_wcs(wcs)
+ self.set_wcs(wcs)
- def _set_wcs(self, wcs):
+ def set_wcs(self, wcs):
"""
Set the WCS coordinate information for all images
with a WCS object *wcs*.
@@ -250,6 +254,50 @@
axis_wcs = [[1,2],[0,2],[0,1]]
+def create_sky_wcs(old_wcs, sky_center, sky_scale,
+ ctype=["RA---TAN","DEC--TAN"], crota=None):
+ """
+ Takes an astropy.wcs.WCS instance created in yt from a
+ simulation that has a Cartesian coordinate system and
+ converts it to one in a celestial coordinate system.
+
+ Parameters
+ ----------
+ old_wcs : astropy.wcs.WCS
+ The original WCS to be converted.
+ sky_center : tuple
+ Reference coordinates of the WCS in degrees.
+ sky_scale : tuple
+ Conversion between an angle unit and a length unit,
+ e.g. (3.0, "arcsec/kpc")
+ ctype : list of strings, optional
+ The type of the coordinate system to create.
+ crota : list of floats, optional
+ Rotation angles between cartesian coordinates and
+ the celestial coordinates.
+ """
+ naxis = old_wcs.naxis
+ crval = [sky_center[0], sky_center[1]]
+ scaleq = YTQuantity(sky_scale[0],sky_scale[1])
+ deltas = old_wcs.wcs.cdelt
+ units = [str(unit) for unit in old_wcs.wcs.cunit]
+ new_dx = (YTQuantity(-deltas[0], units[0])*scaleq).in_units("deg")
+ new_dy = (YTQuantity(deltas[1], units[1])*scaleq).in_units("deg")
+ new_wcs = pywcs.WCS(naxis=naxis)
+ cdelt = [new_dx.v, new_dy.v]
+ if naxis == 3:
+ crval.append(old_wcs.wcs.crval[2])
+ cdelt.append(old_wcs.wcs.cdelt[2])
+ ctype.append(old_wcs.wcs.ctype[2])
+ new_wcs.wcs.crpix = old_wcs.wcs.crpix
+ new_wcs.wcs.cdelt = cdelt
+ new_wcs.wcs.crval = crval
+ new_wcs.wcs.cunit = ["deg"]*naxis
+ new_wcs.wcs.ctype = ctype
+ if crota is not None:
+ new_wcs.wcs.crota = crota
+ return new_wcs
+
def sanitize_fits_unit(unit):
if unit == "Mpc":
mylog.info("Changing FITS file unit to kpc.")
diff -r 02ef23fcc82e96c58e4786027e42534cc6d09be1 -r ee4276f9be712f57aa6dc829d7990fe6f8511f83 yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -775,7 +775,7 @@
"%s_axis" % ("xy"[i]))[axis_index]
unn = self.ds.coordinates.default_unit_label.get(axax, "")
if unn != "":
- axes_unit_labels[i] = '\/\/\left('+unn+'\right)'
+ axes_unit_labels[i] = r'\/\/\left('+unn+r'\right)'
continue
# Use sympy to factor h out of the unit. In this context 'un'
# is a string, so we call the Unit constructor.
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