[yt-svn] commit/yt: xarthisius: Merged in jzuhone/yt (pull request #1603)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Jun 18 06:17:38 PDT 2015
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/e4bc6b4f239b/
Changeset: e4bc6b4f239b
Branch: yt
User: xarthisius
Date: 2015-06-18 13:17:28+00:00
Summary: Merged in jzuhone/yt (pull request #1603)
FITS image writing refactor
Affected #: 12 files
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f doc/source/visualizing/FITSImageBuffer.ipynb
--- a/doc/source/visualizing/FITSImageBuffer.ipynb
+++ /dev/null
@@ -1,205 +0,0 @@
-{
- "metadata": {
- "name": "",
- "signature": "sha256:872f7525edd3c1ee09c67f6ecdd8552218df05ebe5ab73bcab55654edf0ac2bb"
- },
- "nbformat": 3,
- "nbformat_minor": 0,
- "worksheets": [
- {
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "yt has capabilities for writing 2D and 3D uniformly gridded data generated from datasets to FITS files. This is via the `FITSImageBuffer` class, which has subclasses `FITSSlice` and `FITSProjection` to write slices and projections directly to FITS. We'll test this out on an Athena dataset."
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "%matplotlib inline\n",
- "import yt\n",
- "from yt.utilities.fits_image import FITSImageBuffer, FITSSlice, FITSProjection"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "ds = yt.load(\"MHDSloshing/virgo_low_res.0054.vtk\", parameters={\"length_unit\":(1.0,\"Mpc\"),\n",
- " \"mass_unit\":(1.0e14,\"Msun\"),\n",
- " \"time_unit\":(1.0,\"Myr\")})"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "To demonstrate a useful example of creating a FITS file, let's first make a `ProjectionPlot`:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "prj = yt.ProjectionPlot(ds, \"z\", [\"temperature\"], weight_field=\"density\", width=(500.,\"kpc\"))\n",
- "prj.show()"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Suppose that we wanted to write this projection to a FITS file for analysis and visualization in other programs, such as ds9. We can do that using `FITSProjection`:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "prj_fits = FITSProjection(ds, \"z\", [\"temperature\"], weight_field=\"density\")"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which took the same parameters as `ProjectionPlot` except the width, because `FITSProjection` and `FITSSlice` always make slices and projections of the width of the domain size, at the finest resolution available in the simulation, in a unit determined to be appropriate for the physical size of the dataset. `prj_fits` is a full-fledged FITS file in memory, specifically an [AstroPy `HDUList`](http://astropy.readthedocs.org/en/latest/io/fits/api/hdulists.html) object. This means that we can use all of the methods inherited from `HDUList`:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "prj_fits.info()"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "`info` shows us the contents of the virtual FITS file. We can also look at the header for the `\"temperature\"` image, like so:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "prj_fits[\"temperature\"].header"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "where we can see that the temperature units are in Kelvin and the cell widths are in kiloparsecs. The projection can be written to disk using the `writeto` method:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "prj_fits.writeto(\"sloshing.fits\", clobber=True)"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Since yt can read FITS image files, it can be loaded up just like any other dataset:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "ds2 = yt.load(\"sloshing.fits\")"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "and we can make a `SlicePlot` of the 2D image, which shows the same data as the previous image:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "slc2 = yt.SlicePlot(ds2, \"z\", [\"temperature\"], width=(500.,\"kpc\"))\n",
- "slc2.set_log(\"temperature\", True)\n",
- "slc2.show()"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "If you want more fine-grained control over what goes into the FITS file, you can call `FITSImageBuffer` directly, with various kinds of inputs. For example, you could use a `FixedResolutionBuffer`, and specify you want the units in parsecs instead:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "slc3 = ds.slice(0, 0.0)\n",
- "frb = slc3.to_frb((500.,\"kpc\"), 800)\n",
- "fib = FITSImageBuffer(frb, fields=[\"density\",\"temperature\"], units=\"pc\")"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Finally, a 3D FITS cube can be created from a covering grid:"
- ]
- },
- {
- "cell_type": "code",
- "collapsed": false,
- "input": [
- "cvg = ds.covering_grid(ds.index.max_level, [-0.5,-0.5,-0.5], [64, 64, 64], fields=[\"density\",\"temperature\"])\n",
- "fib = FITSImageBuffer(cvg, fields=[\"density\",\"temperature\"], units=\"Mpc\")"
- ],
- "language": "python",
- "metadata": {},
- "outputs": []
- }
- ],
- "metadata": {}
- }
- ]
-}
\ No newline at end of file
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f doc/source/visualizing/FITSImageData.ipynb
--- /dev/null
+++ b/doc/source/visualizing/FITSImageData.ipynb
@@ -0,0 +1,409 @@
+{
+ "metadata": {
+ "name": "",
+ "signature": "sha256:c7de5ef190feaa2289595aec7eaa05db02fd535e408e0d04aa54088b0bd3ebae"
+ },
+ "nbformat": 3,
+ "nbformat_minor": 0,
+ "worksheets": [
+ {
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "yt has capabilities for writing 2D and 3D uniformly gridded data generated from datasets to FITS files. This is via the `FITSImageData` class. We'll test these capabilities out on an Athena dataset."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "import yt\n",
+ "from yt.utilities.fits_image import FITSImageData, FITSProjection"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "ds = yt.load(\"MHDSloshing/virgo_low_res.0054.vtk\", parameters={\"length_unit\":(1.0,\"Mpc\"),\n",
+ " \"mass_unit\":(1.0e14,\"Msun\"),\n",
+ " \"time_unit\":(1.0,\"Myr\")})"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "heading",
+ "level": 2,
+ "metadata": {},
+ "source": [
+ "Creating FITS images from Slices and Projections"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "There are several ways to make a `FITSImageData` instance. The most intuitive ways are to use the `FITSSlice`, `FITSProjection`, `FITSOffAxisSlice`, and `FITSOffAxisProjection` classes to write slices and projections directly to FITS. To demonstrate a useful example of creating a FITS file, let's first make a `ProjectionPlot`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj = yt.ProjectionPlot(ds, \"z\", [\"temperature\"], weight_field=\"density\", width=(500.,\"kpc\"))\n",
+ "prj.show()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Suppose that we wanted to write this projection to a FITS file for analysis and visualization in other programs, such as ds9. We can do that using `FITSProjection`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits = FITSProjection(ds, \"z\", [\"temperature\"], weight_field=\"density\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "which took the same parameters as `ProjectionPlot` except the width, because `FITSProjection` and `FITSSlice` always make slices and projections of the width of the domain size, at the finest resolution available in the simulation, in a unit determined to be appropriate for the physical size of the dataset."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Because `FITSImageData` inherits from the [AstroPy `HDUList`](http://astropy.readthedocs.org/en/latest/io/fits/api/hdulists.html) class, we can call its methods. For example, `info` shows us the contents of the virtual FITS file:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits.info()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can also look at the header for a particular field:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits[\"temperature\"].header"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where we can see that the temperature units are in Kelvin and the cell widths are in kiloparsecs. If we want the raw image data with units, we can call `get_data`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits.get_data(\"temperature\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can use the `set_unit` method to change the units of a particular field:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits.set_unit(\"temperature\",\"R\")\n",
+ "prj_fits.get_data(\"temperature\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The image can be written to disk using the `writeto` method:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits.writeto(\"sloshing.fits\", clobber=True)"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Since yt can read FITS image files, it can be loaded up just like any other dataset:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "ds2 = yt.load(\"sloshing.fits\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and we can make a `SlicePlot` of the 2D image, which shows the same data as the previous image:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "slc2 = yt.SlicePlot(ds2, \"z\", [\"temperature\"], width=(500.,\"kpc\"))\n",
+ "slc2.set_log(\"temperature\", True)\n",
+ "slc2.show()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "heading",
+ "level": 2,
+ "metadata": {},
+ "source": [
+ "Using `FITSImageData` directly"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "If you want more fine-grained control over what goes into the FITS file, you can call `FITSImageData` directly, with various kinds of inputs. For example, you could use a `FixedResolutionBuffer`, and specify you want the units in parsecs instead:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "slc3 = ds.slice(0, 0.0)\n",
+ "frb = slc3.to_frb((500.,\"kpc\"), 800)\n",
+ "fid_frb = FITSImageData(frb, fields=[\"density\",\"temperature\"], units=\"pc\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "A 3D FITS cube can also be created from a covering grid:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "cvg = ds.covering_grid(ds.index.max_level, [-0.5,-0.5,-0.5], [64, 64, 64], fields=[\"density\",\"temperature\"])\n",
+ "fid_cvg = FITSImageData(cvg, fields=[\"density\",\"temperature\"], units=\"Mpc\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "heading",
+ "level": 2,
+ "metadata": {},
+ "source": [
+ "Other `FITSImageData` Methods"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "A `FITSImageData` instance can be generated from one previously written to disk using the `from_file` classmethod:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "fid = FITSImageData.from_file(\"sloshing.fits\")\n",
+ "fid.info()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Multiple `FITSImageData` can be combined to create a new one, provided that the coordinate information is the same:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits2 = FITSProjection(ds, \"z\", [\"density\"])\n",
+ "prj_fits3 = FITSImageData.from_images([prj_fits, prj_fits2])\n",
+ "prj_fits3.info()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Alternatively, individual fields can be popped as well:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "dens_fits = prj_fits3.pop(\"density\")"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "dens_fits.info()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits3.info()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "So far, the FITS images we have shown have linear spatial coordinates. One may want to take a projection of an object and make a crude mock observation out of it, with celestial coordinates. For this, we can use the `create_sky_wcs` method. Specify a center (RA, Dec) coordinate in degrees, as well as a linear scale in terms of angle per distance:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "sky_center = [30.,45.] # in degrees\n",
+ "sky_scale = (2.5, \"arcsec/kpc\") # could also use a YTQuantity\n",
+ "prj_fits.create_sky_wcs(sky_center, sky_scale, ctype=[\"RA---TAN\",\"DEC--TAN\"])"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "By the default, a tangent RA/Dec projection is used, but one could also use another projection using the `ctype` keyword. We can now look at the header and see it has the appropriate WCS:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits[\"temperature\"].header"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Finally, we can add header keywords to a single field or for all fields in the FITS image using `update_header`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "fid_frb.update_header(\"all\", \"time\", 0.1) # Update all the fields\n",
+ "fid_frb.update_header(\"temperature\", \"scale\", \"Rankine\") # Update just one field"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "print fid_frb[\"density\"].header[\"time\"]\n",
+ "print fid_frb[\"temperature\"].header[\"scale\"]"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ }
+ ],
+ "metadata": {}
+ }
+ ]
+}
\ No newline at end of file
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f doc/source/visualizing/writing_fits_images.rst
--- a/doc/source/visualizing/writing_fits_images.rst
+++ b/doc/source/visualizing/writing_fits_images.rst
@@ -3,4 +3,4 @@
Writing FITS Images
==========================
-.. notebook:: FITSImageBuffer.ipynb
\ No newline at end of file
+.. notebook:: FITSImageData.ipynb
\ No newline at end of file
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f 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,8 +13,7 @@
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, \
- create_sky_wcs
+from yt.utilities.fits_image import FITSImageData, 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
@@ -89,6 +88,8 @@
dims : integer, optional
The spatial resolution of the cube. Implies nx = ny, e.g. the
aspect ratio of the PPVCube's spatial dimensions is 1.
+ thermal_broad : boolean, optional
+ Whether or not to broaden the line using the gas temperature. Default: False.
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).
@@ -152,9 +153,7 @@
"methods are supported in PPVCube.")
dd = ds.all_data()
-
fd = dd._determine_fields(field)[0]
-
self.field_units = ds._get_field_info(fd).units
self.vbins = ds.arr(np.linspace(velocity_bounds[0],
@@ -172,7 +171,7 @@
_vlos = create_vlos(normal, self.no_shifting)
self.ds.add_field(("gas","v_los"), function=_vlos, units="cm/s")
- _intensity = self.create_intensity()
+ _intensity = self._create_intensity()
self.ds.add_field(("gas","intensity"), function=_intensity, units=self.field_units)
if method == "integrate" and weight_field is None:
@@ -214,16 +213,6 @@
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
- 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
@@ -259,17 +248,18 @@
self.dv = self.vbins[1]-self.vbins[0]
@parallel_root_only
- def write_fits(self, filename, clobber=True, length_unit=None,
+ def write_fits(self, filename, clobber=False, length_unit=None,
sky_scale=None, sky_center=None):
r""" Write the PPVCube to a FITS file.
Parameters
----------
filename : string
- The name of the file to write.
- clobber : boolean
- Whether or not to clobber an existing file with the same name.
- length_unit : string
+ The name of the file to write to.
+ clobber : boolean, optional
+ Whether to overwrite a file with the same name that already
+ exists. Default False.
+ length_unit : string, optional
The units to convert the coordinates to in the file.
sky_scale : tuple, optional
Conversion between an angle unit and a length unit, if sky
@@ -280,7 +270,8 @@
Examples
--------
- >>> cube.write_fits("my_cube.fits", clobber=False, sky_scale=(1.0,"arcsec/kpc"))
+ >>> cube.write_fits("my_cube.fits", clobber=False,
+ ... sky_scale=(1.0,"arcsec/kpc"), sky_center=(30.,45.))
"""
vunit = fits_info[self.axis_type][0]
vtype = fits_info[self.axis_type][1]
@@ -303,13 +294,11 @@
w.wcs.cunit = [units,units,vunit]
w.wcs.ctype = ["LINEAR","LINEAR",vtype]
+ fib = FITSImageData(self.data.transpose(), fields=self.field, wcs=w)
+ fib.update_all_headers("bunit", re.sub('()', '', str(self.proj_units)))
+ fib.update_all_headers("btype", self.field)
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))
- fib[0].header["btype"] = self.field
-
+ fib.create_sky_wcs(sky_center, sky_scale)
fib.writeto(filename, clobber=clobber)
def __repr__(self):
@@ -320,3 +309,13 @@
def __getitem__(self, item):
return self.data[item]
+
+ def _create_intensity(self):
+ def _intensity(field, data):
+ v = self.current_v-data["v_los"].in_cgs().v
+ T = (data["temperature"]).in_cgs().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
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f yt/analysis_modules/ppv_cube/tests/test_ppv.py
--- a/yt/analysis_modules/ppv_cube/tests/test_ppv.py
+++ b/yt/analysis_modules/ppv_cube/tests/test_ppv.py
@@ -1,5 +1,5 @@
"""
-Unit test the sunyaev_zeldovich analysis module.
+Unit test the PPVCube analysis module.
"""
#-----------------------------------------------------------------------------
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f yt/analysis_modules/sunyaev_zeldovich/projection.py
--- a/yt/analysis_modules/sunyaev_zeldovich/projection.py
+++ b/yt/analysis_modules/sunyaev_zeldovich/projection.py
@@ -316,7 +316,7 @@
>>> 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, create_sky_wcs
+ from yt.utilities.fits_image import FITSImageData
dx = self.dx.in_units("kpc")
dy = dx
@@ -328,10 +328,9 @@
w.wcs.cunit = ["kpc"]*2
w.wcs.ctype = ["LINEAR"]*2
+ fib = FITSImageData(self.data, fields=self.data.keys(), wcs=w)
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.create_sky_wcs(sky_center, sky_scale)
fib.writeto(filename, clobber=clobber)
@parallel_root_only
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f yt/frontends/fits/misc.py
--- a/yt/frontends/fits/misc.py
+++ b/yt/frontends/fits/misc.py
@@ -18,7 +18,7 @@
from yt.funcs import mylog, get_image_suffix
from yt.visualization._mpl_imports import FigureCanvasAgg
from yt.units.yt_array import YTQuantity, YTArray
-from yt.utilities.fits_image import FITSImageBuffer
+from yt.utilities.fits_image import FITSImageData
import os
@@ -127,7 +127,7 @@
w = subcube.wcs.copy()
w.wcs.crpix[-1] = 0.5
w.wcs.crval[-1] = -0.5*width
- fid = FITSImageBuffer(slab_data, wcs=w)
+ fid = FITSImageData(slab_data, wcs=w)
for hdu in fid:
hdu.header.pop("RESTFREQ", None)
hdu.header.pop("RESTFRQ", None)
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -1,5 +1,5 @@
"""
-FITSImageBuffer Class
+FITSImageData Class
"""
#-----------------------------------------------------------------------------
@@ -16,6 +16,10 @@
from yt.data_objects.construction_data_containers import YTCoveringGridBase
from yt.utilities.on_demand_imports import _astropy, NotAModule
from yt.units.yt_array import YTQuantity, YTArray
+from yt.units import dimensions
+from yt.utilities.parallel_tools.parallel_analysis_interface import \
+ parallel_root_only
+from yt.visualization.volume_rendering.camera import off_axis_projection
import re
pyfits = _astropy.pyfits
@@ -26,18 +30,17 @@
else:
HDUList = pyfits.HDUList
-class FITSImageBuffer(HDUList):
+class FITSImageData(HDUList):
- def __init__(self, data, fields=None, units=None, pixel_scale=None, wcs=None):
- r""" Initialize a FITSImageBuffer object.
+ def __init__(self, data, fields=None, units=None, width=None, wcs=None):
+ r""" Initialize a FITSImageData object.
- FITSImageBuffer contains a list of FITS ImageHDU instances, and
- optionally includes WCS information. It inherits from HDUList, so
- operations such as `writeto` are enabled. Images can be constructed
- from ImageArrays, NumPy arrays, dicts of such arrays,
- FixedResolutionBuffers, and YTCoveringGrids. The latter two are the
- most powerful because WCS information can be constructed from their
- coordinates.
+ FITSImageData contains a collection of FITS ImageHDU instances and
+ WCS information, along with units for each of the images. FITSImageData
+ instances can be constructed from ImageArrays, NumPy arrays, dicts
+ of such arrays, FixedResolutionBuffers, and YTCoveringGrids. The latter
+ two are the most powerful because WCS information can be constructed
+ automatically from their coordinates.
Parameters
----------
@@ -50,10 +53,10 @@
single array one field name must be specified.
units : string
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*.
+ width : float or YTQuantity
+ The width of the image. Either a single value or iterable of values.
+ If a float, assumed to be in *units*. 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.
@@ -66,7 +69,7 @@
>>> 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")
+ >>> f_kpc = FITSImageData(frb, fields="kT", units="kpc")
>>> # This example specifies a specific WCS.
>>> from astropy.wcs import WCS
>>> w = WCS(naxis=self.dimensionality)
@@ -77,46 +80,52 @@
>>> 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 = FITSImageData(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
+ if units is None:
+ units = "cm"
+ if width is None:
+ width = 1.0
- super(FITSImageBuffer, self).__init__()
+ exclude_fields = ['x','y','z','px','py','pz',
+ 'pdx','pdy','pdz','weight_field']
+
+ super(FITSImageData, self).__init__()
if isinstance(fields, string_types):
fields = [fields]
- exclude_fields = ['x', 'y', 'z', 'px', 'py', 'pz',
- 'pdx', 'pdy', 'pdz', 'weight_field']
-
if hasattr(data, 'keys'):
img_data = data
- else:
- img_data = {}
if fields is None:
- mylog.error("Please specify a field name for this array.")
- raise KeyError("Please specify a field name for this array.")
- img_data[fields[0]] = data
+ fields = list(img_data.keys())
+ elif isinstance(data, np.ndarray):
+ if fields is None:
+ mylog.warning("No field name given for this array. Calling it 'image_data'.")
+ fn = 'image_data'
+ fields = [fn]
+ else:
+ fn = fields[0]
+ img_data = {fn: 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("Please specify one or more fields to write.")
+ self.fields = []
+ for fd in fields:
+ if isinstance(fd, tuple):
+ self.fields.append(fd[1])
+ else:
+ self.fields.append(fd)
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))
+ mylog.info("Making a FITS image of field %s" % key)
if first:
hdu = pyfits.PrimaryHDU(np.array(img_data[key]))
first = False
@@ -128,12 +137,8 @@
hdu.header["bunit"] = re.sub('()', '', str(img_data[key].units))
self.append(hdu)
- self.dimensionality = len(self[0].data.shape)
-
- if self.dimensionality == 2:
- self.nx, self.ny = self[0].data.shape
- elif self.dimensionality == 3:
- self.nx, self.ny, self.nz = self[0].data.shape
+ self.shape = self[0].shape
+ self.dimensionality = len(self.shape)
if wcs is None:
w = pywcs.WCS(header=self[0].header, naxis=self.dimensionality)
@@ -141,22 +146,24 @@
# FRBs are a special case where we have coordinate
# information, so we take advantage of this and
# construct the WCS object
- dx = (img_data.bounds[1]-img_data.bounds[0]).in_units(units)/self.nx
- dy = (img_data.bounds[3]-img_data.bounds[2]).in_units(units)/self.ny
- 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)
+ dx = (img_data.bounds[1]-img_data.bounds[0]).in_units(units).v/self.shape[0]
+ dy = (img_data.bounds[3]-img_data.bounds[2]).in_units(units).v/self.shape[1]
+ xctr = 0.5*(img_data.bounds[1]+img_data.bounds[0]).in_units(units).v
+ yctr = 0.5*(img_data.bounds[3]+img_data.bounds[2]).in_units(units).v
center = [xctr, yctr]
cdelt = [dx,dy]
elif isinstance(img_data, YTCoveringGridBase):
cdelt = img_data.dds.in_units(units).v
- center = 0.5*(img_data.left_edge+img_data.right_edge).in_units(units)
+ center = 0.5*(img_data.left_edge+img_data.right_edge).in_units(units).v
else:
# 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
+ # and use the image width to determine the cell widths
+ if not iterable(width):
+ width = [width]*self.dimensionality
+ if isinstance(width[0], YTQuantity):
+ cdelt = [wh.in_units(units).v/n for wh, n in zip(width, self.shape)]
else:
- cdelt = [pixel_scale]*self.dimensionality
+ cdelt = [float(wh)/n for wh, n in zip(width, self.shape)]
center = [0.0]*self.dimensionality
w.wcs.crpix = 0.5*(np.array(self.shape)+1)
w.wcs.crval = center
@@ -178,53 +185,79 @@
for k, v in h.items():
img.header[k] = v
+ def update_header(self, field, key, value):
+ """
+ Update the FITS header for *field* with a
+ *key*, *value* pair. If *field* == "all", all
+ headers will be updated.
+ """
+ if field == "all":
+ for img in self:
+ img.header[key] = value
+ else:
+ if field not in self.keys():
+ raise KeyError("%s not an image!" % field)
+ idx = self.fields.index(field)
+ self[idx].header[key] = value
+
def update_all_headers(self, key, value):
- """
- Update the FITS headers for all images with the
- same *key*, *value* pair.
- """
- for img in self: img.header[key] = value
+ mylog.warning("update_all_headers is deprecated. "+
+ "Use update_header('all', key, value) instead.")
+ self.update_header("all", key, value)
def keys(self):
- return [f.name.lower() for f in self]
+ return self.fields
def has_key(self, key):
- return key in self.keys()
+ return key in self.fields
def values(self):
- return [self[k] for k in self.keys()]
+ return [self[k] for k in self.fields]
def items(self):
- return [(k, self[k]) for k in self.keys()]
+ return [(k, self[k]) for k in self.fields]
- def writeto(self, fileobj, **kwargs):
- pyfits.HDUList(self).writeto(fileobj, **kwargs)
+ @parallel_root_only
+ def writeto(self, fileobj, fields=None, clobber=False, **kwargs):
+ r"""
+ Write all of the fields or a subset of them to a FITS file.
- @property
- def shape(self):
- if self.dimensionality == 2:
- return self.nx, self.ny
- elif self.dimensionality == 3:
- return self.nx, self.ny, self.nz
+ Parameters
+ ----------
+ fileobj : string
+ The name of the file to write to.
+ fields : list of strings, optional
+ The fields to write to the file. If not specified
+ all of the fields in the buffer will be written.
+ clobber : boolean, optional
+ Whether or not to overwrite a previously existing file.
+ Default: False
+ All other keyword arguments are passed to the `writeto`
+ method of `astropy.io.fits.HDUList`.
+ """
+ if fields is None:
+ hdus = pyfits.HDUList(self)
+ else:
+ hdus = pyfits.HDUList()
+ for field in fields:
+ hdus.append(self[field])
+ hdus.writeto(fileobj, clobber=clobber, **kwargs)
def to_glue(self, label="yt", data_collection=None):
"""
- Takes the data in the FITSImageBuffer and exports it to
- Glue (http://www.glueviz.org) for interactive
- analysis. Optionally add a *label*. If you are already within
- the Glue environment, you can pass a *data_collection* object,
- otherwise Glue will be started.
+ Takes the data in the FITSImageData instance and exports it to
+ Glue (http://www.glueviz.org) for interactive analysis. Optionally
+ add a *label*. If you are already within the Glue environment, you
+ can pass a *data_collection* object, otherwise Glue will be started.
"""
from glue.core import DataCollection, Data
from glue.core.coordinates import coordinates_from_header
from glue.qt.glue_application import GlueApplication
- field_dict = dict((key,self[key].data) for key in self.keys())
-
image = Data(label=label)
image.coords = coordinates_from_header(self.wcs.to_header())
- for k,v in field_dict.items():
- image.add_component(v, k)
+ for k,f in self.items():
+ image.add_component(f.data, k)
if data_collection is None:
dc = DataCollection([image])
app = GlueApplication(dc)
@@ -242,64 +275,139 @@
return aplpy.FITSFigure(self, **kwargs)
def get_data(self, field):
+ """
+ Return the data array of the image corresponding to *field*
+ with units attached.
+ """
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
+ if field not in self.keys():
+ raise KeyError("%s not an image!" % field)
+ idx = self.fields.index(field)
+ new_data = YTArray(self[idx].data, self.field_units[field]).in_units(units)
+ self[idx].data = new_data.v
+ self[idx].header["bunit"] = units
self.field_units[field] = units
-axis_wcs = [[1,2],[0,2],[0,1]]
+ def pop(self, key):
+ """
+ Remove a field with name *key*
+ and return it as a new FITSImageData
+ instance.
+ """
+ if key not in self.keys():
+ raise KeyError("%s not an image!" % key)
+ idx = self.fields.index(key)
+ im = super(FITSImageData, self).pop(idx)
+ data = YTArray(im.data, self.field_units[key])
+ self.field_units.pop(key)
+ self.fields.remove(key)
+ return FITSImageData(data, fields=key, wcs=self.wcs)
-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.
+ @classmethod
+ def from_file(cls, filename):
+ """
+ Generate a FITSImageData instance from one previously written to
+ disk.
- 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]
- cunit = ["deg"]*2
- if naxis == 3:
- crval.append(old_wcs.wcs.crval[2])
- cdelt.append(old_wcs.wcs.cdelt[2])
- ctype.append(old_wcs.wcs.ctype[2])
- cunit.append(old_wcs.wcs.cunit[2])
- new_wcs.wcs.crpix = old_wcs.wcs.crpix
- new_wcs.wcs.cdelt = cdelt
- new_wcs.wcs.crval = crval
- new_wcs.wcs.cunit = cunit
- new_wcs.wcs.ctype = ctype
- if crota is not None:
- new_wcs.wcs.crota = crota
- return new_wcs
+ Parameters
+ ----------
+ filename : string
+ The name of the file to open.
+ """
+ f = pyfits.open(filename)
+ data = {}
+ for hdu in f:
+ data[hdu.header["btype"]] = YTArray(hdu.data, hdu.header["bunit"])
+ f.close()
+ return cls(data, wcs=pywcs.WCS(header=hdu.header))
+
+ @classmethod
+ def from_images(cls, image_list):
+ """
+ Generate a new FITSImageData instance from a list of FITSImageData
+ instances.
+
+ Parameters
+ ----------
+ image_list : list of FITSImageData instances
+ The images to be combined.
+ """
+ w = image_list[0].wcs
+ img_shape = image_list[0].shape
+ data = {}
+ for image in image_list:
+ assert_same_wcs(w, image.wcs)
+ if img_shape != image.shape:
+ raise RuntimeError("Images do not have the same shape!")
+ for key in image.keys():
+ data[key] = image.get_data(key)
+ return cls(data, wcs=w)
+
+ def create_sky_wcs(self, sky_center, sky_scale,
+ ctype=["RA---TAN","DEC--TAN"],
+ crota=None, cd=None, pc=None):
+ """
+ Takes a Cartesian WCS and converts it to one in a
+ celestial coordinate system.
+
+ Parameters
+ ----------
+ sky_center : iterable of floats
+ Reference coordinates of the WCS in degrees.
+ sky_scale : tuple or YTQuantity
+ 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 : 2-element ndarray, optional
+ Rotation angles between cartesian coordinates and
+ the celestial coordinates.
+ cd : 2x2-element ndarray, optional
+ Dimensioned coordinate transformation matrix.
+ pc : 2x2-element ndarray, optional
+ Coordinate transformation matrix.
+ """
+ old_wcs = self.wcs
+ naxis = old_wcs.naxis
+ crval = [sky_center[0], sky_center[1]]
+ if isinstance(sky_scale, YTQuantity):
+ scaleq = sky_scale
+ else:
+ scaleq = YTQuantity(sky_scale[0],sky_scale[1])
+ if scaleq.units.dimensions != dimensions.angle/dimensions.length:
+ raise RuntimeError("sky_scale %s not in correct dimensions of angle/length!" % sky_scale)
+ 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]
+ cunit = ["deg"]*2
+ if naxis == 3:
+ crval.append(old_wcs.wcs.crval[2])
+ cdelt.append(old_wcs.wcs.cdelt[2])
+ ctype.append(old_wcs.wcs.ctype[2])
+ cunit.append(old_wcs.wcs.cunit[2])
+ new_wcs.wcs.crpix = old_wcs.wcs.crpix
+ new_wcs.wcs.cdelt = cdelt
+ new_wcs.wcs.crval = crval
+ new_wcs.wcs.cunit = cunit
+ new_wcs.wcs.ctype = ctype
+ if crota is not None:
+ new_wcs.wcs.crota = crota
+ if cd is not None:
+ new_wcs.wcs.cd = cd
+ if pc is not None:
+ new_wcs.wcs.cd = pc
+ self.set_wcs(new_wcs)
+
+class FITSImageBuffer(FITSImageData):
+ pass
def sanitize_fits_unit(unit):
if unit == "Mpc":
@@ -309,11 +417,9 @@
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]]
+axis_wcs = [[1,2],[0,2],[0,1]]
+
+def construct_image(ds, axis, data_source, center, width=None, image_res=None):
if width is None:
width = ds.domain_width[axis_wcs[axis]]
unit = ds.get_smallest_appropriate_unit(width[0])
@@ -323,28 +429,28 @@
width = ds.coordinates.sanitize_width(axis, width, None)
unit = str(width[0].units)
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"))
+ ddims = ds.domain_dimensions*ds.refine_by**ds.index.max_level
+ if iterable(axis):
+ nx = ddims.max()
+ ny = ddims.max()
+ else:
+ nx, ny = [ddims[idx] for idx in axis_wcs[axis]]
else:
if iterable(image_res):
nx, ny = image_res
else:
nx, ny = image_res, image_res
- dx, dy = width[0]/nx, width[1]/ny
+ dx, dy = width[0]/nx, width[1]/ny
crpix = [0.5*(nx+1), 0.5*(ny+1)]
- if hasattr(ds, "wcs"):
+ if hasattr(ds, "wcs") and not iterable(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")[:ds.dimensionality].v
- crval = ds.wcs.wcs_pix2world(ctr_pix.reshape(1,ds.dimensionality))[0]
+ crval = ds.wcs.wcs_pix2world(ctr_pix.reshape(1, ds.dimensionality))[0]
crval = [crval[idx] for idx in axis_wcs[axis]]
else:
- # This is some other kind of dataset
if unit == "unitary":
unit = ds.get_smallest_appropriate_unit(ds.domain_width.max())
elif unit == "code_length":
@@ -353,8 +459,17 @@
cunit = [unit]*2
ctype = ["LINEAR"]*2
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])
+ if iterable(axis):
+ crval = center.in_units(unit)
+ else:
+ crval = [center[idx].in_units(unit) for idx in axis_wcs[axis]]
+ if hasattr(data_source, 'to_frb'):
+ if iterable(axis):
+ frb = data_source.to_frb(width[0], (nx, ny), height=width[1])
+ else:
+ frb = data_source.to_frb(width[0], (nx, ny), center=center, height=width[1])
+ else:
+ frb = None
w = pywcs.WCS(naxis=2)
w.wcs.crpix = crpix
w.wcs.cdelt = cdelt
@@ -363,14 +478,42 @@
w.wcs.ctype = ctype
return w, frb
-class FITSSlice(FITSImageBuffer):
+def assert_same_wcs(wcs1, wcs2):
+ from numpy.testing import assert_allclose
+ assert wcs1.naxis == wcs2.naxis
+ for i in range(wcs1.naxis):
+ assert wcs1.wcs.cunit[i] == wcs2.wcs.cunit[i]
+ assert wcs1.wcs.ctype[i] == wcs2.wcs.ctype[i]
+ assert_allclose(wcs1.wcs.crpix, wcs2.wcs.crpix)
+ assert_allclose(wcs1.wcs.cdelt, wcs2.wcs.cdelt)
+ assert_allclose(wcs1.wcs.crval, wcs2.wcs.crval)
+ crota1 = getattr(wcs1.wcs, "crota", None)
+ crota2 = getattr(wcs2.wcs, "crota", None)
+ if crota1 is None or crota2 is None:
+ assert crota1 == crota2
+ else:
+ assert_allclose(wcs1.wcs.crota, wcs2.wcs.crota)
+ cd1 = getattr(wcs1.wcs, "cd", None)
+ cd2 = getattr(wcs2.wcs, "cd", None)
+ if cd1 is None or cd2 is None:
+ assert cd1 == cd2
+ else:
+ assert_allclose(wcs1.wcs.cd, wcs2.wcs.cd)
+ pc1 = getattr(wcs1.wcs, "pc", None)
+ pc2 = getattr(wcs2.wcs, "pc", None)
+ if pc1 is None or pc2 is None:
+ assert pc1 == pc2
+ else:
+ assert_allclose(wcs1.wcs.pc, wcs2.wcs.pc)
+
+class FITSSlice(FITSImageData):
r"""
- Generate a FITSImageBuffer of an on-axis slice.
+ Generate a FITSImageData of an on-axis slice.
Parameters
----------
- ds : FITSDataset
- The FITS dataset object.
+ ds : :class:`yt.data_objects.api.Dataset`
+ The dataset object.
axis : character or integer
The axis of the slice. One of "x","y","z", or 0,1,2.
fields : string or list of strings
@@ -414,20 +557,18 @@
axis = fix_axis(axis, ds)
center, dcenter = ds.coordinates.sanitize_center(center, axis)
slc = ds.slice(axis, center[axis], **kwargs)
- w, frb = construct_image(slc, center=dcenter, width=width,
- image_res=image_res)
+ w, frb = construct_image(ds, axis, slc, 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)
-class FITSProjection(FITSImageBuffer):
+
+class FITSProjection(FITSImageData):
r"""
- Generate a FITSImageBuffer of an on-axis projection.
+ Generate a FITSImageData of an on-axis projection.
Parameters
----------
- ds : FITSDataset
- The FITS dataset object.
+ ds : :class:`yt.data_objects.api.Dataset`
+ The dataset object.
axis : character or integer
The axis along which to project. One of "x","y","z", or 0,1,2.
fields : string or list of strings
@@ -468,14 +609,161 @@
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", width=None,
+ 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, dcenter = ds.coordinates.sanitize_center(center, axis)
prj = ds.proj(fields[0], axis, weight_field=weight_field, **kwargs)
- w, frb = construct_image(prj, center=dcenter, width=width,
- image_res=image_res)
+ w, frb = construct_image(ds, axis, prj, 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)
+
+class FITSOffAxisSlice(FITSImageData):
+ r"""
+ Generate a FITSImageData of an off-axis slice.
+
+ Parameters
+ ----------
+ ds : :class:`yt.data_objects.api.Dataset`
+ The dataset object.
+ normal : a sequence of floats
+ The vector normal to the projection plane.
+ fields : string or list of strings
+ The fields to slice
+ 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.
+ north_vector : a sequence of floats
+ A vector defining the 'up' direction in the plot. This
+ option sets the orientation of the slicing plane. If not
+ set, an arbitrary grid-aligned north-vector is chosen.
+ """
+ def __init__(self, ds, normal, fields, center='c', width=None, image_res=512,
+ north_vector=None):
+ fields = ensure_list(fields)
+ center, dcenter = ds.coordinates.sanitize_center(center, 4)
+ cut = ds.cutting(normal, center, north_vector=north_vector)
+ center = ds.arr([0.0] * 2, 'code_length')
+ w, frb = construct_image(ds, normal, cut, center, width=width, image_res=image_res)
+ super(FITSOffAxisSlice, self).__init__(frb, fields=fields, wcs=w)
+
+
+class FITSOffAxisProjection(FITSImageData):
+ r"""
+ Generate a FITSImageData of an off-axis projection.
+
+ Parameters
+ ----------
+ ds : :class:`yt.data_objects.api.Dataset`
+ This is the dataset object corresponding to the
+ simulation output to be plotted.
+ normal : a sequence of floats
+ The vector normal to the projection plane.
+ fields : string, list of strings
+ The name of the field(s) to be plotted.
+ 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.
+ depth : A tuple or a float
+ A tuple containing the depth to project through and the string
+ key of the unit: (width, 'unit'). If set to a float, code units
+ are assumed
+ weight_field : string
+ The name of the weighting field. Set to None for no weight.
+ image_res : an int or 2-tuple of ints
+ Specify the resolution of the resulting image.
+ depth_res : an int
+ Specify the resolution of the depth of the projection.
+ north_vector : a sequence of floats
+ A vector defining the 'up' direction in the plot. This
+ option sets the orientation of the slicing plane. If not
+ set, an arbitrary grid-aligned north-vector is chosen.
+ 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.
+ """
+ def __init__(self, ds, normal, fields, center='c', width=(1.0, 'unitary'),
+ weight_field=None, image_res=512, depth_res=256,
+ north_vector=None, depth=(1.0,"unitary"), no_ghost=False, method='integrate'):
+ fields = ensure_list(fields)
+ center, dcenter = ds.coordinates.sanitize_center(center, 4)
+ buf = {}
+ width = ds.coordinates.sanitize_width(normal, width, depth)
+ wd = tuple(el.in_units('code_length').v for el in width)
+ if not iterable(image_res):
+ image_res = (image_res, image_res)
+ res = (image_res[0], image_res[1], depth_res)
+ for field in fields:
+ buf[field] = off_axis_projection(ds, center, normal, wd, res, field,
+ no_ghost=no_ghost, north_vector=north_vector,
+ method=method, weight=weight_field).swapaxes(0, 1)
+ center = ds.arr([0.0] * 2, 'code_length')
+ w, not_an_frb = construct_image(ds, normal, buf, center, width=width, image_res=image_res)
+ super(FITSOffAxisProjection, self).__init__(buf, fields=fields, wcs=w)
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f yt/utilities/tests/test_fits_image.py
--- /dev/null
+++ b/yt/utilities/tests/test_fits_image.py
@@ -0,0 +1,128 @@
+"""
+Unit test FITS image creation in yt.
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# 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.
+#-----------------------------------------------------------------------------
+
+import tempfile
+import os
+import numpy as np
+import shutil
+from yt.testing import fake_random_ds
+from yt.convenience import load
+from numpy.testing import \
+ assert_equal
+from yt.utilities.fits_image import \
+ FITSImageData, FITSProjection, \
+ FITSSlice, FITSOffAxisSlice, \
+ FITSOffAxisProjection, \
+ assert_same_wcs
+from yt.visualization.volume_rendering.camera import \
+ off_axis_projection
+
+def test_fits_image():
+ tmpdir = tempfile.mkdtemp()
+ curdir = os.getcwd()
+ os.chdir(tmpdir)
+
+ fields = ("density", "temperature")
+ units = ('g/cm**3', 'K',)
+ ds = fake_random_ds(64, fields=fields, units=units, nprocs=16,
+ length_unit=100.0)
+
+ prj = ds.proj("density", 2)
+ prj_frb = prj.to_frb((0.5, "unitary"), 128)
+
+ fid1 = FITSImageData(prj_frb, fields=["density","temperature"], units="cm")
+ fits_prj = FITSProjection(ds, "z", ["density","temperature"], image_res=128,
+ width=(0.5,"unitary"))
+
+ yield assert_equal, fid1.get_data("density"), fits_prj.get_data("density")
+ yield assert_equal, fid1.get_data("temperature"), fits_prj.get_data("temperature")
+
+ fid1.writeto("fid1.fits", clobber=True)
+ new_fid1 = FITSImageData.from_file("fid1.fits")
+
+ yield assert_equal, fid1.get_data("density"), new_fid1.get_data("density")
+ yield assert_equal, fid1.get_data("temperature"), new_fid1.get_data("temperature")
+
+ ds2 = load("fid1.fits")
+ ds2.index
+
+ assert ("fits","density") in ds2.field_list
+ assert ("fits","temperature") in ds2.field_list
+
+ dw_cm = ds2.domain_width.in_units("cm")
+
+ assert dw_cm[0].v == 50.
+ assert dw_cm[1].v == 50.
+
+ slc = ds.slice(2, 0.5)
+ slc_frb = slc.to_frb((0.5, "unitary"), 128)
+
+ fid2 = FITSImageData(slc_frb, fields=["density","temperature"], units="cm")
+ fits_slc = FITSSlice(ds, "z", ["density","temperature"], image_res=128,
+ width=(0.5,"unitary"))
+
+ yield assert_equal, fid2.get_data("density"), fits_slc.get_data("density")
+ yield assert_equal, fid2.get_data("temperature"), fits_slc.get_data("temperature")
+
+ dens_img = fid2.pop("density")
+ temp_img = fid2.pop("temperature")
+
+ # This already has some assertions in it, so we don't need to do anything
+ # with it other can just make one
+ fid_comb = FITSImageData.from_images([dens_img, temp_img])
+
+ cut = ds.cutting([0.1, 0.2, -0.9], [0.5, 0.42, 0.6])
+ cut_frb = cut.to_frb((0.5, "unitary"), 128)
+
+ fid3 = FITSImageData(cut_frb, fields=["density","temperature"], units="cm")
+ fits_cut = FITSOffAxisSlice(ds, [0.1, 0.2, -0.9], ["density","temperature"],
+ image_res=128, center=[0.5, 0.42, 0.6],
+ width=(0.5,"unitary"))
+
+ yield assert_equal, fid3.get_data("density"), fits_cut.get_data("density")
+ yield assert_equal, fid3.get_data("temperature"), fits_cut.get_data("temperature")
+
+ fid3.create_sky_wcs([30.,45.], (1.0,"arcsec/kpc"))
+ fid3.writeto("fid3.fits", clobber=True)
+ new_fid3 = FITSImageData.from_file("fid3.fits")
+ assert_same_wcs(fid3.wcs, new_fid3.wcs)
+ assert new_fid3.wcs.wcs.cunit[0] == "deg"
+ assert new_fid3.wcs.wcs.cunit[1] == "deg"
+ assert new_fid3.wcs.wcs.ctype[0] == "RA---TAN"
+ assert new_fid3.wcs.wcs.ctype[1] == "DEC--TAN"
+
+ buf = off_axis_projection(ds, ds.domain_center, [0.1, 0.2, -0.9],
+ 0.5, 128, "density").swapaxes(0, 1)
+ fid4 = FITSImageData(buf, fields="density", width=100.0)
+ fits_oap = FITSOffAxisProjection(ds, [0.1, 0.2, -0.9], "density",
+ width=(0.5,"unitary"), image_res=128,
+ depth_res=128, depth=(0.5,"unitary"))
+
+ yield assert_equal, fid4.get_data("density"), fits_oap.get_data("density")
+
+ cvg = ds.covering_grid(ds.index.max_level, [0.25,0.25,0.25],
+ [32, 32, 32], fields=["density","temperature"])
+ fid5 = FITSImageData(cvg, fields=["density","temperature"])
+ assert fid5.dimensionality == 3
+
+ fid5.update_header("density", "time", 0.1)
+ fid5.update_header("all", "units", "cgs")
+
+ assert fid5["density"].header["time"] == 0.1
+ assert fid5["temperature"].header["units"] == "cgs"
+ assert fid5["density"].header["units"] == "cgs"
+
+ os.chdir(curdir)
+ shutil.rmtree(tmpdir)
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f yt/visualization/fixed_resolution.py
--- a/yt/visualization/fixed_resolution.py
+++ b/yt/visualization/fixed_resolution.py
@@ -327,11 +327,11 @@
the length units that the coordinates are written in, default 'cm'.
"""
- from yt.utilities.fits_image import FITSImageBuffer
+ from yt.utilities.fits_image import FITSImageData
if fields is None: fields = list(self.data.keys())
- fib = FITSImageBuffer(self, fields=fields, units=units)
+ fib = FITSImageData(self, fields=fields, units=units)
if other_keys is not None:
for k,v in other_keys.items():
fib.update_all_headers(k,v)
@@ -470,7 +470,7 @@
def __getitem__(self, item):
if item in self.data: return self.data[item]
- mylog.info("Making a fixed resolutuion buffer of (%s) %d by %d" % \
+ mylog.info("Making a fixed resolution buffer of (%s) %d by %d" % \
(item, self.buff_size[0], self.buff_size[1]))
dd = self.data_source
width = self.ds.arr((self.bounds[1] - self.bounds[0],
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py
+++ b/yt/visualization/volume_rendering/camera.py
@@ -2252,7 +2252,6 @@
def temp_weightfield(a, b):
tr = b[f].astype("float64") * b[w]
return b.apply_units(tr, a.units)
- return tr
return temp_weightfield
ds.field_info.add_field(weightfield,
function=_make_wf(field, weight))
diff -r 73ca90e953997cb93a5b04ff33090839f55e2feb -r e4bc6b4f239bd0a557e3819586330a75028f4c4f yt/visualization/volume_rendering/image_handling.py
--- a/yt/visualization/volume_rendering/image_handling.py
+++ b/yt/visualization/volume_rendering/image_handling.py
@@ -33,14 +33,14 @@
f.create_dataset("A", data=image[:,:,3])
f.close()
if fits:
- from yt.utilities.fits_image import FITSImageBuffer
+ from yt.utilities.fits_image import FITSImageData
data = {}
data["r"] = image[:,:,0]
data["g"] = image[:,:,1]
data["b"] = image[:,:,2]
data["a"] = image[:,:,3]
nx, ny = data["r"].shape
- fib = FITSImageBuffer(data)
+ fib = FITSImageData(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