[yt-svn] commit/yt: 28 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Jun 18 06:17:37 PDT 2015
28 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/7d409e929b96/
Changeset: 7d409e929b96
Branch: yt
User: jzuhone
Date: 2015-05-15 13:38:06+00:00
Summary: Docstring should reflect the correct module name
Affected #: 1 file
diff -r 4b4b61882407d5a5c909647b24a0147f67eee5ac -r 7d409e929b9677b2fd99232163d6941911468a82 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.
"""
#-----------------------------------------------------------------------------
https://bitbucket.org/yt_analysis/yt/commits/8378ed9f4cdf/
Changeset: 8378ed9f4cdf
Branch: yt
User: jzuhone
Date: 2015-05-15 13:38:37+00:00
Summary: Remove second return statement
Affected #: 1 file
diff -r 7d409e929b9677b2fd99232163d6941911468a82 -r 8378ed9f4cdff934c728b57c2da266f0b2b696d2 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))
https://bitbucket.org/yt_analysis/yt/commits/dbbc3defc965/
Changeset: dbbc3defc965
Branch: yt
User: jzuhone
Date: 2015-05-15 13:51:29+00:00
Summary: FITSImageBuffer --> FITSImageData. Fixing some small issues.
Affected #: 1 file
diff -r 8378ed9f4cdff934c728b57c2da266f0b2b696d2 -r dbbc3defc9655d675e02487a50a2a697755f1d6a 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
"""
#-----------------------------------------------------------------------------
@@ -66,7 +66,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)
@@ -141,15 +141,15 @@
# 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
@@ -209,22 +209,19 @@
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,v in self.items():
+ image.add_component(v.v, k)
if data_collection is None:
dc = DataCollection([image])
app = GlueApplication(dc)
@@ -365,12 +362,12 @@
class FITSSlice(FITSImageBuffer):
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
https://bitbucket.org/yt_analysis/yt/commits/f165aa37acd5/
Changeset: f165aa37acd5
Branch: yt
User: jzuhone
Date: 2015-05-15 14:13:38+00:00
Summary: We won't subclass from HDUList anymore, so that we can implement __getitem__ to give back YTArrays. Instead, we will use HDUList under the hood and replicate some of its methods
Affected #: 1 file
diff -r dbbc3defc9655d675e02487a50a2a697755f1d6a -r f165aa37acd5f0b3e56f3b9aa2b9cfa6ce94d81b yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -126,17 +126,13 @@
hdu.header["btype"] = key
if hasattr(img_data[key], "units"):
hdu.header["bunit"] = re.sub('()', '', str(img_data[key].units))
- self.append(hdu)
+ self.hdulist.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.hdulist[0].shape
+ self.dimensionality = len(self.shape)
if wcs is None:
- w = pywcs.WCS(header=self[0].header, naxis=self.dimensionality)
+ w = pywcs.WCS(header=self.hdulist[0].header, naxis=self.dimensionality)
if isinstance(img_data, FixedResolutionBuffer):
# FRBs are a special case where we have coordinate
# information, so we take advantage of this and
@@ -174,7 +170,7 @@
"""
self.wcs = wcs
h = self.wcs.to_header()
- for img in self:
+ for img in self.hdulist:
for k, v in h.items():
img.header[k] = v
@@ -186,26 +182,68 @@
for img in self: img.header[key] = value
def keys(self):
- return [f.name.lower() for f in self]
+ return self.fields
+
+ def __getitem__(self, field):
+ if field not in self.keys():
+ raise KeyError("%s not an image!" % field)
+ idx = self.fields.index(field)
+ return YTArray(self.hdulist[idx].data, self.field_units[field])
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)
+ def get_header(self, field):
+ """
+ Get the FITS header for a specific field.
- @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
+ ----------
+ field : string
+ The field for which to get the corresponding header.
+ """
+ if field not in self.keys():
+ raise KeyError("%s not an image!" % field)
+ idx = self.fields.index(field)
+ return self.hdulist[idx].header
+
+ @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.
+
+ 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 = self.hdulist
+ else:
+ hdus = pyfits.HDUList()
+ for field in fields:
+ hdus.append(self.hdulist[field])
+ hdus.writeto(fileobj, clobber=clobber, **kwargs)
+
+ def info(self):
+ """
+ Display information about the underlying FITS file.
+ """
+ self.hdulist.info()
def to_glue(self, label="yt", data_collection=None):
"""
@@ -236,18 +274,18 @@
`aplpy.FITSFigure` constructor.
"""
import aplpy
- return aplpy.FITSFigure(self, **kwargs)
-
- def get_data(self, field):
- return YTArray(self[field].data, self.field_units[field])
+ return aplpy.FITSFigure(self.hdulist, **kwargs)
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)
+ new_data = self[field].in_units(units)
+ idx = self.fields.index(field)
+ self.hdulist[idx].data = new_data.v
+ self.hdulist[idx].header["bunit"] = units
self.field_units[field] = units
axis_wcs = [[1,2],[0,2],[0,1]]
https://bitbucket.org/yt_analysis/yt/commits/4b6df9381fc4/
Changeset: 4b6df9381fc4
Branch: yt
User: jzuhone
Date: 2015-05-15 14:19:00+00:00
Summary: Specify a width for the entire FITS image instead of the pixel width, and figure out the latter that way.
Affected #: 1 file
diff -r f165aa37acd5f0b3e56f3b9aa2b9cfa6ce94d81b -r 4b6df9381fc45e4ba5fe63e1fe5e1d495c535b3d yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -50,10 +50,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.
@@ -148,11 +148,13 @@
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 = [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
https://bitbucket.org/yt_analysis/yt/commits/916129ccf00b/
Changeset: 916129ccf00b
Branch: yt
User: jzuhone
Date: 2015-05-15 14:38:17+00:00
Summary: FITSImageBuffer --> FITSImageData and cleaning up some other things
Affected #: 1 file
diff -r 4b6df9381fc45e4ba5fe63e1fe5e1d495c535b3d -r 916129ccf00bed605526e2a1595884ff2e12ba45 yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -21,23 +21,17 @@
pyfits = _astropy.pyfits
pywcs = _astropy.pywcs
-if isinstance(pyfits, NotAModule):
- HDUList = object
-else:
- HDUList = pyfits.HDUList
+class FITSImageData(object):
-class FITSImageBuffer(HDUList):
+ def __init__(self, data, fields=None, units=None, width=None, wcs=None):
+ r""" Initialize a FITSImageData object.
- 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
- 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
----------
@@ -77,46 +71,47 @@
>>> 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']
- if isinstance(fields, string_types):
+ self.hdulist = pyfits.HDUList()
+
+ 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 = fields
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
https://bitbucket.org/yt_analysis/yt/commits/f49b7db74d04/
Changeset: f49b7db74d04
Branch: yt
User: jzuhone
Date: 2015-05-15 14:39:50+00:00
Summary: Use one method to update headers, either just one or all of them at once
Affected #: 1 file
diff -r 916129ccf00bed605526e2a1595884ff2e12ba45 -r f49b7db74d04fc570b363d90b96d953d8389f45f yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -171,12 +171,20 @@
for k, v in h.items():
img.header[k] = v
- def update_all_headers(self, key, value):
+ def update_header(self, field, key, value):
"""
- Update the FITS headers for all images with the
- same *key*, *value* pair.
+ Update the FITS header for *field* with a
+ *key*, *value* pair. If *field* == "all", all
+ headers will be updated.
"""
- for img in self: img.header[key] = value
+ if field == "all":
+ for img in self.hdulist:
+ img.header[key] = value
+ else:
+ if field not in self.keys():
+ raise KeyError("%s not an image!" % field)
+ idx = self.fields.index(field)
+ self.hdulist[idx].header[key] = value
def keys(self):
return self.fields
https://bitbucket.org/yt_analysis/yt/commits/bfd2f0e2f7e9/
Changeset: bfd2f0e2f7e9
Branch: yt
User: jzuhone
Date: 2015-05-15 14:41:16+00:00
Summary: Pop a FITS image from a FITSImageData instance and create a new instance from it
Affected #: 1 file
diff -r f49b7db74d04fc570b363d90b96d953d8389f45f -r bfd2f0e2f7e9718961948bf5313cde707b81dc34 yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -293,7 +293,20 @@
self.hdulist[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)
+ data = self[key]
+ idx = self.fields.index(key)
+ self.hdulist.pop(idx)
+ 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):
https://bitbucket.org/yt_analysis/yt/commits/059a6449789d/
Changeset: 059a6449789d
Branch: yt
User: jzuhone
Date: 2015-05-15 14:43:49+00:00
Summary: Add classmethods to create new FITSImageData from files or from lists of FITSImageData. create_sky_wcs should be a method of FITSImageData.
Affected #: 1 file
diff -r bfd2f0e2f7e9718961948bf5313cde707b81dc34 -r 059a6449789dbf18dbdc534f148e3e6fc5e04e54 yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -308,51 +308,103 @@
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 k,v in image.items():
+ data[k] = v
+ 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)
def sanitize_fits_unit(unit):
if unit == "Mpc":
https://bitbucket.org/yt_analysis/yt/commits/dd335f94afe6/
Changeset: dd335f94afe6
Branch: yt
User: jzuhone
Date: 2015-05-15 14:44:49+00:00
Summary: Adding FITSOffAxisSlice and FITSOffAxisProjection classes. Making some other necessary changes.
Affected #: 1 file
diff -r 059a6449789dbf18dbdc534f148e3e6fc5e04e54 -r dd335f94afe6d2b0bfed215629983aa9867dc0a2 yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -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
@@ -414,11 +418,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])
@@ -428,28 +430,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*2**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":
@@ -458,8 +460,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
@@ -468,7 +479,35 @@
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 FITSImageData of an on-axis slice.
@@ -519,20 +558,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
@@ -579,8 +616,155 @@
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)
https://bitbucket.org/yt_analysis/yt/commits/a213dfc2cff8/
Changeset: a213dfc2cff8
Branch: yt
User: jzuhone
Date: 2015-05-15 14:45:31+00:00
Summary: Catch sunyaev_zeldovich up with the new FITSImageData changes.
Affected #: 1 file
diff -r dd335f94afe6d2b0bfed215629983aa9867dc0a2 -r a213dfc2cff8d3dced45dba2b2bc31073fb052b8 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
https://bitbucket.org/yt_analysis/yt/commits/47efdc422932/
Changeset: 47efdc422932
Branch: yt
User: jzuhone
Date: 2015-05-15 14:46:01+00:00
Summary: Implement the new FITSImageData changes in ppv_cube, and some small cosmetic changes to the module
Affected #: 1 file
diff -r a213dfc2cff8d3dced45dba2b2bc31073fb052b8 -r 47efdc422932ca439f58244bcef87ceebadbffb2 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
@@ -89,6 +89,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 +154,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],
@@ -214,16 +214,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 +249,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 +271,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 +295,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 +310,13 @@
def __getitem__(self, item):
return self.data[item]
+
+ 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
https://bitbucket.org/yt_analysis/yt/commits/082397fb64d1/
Changeset: 082397fb64d1
Branch: yt
User: jzuhone
Date: 2015-05-15 14:46:20+00:00
Summary: Catching up fixed_resolution with the FITSImageData changes
Affected #: 1 file
diff -r 47efdc422932ca439f58244bcef87ceebadbffb2 -r 082397fb64d193f35f0bf485bed9bc46e15cc769 yt/visualization/fixed_resolution.py
--- a/yt/visualization/fixed_resolution.py
+++ b/yt/visualization/fixed_resolution.py
@@ -318,14 +318,14 @@
requested.
"""
- from yt.utilities.fits_image import FITSImageBuffer
+ from yt.utilities.fits_image import FITSImageData
extra_fields = ['x','y','z','px','py','pz','pdx','pdy','pdz','weight_field']
if fields is None:
fields = [field[-1] for field in self.data_source.field_data
if field not in extra_fields]
- 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)
@@ -410,7 +410,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],
https://bitbucket.org/yt_analysis/yt/commits/a57cf2c79d41/
Changeset: a57cf2c79d41
Branch: yt
User: jzuhone
Date: 2015-05-15 14:46:37+00:00
Summary: Catching up image_handling with the FITSImageData changes
Affected #: 1 file
diff -r 082397fb64d193f35f0bf485bed9bc46e15cc769 -r a57cf2c79d41f1385c17f5157e323de37824f7a9 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):
https://bitbucket.org/yt_analysis/yt/commits/51fe799a7744/
Changeset: 51fe799a7744
Branch: yt
User: jzuhone
Date: 2015-05-15 14:46:52+00:00
Summary: FITSImageData tests
Affected #: 1 file
diff -r a57cf2c79d41f1385c17f5157e323de37824f7a9 -r 51fe799a7744a407967da2f116635437203a4b80 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["density"], fits_prj["density"]
+ yield assert_equal, fid1["temperature"], fits_prj["temperature"]
+
+ fid1.writeto("fid1.fits", clobber=True)
+ new_fid1 = FITSImageData.from_file("fid1.fits")
+
+ yield assert_equal, fid1["density"], new_fid1["density"]
+ yield assert_equal, fid1["temperature"], new_fid1["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["density"], fits_slc["density"]
+ yield assert_equal, fid2["temperature"], fits_slc["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["density"], fits_cut["density"]
+ yield assert_equal, fid3["temperature"], fits_cut["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["density"], fits_oap["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.get_header("density")["time"] == 0.1
+ assert fid5.get_header("temperature")["units"] == "cgs"
+ assert fid5.get_header("density")["units"] == "cgs"
+
+ os.chdir(curdir)
+ shutil.rmtree(tmpdir)
https://bitbucket.org/yt_analysis/yt/commits/791fa17f0474/
Changeset: 791fa17f0474
Branch: yt
User: jzuhone
Date: 2015-05-29 02:30:05+00:00
Summary: Updating docs for FITS images.
Affected #: 3 files
diff -r 51fe799a7744a407967da2f116635437203a4b80 -r 791fa17f047401a9c4db12bb60c61a7388de109f 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 51fe799a7744a407967da2f116635437203a4b80 -r 791fa17f047401a9c4db12bb60c61a7388de109f doc/source/visualizing/FITSImageData.ipynb
--- /dev/null
+++ b/doc/source/visualizing/FITSImageData.ipynb
@@ -0,0 +1,392 @@
+{
+ "metadata": {
+ "name": "",
+ "signature": "sha256:904e7ea07d0e4df6a73f57089fc9cd8347b056f0e763d6a8c4782182d1b6e2bb"
+ },
+ "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": [
+ " `prj_fits` implements some methods from the [AstroPy `HDUList`](http://astropy.readthedocs.org/en/latest/io/fits/api/hdulists.html) class. `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, using `get_header`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits.get_header(\"temperature\")"
+ ],
+ "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. 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",
+ "print prj_fits[\"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 now:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "prj_fits.get_header(\"temperature\")"
+ ],
+ "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.get_header(\"density\")[\"time\"]\n",
+ "print fid_frb.get_header(\"temperature\")[\"scale\"]"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ }
+ ],
+ "metadata": {}
+ }
+ ]
+}
\ No newline at end of file
diff -r 51fe799a7744a407967da2f116635437203a4b80 -r 791fa17f047401a9c4db12bb60c61a7388de109f 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
https://bitbucket.org/yt_analysis/yt/commits/0e1aba8962f2/
Changeset: 0e1aba8962f2
Branch: yt
User: jzuhone
Date: 2015-05-29 05:07:19+00:00
Summary: Forgot this in here
Affected #: 1 file
diff -r 791fa17f047401a9c4db12bb60c61a7388de109f -r 0e1aba8962f244487268becafc6ee53691f70330 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
https://bitbucket.org/yt_analysis/yt/commits/19aefe028c92/
Changeset: 19aefe028c92
Branch: yt
User: jzuhone
Date: 2015-05-29 05:08:47+00:00
Summary: Decided to revert back to some original behaviors. Inheriting directly from HDUList. No longer overloading __getitem__.
Affected #: 1 file
diff -r 0e1aba8962f244487268becafc6ee53691f70330 -r 19aefe028c92212445658c43adab9cec771a4001 yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -25,7 +25,12 @@
pyfits = _astropy.pyfits
pywcs = _astropy.pywcs
-class FITSImageData(object):
+if isinstance(pyfits, NotAModule):
+ HDUList = object
+else:
+ HDUList = pyfits.HDUList
+
+class FITSImageData(HDUList):
def __init__(self, data, fields=None, units=None, width=None, wcs=None):
r""" Initialize a FITSImageData object.
@@ -79,17 +84,17 @@
>>> f_deg.writeto("temp.fits")
"""
- if units is None:
+ if units is None:
units = "cm"
- if width is None:
+ if width is None:
width = 1.0
exclude_fields = ['x','y','z','px','py','pz',
'pdx','pdy','pdz','weight_field']
- self.hdulist = pyfits.HDUList()
+ super(FITSImageData, self).__init__()
- if isinstance(fields, string_types):
+ if isinstance(fields, string_types):
fields = [fields]
if hasattr(data, 'keys'):
@@ -125,13 +130,13 @@
hdu.header["btype"] = key
if hasattr(img_data[key], "units"):
hdu.header["bunit"] = re.sub('()', '', str(img_data[key].units))
- self.hdulist.append(hdu)
+ self.append(hdu)
- self.shape = self.hdulist[0].shape
+ self.shape = self[0].shape
self.dimensionality = len(self.shape)
if wcs is None:
- w = pywcs.WCS(header=self.hdulist[0].header, naxis=self.dimensionality)
+ w = pywcs.WCS(header=self[0].header, naxis=self.dimensionality)
if isinstance(img_data, FixedResolutionBuffer):
# FRBs are a special case where we have coordinate
# information, so we take advantage of this and
@@ -171,7 +176,7 @@
"""
self.wcs = wcs
h = self.wcs.to_header()
- for img in self.hdulist:
+ for img in self:
for k, v in h.items():
img.header[k] = v
@@ -182,23 +187,22 @@
headers will be updated.
"""
if field == "all":
- for img in self.hdulist:
+ 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.hdulist[idx].header[key] = value
+ self[idx].header[key] = value
+
+ def update_all_headers(self, 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 self.fields
- def __getitem__(self, field):
- if field not in self.keys():
- raise KeyError("%s not an image!" % field)
- idx = self.fields.index(field)
- return YTArray(self.hdulist[idx].data, self.field_units[field])
-
def has_key(self, key):
return key in self.fields
@@ -208,20 +212,6 @@
def items(self):
return [(k, self[k]) for k in self.fields]
- def get_header(self, field):
- """
- Get the FITS header for a specific field.
-
- Parameters
- ----------
- field : string
- The field for which to get the corresponding header.
- """
- if field not in self.keys():
- raise KeyError("%s not an image!" % field)
- idx = self.fields.index(field)
- return self.hdulist[idx].header
-
@parallel_root_only
def writeto(self, fileobj, fields=None, clobber=False, **kwargs):
r"""
@@ -241,19 +231,13 @@
method of `astropy.io.fits.HDUList`.
"""
if fields is None:
- hdus = self.hdulist
+ hdus = pyfits.HDUList(self)
else:
hdus = pyfits.HDUList()
for field in fields:
- hdus.append(self.hdulist[field])
+ hdus.append(self[field])
hdus.writeto(fileobj, clobber=clobber, **kwargs)
- def info(self):
- """
- Display information about the underlying FITS file.
- """
- self.hdulist.info()
-
def to_glue(self, label="yt", data_collection=None):
"""
Takes the data in the FITSImageData instance and exports it to
@@ -267,8 +251,8 @@
image = Data(label=label)
image.coords = coordinates_from_header(self.wcs.to_header())
- for k,v in self.items():
- image.add_component(v.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)
@@ -283,7 +267,14 @@
`aplpy.FITSFigure` constructor.
"""
import aplpy
- return aplpy.FITSFigure(self.hdulist, **kwargs)
+ 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):
"""
@@ -291,10 +282,10 @@
"""
if field not in self.keys():
raise KeyError("%s not an image!" % field)
- new_data = self[field].in_units(units)
idx = self.fields.index(field)
- self.hdulist[idx].data = new_data.v
- self.hdulist[idx].header["bunit"] = units
+ 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
def pop(self, key):
@@ -305,9 +296,8 @@
"""
if key not in self.keys():
raise KeyError("%s not an image!" % key)
- data = self[key]
idx = self.fields.index(key)
- self.hdulist.pop(idx)
+ data = YTArray(super(FITSImageData, self).pop(idx), self.field_units[key])
self.field_units.pop(key)
self.fields.remove(key)
return FITSImageData(data, fields=key, wcs=self.wcs)
@@ -353,10 +343,10 @@
return cls(data, wcs=w)
def create_sky_wcs(self, sky_center, sky_scale,
- ctype=["RA---TAN","DEC--TAN"],
+ ctype=["RA---TAN","DEC--TAN"],
crota=None, cd=None, pc=None):
"""
- Takes a Cartesian WCS and converts it to one in a
+ Takes a Cartesian WCS and converts it to one in a
celestial coordinate system.
Parameters
@@ -410,6 +400,9 @@
new_wcs.wcs.cd = pc
self.set_wcs(new_wcs)
+class FITSImageBuffer(FITSImageData):
+ pass
+
def sanitize_fits_unit(unit):
if unit == "Mpc":
mylog.info("Changing FITS file unit to kpc.")
@@ -610,7 +603,7 @@
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)
@@ -668,7 +661,7 @@
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,
+ 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)
@@ -750,8 +743,8 @@
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,
+ 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)
@@ -762,8 +755,8 @@
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,
+ 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)
https://bitbucket.org/yt_analysis/yt/commits/04b4e2efd6da/
Changeset: 04b4e2efd6da
Branch: yt
User: jzuhone
Date: 2015-05-29 05:26:53+00:00
Summary: Fixed some bugs
Affected #: 1 file
diff -r 19aefe028c92212445658c43adab9cec771a4001 -r 04b4e2efd6da316bb24d2a46ade57e781822f55e yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -297,7 +297,8 @@
if key not in self.keys():
raise KeyError("%s not an image!" % key)
idx = self.fields.index(key)
- data = YTArray(super(FITSImageData, self).pop(idx), self.field_units[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)
@@ -338,8 +339,8 @@
assert_same_wcs(w, image.wcs)
if img_shape != image.shape:
raise RuntimeError("Images do not have the same shape!")
- for k,v in image.items():
- data[k] = v
+ 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,
https://bitbucket.org/yt_analysis/yt/commits/6f7e872d6305/
Changeset: 6f7e872d6305
Branch: yt
User: jzuhone
Date: 2015-05-29 05:31:39+00:00
Summary: Fixing tests
Affected #: 1 file
diff -r 04b4e2efd6da316bb24d2a46ade57e781822f55e -r 6f7e872d6305861d67ea0d7032e98825dc5159a2 yt/utilities/tests/test_fits_image.py
--- a/yt/utilities/tests/test_fits_image.py
+++ b/yt/utilities/tests/test_fits_image.py
@@ -46,14 +46,14 @@
fits_prj = FITSProjection(ds, "z", ["density","temperature"], image_res=128,
width=(0.5,"unitary"))
- yield assert_equal, fid1["density"], fits_prj["density"]
- yield assert_equal, fid1["temperature"], fits_prj["temperature"]
+ 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["density"], new_fid1["density"]
- yield assert_equal, fid1["temperature"], new_fid1["temperature"]
+ 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
@@ -73,8 +73,8 @@
fits_slc = FITSSlice(ds, "z", ["density","temperature"], image_res=128,
width=(0.5,"unitary"))
- yield assert_equal, fid2["density"], fits_slc["density"]
- yield assert_equal, fid2["temperature"], fits_slc["temperature"]
+ 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")
@@ -91,8 +91,8 @@
image_res=128, center=[0.5, 0.42, 0.6],
width=(0.5,"unitary"))
- yield assert_equal, fid3["density"], fits_cut["density"]
- yield assert_equal, fid3["temperature"], fits_cut["temperature"]
+ 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)
@@ -110,7 +110,7 @@
width=(0.5,"unitary"), image_res=128,
depth_res=128, depth=(0.5,"unitary"))
- yield assert_equal, fid4["density"], fits_oap["density"]
+ 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"])
@@ -120,9 +120,9 @@
fid5.update_header("density", "time", 0.1)
fid5.update_header("all", "units", "cgs")
- assert fid5.get_header("density")["time"] == 0.1
- assert fid5.get_header("temperature")["units"] == "cgs"
- assert fid5.get_header("density")["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)
https://bitbucket.org/yt_analysis/yt/commits/e05cffc36d98/
Changeset: e05cffc36d98
Branch: yt
User: jzuhone
Date: 2015-05-29 05:34:19+00:00
Summary: Updated docs
Affected #: 1 file
diff -r 6f7e872d6305861d67ea0d7032e98825dc5159a2 -r e05cffc36d9812ba53fbfbc3bb5445929f1a3587 doc/source/visualizing/FITSImageData.ipynb
--- a/doc/source/visualizing/FITSImageData.ipynb
+++ b/doc/source/visualizing/FITSImageData.ipynb
@@ -1,7 +1,7 @@
{
"metadata": {
"name": "",
- "signature": "sha256:904e7ea07d0e4df6a73f57089fc9cd8347b056f0e763d6a8c4782182d1b6e2bb"
+ "signature": "sha256:c7de5ef190feaa2289595aec7eaa05db02fd535e408e0d04aa54088b0bd3ebae"
},
"nbformat": 3,
"nbformat_minor": 0,
@@ -92,7 +92,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- " `prj_fits` implements some methods from the [AstroPy `HDUList`](http://astropy.readthedocs.org/en/latest/io/fits/api/hdulists.html) class. `info` shows us the contents of the virtual FITS file:"
+ "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:"
]
},
{
@@ -109,14 +109,14 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "We can also look at the header for a particular field, using `get_header`:"
+ "We can also look at the header for a particular field:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
- "prj_fits.get_header(\"temperature\")"
+ "prj_fits[\"temperature\"].header"
],
"language": "python",
"metadata": {},
@@ -126,7 +126,24 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "where we can see that the temperature units are in Kelvin and the cell widths are in kiloparsecs. We can use the `set_unit` method to change the units of a particular field:"
+ "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:"
]
},
{
@@ -134,7 +151,7 @@
"collapsed": false,
"input": [
"prj_fits.set_unit(\"temperature\",\"R\")\n",
- "print prj_fits[\"temperature\"]"
+ "prj_fits.get_data(\"temperature\")"
],
"language": "python",
"metadata": {},
@@ -343,14 +360,14 @@
"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 now:"
+ "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.get_header(\"temperature\")"
+ "prj_fits[\"temperature\"].header"
],
"language": "python",
"metadata": {},
@@ -378,8 +395,8 @@
"cell_type": "code",
"collapsed": false,
"input": [
- "print fid_frb.get_header(\"density\")[\"time\"]\n",
- "print fid_frb.get_header(\"temperature\")[\"scale\"]"
+ "print fid_frb[\"density\"].header[\"time\"]\n",
+ "print fid_frb[\"temperature\"].header[\"scale\"]"
],
"language": "python",
"metadata": {},
https://bitbucket.org/yt_analysis/yt/commits/fbbc537ec75c/
Changeset: fbbc537ec75c
Branch: yt
User: jzuhone
Date: 2015-05-29 05:40:04+00:00
Summary: Sanity check against tuple field keys
Affected #: 1 file
diff -r e05cffc36d9812ba53fbfbc3bb5445929f1a3587 -r fbbc537ec75c619b0351f339231f964d7e80262b yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -110,7 +110,12 @@
fn = fields[0]
img_data = {fn: data}
- self.fields = fields
+ 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 = {}
https://bitbucket.org/yt_analysis/yt/commits/696a6134f3f0/
Changeset: 696a6134f3f0
Branch: yt
User: jzuhone
Date: 2015-05-29 13:07:45+00:00
Summary: Bugfix
Affected #: 1 file
diff -r fbbc537ec75c619b0351f339231f964d7e80262b -r 696a6134f3f070d4b55fcfc4c2554db563aab11e 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
@@ -171,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:
https://bitbucket.org/yt_analysis/yt/commits/16985637841f/
Changeset: 16985637841f
Branch: yt
User: jzuhone
Date: 2015-05-29 19:45:01+00:00
Summary: Merge
Affected #: 29 files
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 .hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -13,6 +13,7 @@
yt/frontends/ramses/_ramses_reader.cpp
yt/geometry/fake_octree.c
yt/geometry/grid_container.c
+yt/geometry/grid_visitors.c
yt/geometry/oct_container.c
yt/geometry/oct_visitors.c
yt/geometry/particle_deposit.c
@@ -25,6 +26,7 @@
yt/utilities/spatial/ckdtree.c
yt/utilities/lib/alt_ray_tracers.c
yt/utilities/lib/amr_kdtools.c
+yt/utilities/lib/bitarray.c
yt/utilities/lib/CICDeposit.c
yt/utilities/lib/ContourFinding.c
yt/utilities/lib/DepthFirstOctree.c
@@ -39,6 +41,7 @@
yt/utilities/lib/misc_utilities.c
yt/utilities/lib/Octree.c
yt/utilities/lib/origami.c
+yt/utilities/lib/pixelization_routines.c
yt/utilities/lib/png_writer.c
yt/utilities/lib/PointsInVolume.c
yt/utilities/lib/QuadTree.c
@@ -59,3 +62,4 @@
doc/source/reference/api/generated/*
doc/_temp/*
doc/source/bootcamp/.ipynb_checkpoints/
+dist
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 .python-version
--- /dev/null
+++ b/.python-version
@@ -0,0 +1,1 @@
+2.7.9
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 README
--- a/README
+++ b/README
@@ -20,4 +20,4 @@
For more information on installation, what to do if you run into problems, or
ways to help development, please visit our website.
-Enjoy!
+Enjoy!
\ No newline at end of file
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 doc/source/analyzing/analysis_modules/halo_finders.rst
--- a/doc/source/analyzing/analysis_modules/halo_finders.rst
+++ b/doc/source/analyzing/analysis_modules/halo_finders.rst
@@ -116,7 +116,7 @@
the width of the smallest grid element in the simulation from the
last data snapshot (i.e. the one where time has evolved the
longest) in the time series:
- ``ds_last.index.get_smallest_dx() * ds_last['mpch']``.
+ ``ds_last.index.get_smallest_dx() * ds_last['Mpch']``.
* ``total_particles``, if supplied, this is a pre-calculated
total number of dark matter
particles present in the simulation. For example, this is useful
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 doc/source/analyzing/time_series_analysis.rst
--- a/doc/source/analyzing/time_series_analysis.rst
+++ b/doc/source/analyzing/time_series_analysis.rst
@@ -79,9 +79,7 @@
Analyzing an Entire Simulation
------------------------------
-.. note:: Currently only implemented for Enzo. Other simulation types coming
- soon. Until then, rely on the above prescription for creating
- ``DatasetSeries`` objects.
+.. note:: Implemented for: Enzo, Gadget, OWLS.
The parameter file used to run a simulation contains all the information
necessary to know what datasets should be available. The ``simulation``
@@ -93,8 +91,7 @@
.. code-block:: python
import yt
- my_sim = yt.simulation('enzo_tiny_cosmology/32Mpc_32.enzo', 'Enzo',
- find_outputs=False)
+ my_sim = yt.simulation('enzo_tiny_cosmology/32Mpc_32.enzo', 'Enzo')
Then, create a ``DatasetSeries`` object with the
:meth:`frontends.enzo.simulation_handling.EnzoSimulation.get_time_series`
@@ -123,10 +120,10 @@
to select a subset of the total data:
* ``time_data`` (*bool*): Whether or not to include time outputs when
- gathering datasets for time series. Default: True.
+ gathering datasets for time series. Default: True. (Enzo only)
* ``redshift_data`` (*bool*): Whether or not to include redshift outputs
- when gathering datasets for time series. Default: True.
+ when gathering datasets for time series. Default: True. (Enzo only)
* ``initial_time`` (*float*): The earliest time for outputs to be included.
If None, the initial time of the simulation is used. This can be used in
@@ -139,15 +136,12 @@
* ``times`` (*list*): A list of times for which outputs will be found.
Default: None.
-* ``time_units`` (*str*): The time units used for requesting outputs by time.
- Default: '1' (code units).
-
* ``initial_redshift`` (*float*): The earliest redshift for outputs to be
included. If None, the initial redshift of the simulation is used. This
can be used in combination with either ``final_time`` or ``final_redshift``.
Default: None.
-* ``final_time`` (*float*): The latest redshift for outputs to be included.
+* ``final_redshift`` (*float*): The latest redshift for outputs to be included.
If None, the final redshift of the simulation is used. This can be used
in combination with either ``initial_time`` or ``initial_redshift``.
Default: None.
@@ -157,11 +151,11 @@
* ``initial_cycle`` (*float*): The earliest cycle for outputs to be
included. If None, the initial cycle of the simulation is used. This can
- only be used with final_cycle. Default: None.
+ only be used with final_cycle. Default: None. (Enzo only)
* ``final_cycle`` (*float*): The latest cycle for outputs to be included.
If None, the final cycle of the simulation is used. This can only be used
- in combination with initial_cycle. Default: None.
+ in combination with initial_cycle. Default: None. (Enzo only)
* ``tolerance`` (*float*): Used in combination with ``times`` or ``redshifts``
keywords, this is the tolerance within which outputs are accepted given
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 doc/source/examining/loading_data.rst
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -469,6 +469,8 @@
first image in the primary file. If this is not the case,
yt will raise a warning and will not load this field.
+.. _additional_fits_options:
+
Additional Options
^^^^^^^^^^^^^^^^^^
@@ -570,6 +572,35 @@
``WCSAxes`` is still in an experimental state, but as its functionality improves it will be
utilized more here.
+``create_spectral_slabs``
+"""""""""""""""""""""""""
+
+.. note::
+
+ The following functionality requires the `spectral-cube <http://spectral-cube.readthedocs.org>`_
+ library to be installed.
+
+If you have a spectral intensity dataset of some sort, and would like to extract emission in
+particular slabs along the spectral axis of a certain width, ``create_spectral_slabs`` can be
+used to generate a dataset with these slabs as different fields. In this example, we use it
+to extract individual lines from an intensity cube:
+
+.. code-block:: python
+
+ slab_centers = {'13CN': (218.03117, 'GHz'),
+ 'CH3CH2CHO': (218.284256, 'GHz'),
+ 'CH3NH2': (218.40956, 'GHz')}
+ slab_width = (0.05, "GHz")
+ ds = create_spectral_slabs("intensity_cube.fits",
+ slab_centers, slab_width,
+ nan_mask=0.0)
+
+All keyword arguments to `create_spectral_slabs` are passed on to `load` when creating the dataset
+(see :ref:`additional_fits_options` above). In the returned dataset, the different slabs will be
+different fields, with the field names taken from the keys in ``slab_centers``. The WCS coordinates
+on the spectral axis are reset so that the center of the domain along this axis is zero, and the
+left and right edges of the domain along this axis are :math:`\pm` ``0.5*slab_width``.
+
Examples of Using FITS Data
^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -635,13 +666,14 @@
import yt
ds = yt.load("snapshot_061.hdf5")
-However, yt cannot detect raw-binary Gadget data, and so you must specify the
-format as being Gadget:
+Gadget data in raw binary format can also be loaded with the ``load`` command.
+This is only supported for snapshots created with the ``SnapFormat`` parameter
+set to 1 (the standard for Gadget-2).
.. code-block:: python
import yt
- ds = yt.GadgetDataset("snapshot_061")
+ ds = yt.load("snapshot_061")
.. _particle-bbox:
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 doc/source/installing.rst
--- a/doc/source/installing.rst
+++ b/doc/source/installing.rst
@@ -213,10 +213,31 @@
++++++++++++++++++++++++++++++++++++++
To install yt from source, you must make sure you have yt's dependencies
-installed on your system. These include: a C compiler, ``HDF5``, ``python``,
-``Cython``, ``NumPy``, ``matplotlib``, ``sympy``, and ``h5py``. From here, you
-can use ``pip`` (which comes with ``Python``) to install the latest stable
-version of yt:
+installed on your system.
+
+If you use a Linux OS, use your distro's package manager to install these yt
+dependencies on your system:
+
+- ``HDF5``
+- ``zeromq``
+- ``sqlite``
+- ``mercurial``
+
+Then install the required Python packages with ``pip``:
+
+.. code-block:: bash
+
+ $ pip install -r requirements.txt
+
+If you're using IPython notebooks, you can install its dependencies
+with ``pip`` as well:
+
+.. code-block:: bash
+
+ $ pip install -r optional-requirements.txt
+
+From here, you can use ``pip`` (which comes with ``Python``) to install the latest
+stable version of yt:
.. code-block:: bash
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 doc/source/visualizing/sketchfab.rst
--- a/doc/source/visualizing/sketchfab.rst
+++ b/doc/source/visualizing/sketchfab.rst
@@ -56,7 +56,7 @@
import yt
ds = yt.load("/data/workshop2012/IsolatedGalaxy/galaxy0030/galaxy0030")
- sphere = ds.sphere("max", (1.0, "mpc"))
+ sphere = ds.sphere("max", (1.0, "Mpc"))
surface = ds.surface(sphere, "density", 1e-27)
This object, ``surface``, can be queried for values on the surface. For
@@ -172,7 +172,7 @@
trans = [1.0, 0.5]
filename = './surfaces'
- sphere = ds.sphere("max", (1.0, "mpc"))
+ sphere = ds.sphere("max", (1.0, "Mpc"))
for i,r in enumerate(rho):
surf = ds.surface(sphere, 'density', r)
surf.export_obj(filename, transparency = trans[i], color_field='temperature', plot_index = i)
@@ -248,7 +248,7 @@
return (data['density']*data['density']*np.sqrt(data['temperature']))
add_field("emissivity", function=_Emissivity, units=r"g*K/cm**6")
- sphere = ds.sphere("max", (1.0, "mpc"))
+ sphere = ds.sphere("max", (1.0, "Mpc"))
for i,r in enumerate(rho):
surf = ds.surface(sphere, 'density', r)
surf.export_obj(filename, transparency = trans[i],
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 optional-requirements.txt
--- /dev/null
+++ b/optional-requirements.txt
@@ -0,0 +1,1 @@
+ipython[notebook]
\ No newline at end of file
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 requirements.txt
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,6 @@
+numpy==1.9.2
+matplotlib==1.4.3
+Cython==0.22
+h5py==2.5.0
+nose==1.3.6
+sympy==0.7.6
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -127,8 +127,8 @@
field_units = {"dl": "cm", "redshift": "", "temperature": "K"}
field_data = {}
if use_peculiar_velocity:
- input_fields.append('los_velocity')
- field_units["los_velocity"] = "cm/s"
+ input_fields.append('velocity_los')
+ field_units["velocity_los"] = "cm/s"
for feature in self.line_list + self.continuum_list:
if not feature['field_name'] in input_fields:
input_fields.append(feature['field_name'])
@@ -171,7 +171,7 @@
if use_peculiar_velocity:
# include factor of (1 + z) because our velocity is in proper frame.
delta_lambda += continuum['wavelength'] * (1 + field_data['redshift']) * \
- field_data['los_velocity'] / speed_of_light_cgs
+ field_data['velocity_los'] / speed_of_light_cgs
this_wavelength = delta_lambda + continuum['wavelength']
right_index = np.digitize(this_wavelength, self.lambda_bins).clip(0, self.n_lambda)
left_index = np.digitize((this_wavelength *
@@ -208,7 +208,7 @@
if use_peculiar_velocity:
# include factor of (1 + z) because our velocity is in proper frame.
delta_lambda += line['wavelength'] * (1 + field_data['redshift']) * \
- field_data['los_velocity'] / speed_of_light_cgs
+ field_data['velocity_los'] / speed_of_light_cgs
thermal_b = km_per_cm * np.sqrt((2 * boltzmann_constant_cgs *
field_data['temperature']) /
(amu_cgs * line['atomic_mass']))
@@ -260,7 +260,7 @@
if line['label_threshold'] is not None and \
column_density[lixel] >= line['label_threshold']:
if use_peculiar_velocity:
- peculiar_velocity = km_per_cm * field_data['los_velocity'][lixel]
+ peculiar_velocity = km_per_cm * field_data['velocity_los'][lixel]
else:
peculiar_velocity = 0.0
self.spectrum_line_list.append({'label': line['label'],
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
--- a/yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
+++ b/yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
@@ -33,6 +33,13 @@
class LightRay(CosmologySplice):
"""
+ LightRay(parameter_filename, simulation_type=None,
+ near_redshift=None, far_redshift=None,
+ use_minimum_datasets=True, deltaz_min=0.0,
+ minimum_coherent_box_fraction=0.0,
+ time_data=True, redshift_data=True,
+ find_outputs=False, load_kwargs=None):
+
Create a LightRay object. A light ray is much like a light cone,
in that it stacks together multiple datasets in order to extend a
redshift interval. Unlike a light cone, which does randomly
@@ -94,6 +101,12 @@
Whether or not to search for datasets in the current
directory.
Default: False.
+ load_kwargs : optional, dict
+ Optional dictionary of kwargs to be passed to the "load"
+ function, appropriate for use of certain frontends. E.g.
+ Tipsy using "bounding_box"
+ Gadget using "unit_base", etc.
+ Default : None
"""
def __init__(self, parameter_filename, simulation_type=None,
@@ -101,7 +114,7 @@
use_minimum_datasets=True, deltaz_min=0.0,
minimum_coherent_box_fraction=0.0,
time_data=True, redshift_data=True,
- find_outputs=False):
+ find_outputs=False, load_kwargs=None):
self.near_redshift = near_redshift
self.far_redshift = far_redshift
@@ -109,13 +122,16 @@
self.deltaz_min = deltaz_min
self.minimum_coherent_box_fraction = minimum_coherent_box_fraction
self.parameter_filename = parameter_filename
-
+ if load_kwargs is None:
+ self.load_kwargs = {}
+ else:
+ self.load_kwargs = load_kwargs
self.light_ray_solution = []
self._data = {}
# Make a light ray from a single, given dataset.
if simulation_type is None:
- ds = load(parameter_filename)
+ ds = load(parameter_filename, **self.load_kwargs)
if ds.cosmological_simulation:
redshift = ds.current_redshift
self.cosmology = Cosmology(
@@ -243,6 +259,12 @@
get_los_velocity=True, redshift=None,
njobs=-1):
"""
+ make_light_ray(seed=None, start_position=None, end_position=None,
+ trajectory=None, fields=None, setup_function=None,
+ solution_filename=None, data_filename=None,
+ get_los_velocity=True, redshift=None,
+ njobs=-1)
+
Create a light ray and get field values for each lixel. A light
ray consists of a list of field values for cells intersected by
the ray and the path length of the ray through those cells.
@@ -343,9 +365,9 @@
all_fields = fields[:]
all_fields.extend(['dl', 'dredshift', 'redshift'])
if get_los_velocity:
- all_fields.extend(['x-velocity', 'y-velocity',
- 'z-velocity', 'los_velocity'])
- data_fields.extend(['x-velocity', 'y-velocity', 'z-velocity'])
+ all_fields.extend(['velocity_x', 'velocity_y',
+ 'velocity_z', 'velocity_los'])
+ data_fields.extend(['velocity_x', 'velocity_y', 'velocity_z'])
all_ray_storage = {}
for my_storage, my_segment in parallel_objects(self.light_ray_solution,
@@ -353,7 +375,7 @@
njobs=njobs):
# Load dataset for segment.
- ds = load(my_segment['filename'])
+ ds = load(my_segment['filename'], **self.load_kwargs)
my_segment['unique_identifier'] = ds.unique_identifier
if redshift is not None:
@@ -364,11 +386,15 @@
if setup_function is not None:
setup_function(ds)
-
- my_segment["start"] = ds.domain_width * my_segment["start"] + \
- ds.domain_left_edge
- my_segment["end"] = ds.domain_width * my_segment["end"] + \
- ds.domain_left_edge
+
+ if start_position is not None:
+ my_segment["start"] = ds.arr(my_segment["start"], "code_length")
+ my_segment["end"] = ds.arr(my_segment["end"], "code_length")
+ else:
+ my_segment["start"] = ds.domain_width * my_segment["start"] + \
+ ds.domain_left_edge
+ my_segment["end"] = ds.domain_width * my_segment["end"] + \
+ ds.domain_left_edge
if not ds.cosmological_simulation:
next_redshift = my_segment["redshift"]
@@ -412,10 +438,10 @@
if get_los_velocity:
line_of_sight = sub_segment[1] - sub_segment[0]
line_of_sight /= ((line_of_sight**2).sum())**0.5
- sub_vel = ds.arr([sub_ray['x-velocity'],
- sub_ray['y-velocity'],
- sub_ray['z-velocity']])
- sub_data['los_velocity'].extend((np.rollaxis(sub_vel, 1) *
+ sub_vel = ds.arr([sub_ray['velocity_x'],
+ sub_ray['velocity_y'],
+ sub_ray['velocity_z']])
+ sub_data['velocity_los'].extend((np.rollaxis(sub_vel, 1) *
line_of_sight).sum(axis=1)[asort])
del sub_vel
@@ -423,7 +449,6 @@
del sub_ray, asort
for key in sub_data:
- if key in "xyz": continue
sub_data[key] = ds.arr(sub_data[key]).in_cgs()
# Get redshift for each lixel. Assume linear relation between l and z.
@@ -461,18 +486,32 @@
@parallel_root_only
def _write_light_ray(self, filename, data):
- "Write light ray data to hdf5 file."
+ """
+ _write_light_ray(filename, data)
+
+ Write light ray data to hdf5 file.
+ """
mylog.info("Saving light ray data to %s." % filename)
output = h5py.File(filename, 'w')
for field in data.keys():
- output.create_dataset(field, data=data[field])
- output[field].attrs["units"] = str(data[field].units)
+ # if the field is a tuple, only use the second part of the tuple
+ # in the hdf5 output (i.e. ('gas', 'density') -> 'density')
+ if isinstance(field, tuple):
+ fieldname = field[1]
+ else:
+ fieldname = field
+ output.create_dataset(fieldname, data=data[field])
+ output[fieldname].attrs["units"] = str(data[field].units)
output.close()
@parallel_root_only
def _write_light_ray_solution(self, filename, extra_info=None):
- "Write light ray solution to a file."
+ """
+ _write_light_ray_solution(filename, extra_info=None)
+
+ Write light ray solution to a file.
+ """
mylog.info("Writing light ray solution to %s." % filename)
f = open(filename, 'w')
@@ -490,7 +529,11 @@
f.close()
def _flatten_dict_list(data, exceptions=None):
- "Flatten the list of dicts into one dict."
+ """
+ _flatten_dict_list(data, exceptions=None)
+
+ Flatten the list of dicts into one dict.
+ """
if exceptions is None: exceptions = []
new_data = {}
@@ -505,12 +548,20 @@
return new_data
def vector_length(start, end):
- "Calculate vector length."
+ """
+ vector_length(start, end)
+
+ Calculate vector length.
+ """
return np.sqrt(np.power((end - start), 2).sum())
def periodic_distance(coord1, coord2):
- "Calculate length of shortest vector between to points in periodic domain."
+ """
+ periodic_distance(coord1, coord2)
+
+ Calculate length of shortest vector between to points in periodic domain.
+ """
dif = coord1 - coord2
dim = np.ones(coord1.shape,dtype=int)
@@ -524,6 +575,8 @@
def periodic_ray(start, end, left=None, right=None):
"""
+ periodic_ray(start, end, left=None, right=None)
+
Break up periodic ray into non-periodic segments.
Accepts start and end points of periodic ray as YTArrays.
Accepts optional left and right edges of periodic volume as YTArrays.
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/analysis_modules/halo_finding/halo_objects.py
--- a/yt/analysis_modules/halo_finding/halo_objects.py
+++ b/yt/analysis_modules/halo_finding/halo_objects.py
@@ -1232,9 +1232,8 @@
fglob = path.join(basedir, 'halos_%d.*.bin' % n)
files = glob.glob(fglob)
halos = self._get_halos_binary(files)
- #Jc = mass_sun_cgs/ ds['mpchcm'] * 1e5
Jc = 1.0
- length = 1.0 / ds['mpchcm']
+ length = 1.0 / ds['Mpchcm']
conv = dict(pos = np.array([length, length, length,
1, 1, 1]), # to unitary
r=1.0/ds['kpchcm'], # to unitary
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -729,7 +729,7 @@
>>> import yt
>>> ds = yt.load("RedshiftOutput0005")
- >>> sp = ds.sphere("max", (1.0, 'mpc'))
+ >>> sp = ds.sphere("max", (1.0, 'Mpc'))
>>> cr = ds.cut_region(sp, ["obj['temperature'] < 1e3"])
"""
_type_name = "cut_region"
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/data_objects/time_series.py
--- a/yt/data_objects/time_series.py
+++ b/yt/data_objects/time_series.py
@@ -406,6 +406,7 @@
self.basename = os.path.basename(parameter_filename)
self.directory = os.path.dirname(parameter_filename)
self.parameters = {}
+ self.key_parameters = []
# Set some parameter defaults.
self._set_parameter_defaults()
@@ -420,6 +421,21 @@
self.print_key_parameters()
+ def _set_parameter_defaults(self):
+ pass
+
+ def _parse_parameter_file(self):
+ pass
+
+ def _set_units(self):
+ pass
+
+ def _calculate_simulation_bounds(self):
+ pass
+
+ def _get_all_outputs(**kwargs):
+ pass
+
def __repr__(self):
return self.parameter_filename
@@ -445,23 +461,78 @@
"""
Print out some key parameters for the simulation.
"""
- for a in ["domain_dimensions", "domain_left_edge",
- "domain_right_edge", "initial_time", "final_time",
- "stop_cycle", "cosmological_simulation"]:
- if not hasattr(self, a):
- mylog.error("Missing %s in dataset definition!", a)
- continue
- v = getattr(self, a)
- mylog.info("Parameters: %-25s = %s", a, v)
- if hasattr(self, "cosmological_simulation") and \
- getattr(self, "cosmological_simulation"):
+ if self.simulation_type == "grid":
+ for a in ["domain_dimensions", "domain_left_edge",
+ "domain_right_edge"]:
+ self._print_attr(a)
+ for a in ["initial_time", "final_time",
+ "cosmological_simulation"]:
+ self._print_attr(a)
+ if getattr(self, "cosmological_simulation", False):
for a in ["box_size", "omega_lambda",
"omega_matter", "hubble_constant",
"initial_redshift", "final_redshift"]:
- if not hasattr(self, a):
- mylog.error("Missing %s in dataset definition!", a)
- continue
- v = getattr(self, a)
- mylog.info("Parameters: %-25s = %s", a, v)
+ self._print_attr(a)
+ for a in self.key_parameters:
+ self._print_attr(a)
mylog.info("Total datasets: %d." % len(self.all_outputs))
+ def _print_attr(self, a):
+ """
+ Print the attribute or warn about it missing.
+ """
+ if not hasattr(self, a):
+ mylog.error("Missing %s in dataset definition!", a)
+ return
+ v = getattr(self, a)
+ mylog.info("Parameters: %-25s = %s", a, v)
+
+ def _get_outputs_by_key(self, key, values, tolerance=None, outputs=None):
+ r"""
+ Get datasets at or near to given values.
+
+ Parameters
+ ----------
+ key: str
+ The key by which to retrieve outputs, usually 'time' or
+ 'redshift'.
+ values: array_like
+ A list of values, given as floats.
+ tolerance : float
+ If not None, do not return a dataset unless the value is
+ within the tolerance value. If None, simply return the
+ nearest dataset.
+ Default: None.
+ outputs : list
+ The list of outputs from which to choose. If None,
+ self.all_outputs is used.
+ Default: None.
+
+ Examples
+ --------
+ >>> datasets = es.get_outputs_by_key('redshift', [0, 1, 2], tolerance=0.1)
+
+ """
+
+ if not isinstance(values, YTArray):
+ if isinstance(values, tuple) and len(values) == 2:
+ values = self.arr(*values)
+ else:
+ values = self.arr(values)
+ values = values.in_cgs()
+
+ if outputs is None:
+ outputs = self.all_outputs
+ my_outputs = []
+ if not outputs:
+ return my_outputs
+ for value in values:
+ outputs.sort(key=lambda obj:np.abs(value - obj[key]))
+ if (tolerance is None or np.abs(value - outputs[0][key]) <= tolerance) \
+ and outputs[0] not in my_outputs:
+ my_outputs.append(outputs[0])
+ else:
+ mylog.error("No dataset added for %s = %f.", key, value)
+
+ outputs.sort(key=lambda obj: obj['time'])
+ return my_outputs
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/fields/field_aliases.py
--- a/yt/fields/field_aliases.py
+++ b/yt/fields/field_aliases.py
@@ -141,12 +141,12 @@
("CellMassCode", "code_mass"),
("TotalMassMsun", "msun"),
("CellVolumeCode", "code_length"),
- ("CellVolumeMpc", "mpc**3"),
- ("ParticleSpecificAngularMomentumXKMSMPC","km/s/mpc"),
- ("ParticleSpecificAngularMomentumYKMSMPC","km/s/mpc"),
- ("ParticleSpecificAngularMomentumZKMSMPC","km/s/mpc"),
- ("RadiusMpc", "mpc"),
- ("ParticleRadiusMpc", "mpc"),
+ ("CellVolumeMpc", "Mpc**3"),
+ ("ParticleSpecificAngularMomentumXKMSMPC","km/s/Mpc"),
+ ("ParticleSpecificAngularMomentumYKMSMPC","km/s/Mpc"),
+ ("ParticleSpecificAngularMomentumZKMSMPC","km/s/Mpc"),
+ ("RadiusMpc", "Mpc"),
+ ("ParticleRadiusMpc", "Mpc"),
("ParticleRadiuskpc", "kpc"),
("Radiuskpc", "kpc"),
("ParticleRadiuskpch", "kpc"),
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/frontends/enzo/data_structures.py
--- a/yt/frontends/enzo/data_structures.py
+++ b/yt/frontends/enzo/data_structures.py
@@ -890,7 +890,7 @@
elif self.dimensionality == 2:
self._setup_2d()
- def set_code_units(self):
+ def _set_code_unit_attributes(self):
if self.cosmological_simulation:
k = self.cosmology_get_units()
# Now some CGS values
@@ -928,17 +928,6 @@
magnetic_unit = np.float64(magnetic_unit.in_cgs())
self.magnetic_unit = self.quan(magnetic_unit, "gauss")
- self._override_code_units()
-
- self.unit_registry.modify("code_magnetic", self.magnetic_unit)
- self.unit_registry.modify("code_length", self.length_unit)
- self.unit_registry.modify("code_mass", self.mass_unit)
- self.unit_registry.modify("code_time", self.time_unit)
- self.unit_registry.modify("code_velocity", self.velocity_unit)
- DW = self.arr(self.domain_right_edge - self.domain_left_edge, "code_length")
- self.unit_registry.add("unitary", float(DW.max() * DW.units.base_value),
- DW.units.dimensions)
-
def cosmology_get_units(self):
"""
Return an Enzo-fortran style dictionary of units to feed into custom
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/frontends/enzo/simulation_handling.py
--- a/yt/frontends/enzo/simulation_handling.py
+++ b/yt/frontends/enzo/simulation_handling.py
@@ -13,36 +13,34 @@
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
-from yt.funcs import *
-
import numpy as np
import glob
import os
from yt.convenience import \
- load
+ load, \
+ only_on_root
from yt.data_objects.time_series import \
SimulationTimeSeries, DatasetSeries
from yt.units import dimensions
from yt.units.unit_registry import \
- UnitRegistry
+ UnitRegistry
from yt.units.yt_array import \
- YTArray, YTQuantity
+ YTArray, YTQuantity
from yt.utilities.cosmology import \
Cosmology
-from yt.utilities.definitions import \
- sec_conversion
from yt.utilities.exceptions import \
InvalidSimulationTimeSeries, \
MissingParameter, \
NoStoppingCondition
+from yt.utilities.logger import ytLogger as \
+ mylog
from yt.utilities.parallel_tools.parallel_analysis_interface import \
parallel_objects
-from yt.utilities.physical_constants import \
- gravitational_constant_cgs as G
-
+
class EnzoSimulation(SimulationTimeSeries):
- r"""Initialize an Enzo Simulation object.
+ r"""
+ Initialize an Enzo Simulation object.
Upon creation, the parameter file is parsed and the time and redshift
are calculated and stored in all_outputs. A time units dictionary is
@@ -63,14 +61,8 @@
Examples
--------
- >>> from yt.mods import *
- >>> es = EnzoSimulation("my_simulation.par")
- >>> es.get_time_series()
- >>> for ds in es:
- ... print ds.current_time
-
- >>> from yt.mods import *
- >>> es = simulation("my_simulation.par", "Enzo")
+ >>> import yt
+ >>> es = yt.simulation("my_simulation.par", "Enzo")
>>> es.get_time_series()
>>> for ds in es:
... print ds.current_time
@@ -78,7 +70,8 @@
"""
def __init__(self, parameter_filename, find_outputs=False):
-
+ self.simulation_type = "grid"
+ self.key_parameters = ["stop_cycle"]
SimulationTimeSeries.__init__(self, parameter_filename,
find_outputs=find_outputs)
@@ -87,14 +80,14 @@
self.unit_registry.lut["code_time"] = (1.0, dimensions.time)
if self.cosmological_simulation:
# Instantiate EnzoCosmology object for units and time conversions.
- self.enzo_cosmology = \
+ self.cosmology = \
EnzoCosmology(self.parameters['CosmologyHubbleConstantNow'],
self.parameters['CosmologyOmegaMatterNow'],
self.parameters['CosmologyOmegaLambdaNow'],
0.0, self.parameters['CosmologyInitialRedshift'],
unit_registry=self.unit_registry)
- self.time_unit = self.enzo_cosmology.time_unit.in_units("s")
+ self.time_unit = self.cosmology.time_unit.in_units("s")
self.unit_registry.modify("h", self.hubble_constant)
# Comoving lengths
for my_unit in ["m", "pc", "AU", "au"]:
@@ -160,7 +153,7 @@
used in combination with either final_time or
final_redshift.
Default: None.
- final_time : float
+ final_redshift : float
The latest redshift for outputs to be included. If None,
the final redshift of the simulation is used. This can be
used in combination with either initial_time or
@@ -197,8 +190,8 @@
Examples
--------
- >>> from yt.mods import *
- >>> es = simulation("my_simulation.par", "Enzo")
+ >>> import yt
+ >>> es = yt.simulation("my_simulation.par", "Enzo")
>>> es.get_time_series(initial_redshift=10, final_time=(13.7, "Gyr"),
redshift_data=False)
@@ -207,8 +200,6 @@
>>> es.get_time_series(final_cycle=100000)
- >>> es.get_time_series(find_outputs=True)
-
>>> # after calling get_time_series
>>> for ds in es.piter():
... p = ProjectionPlot(ds, 'x', "density")
@@ -226,7 +217,9 @@
if (initial_redshift is not None or \
final_redshift is not None) and \
not self.cosmological_simulation:
- raise InvalidSimulationTimeSeries('An initial or final redshift has been given for a noncosmological simulation.')
+ raise InvalidSimulationTimeSeries(
+ "An initial or final redshift has been given for a " +
+ "noncosmological simulation.")
if time_data and redshift_data:
my_all_outputs = self.all_outputs
@@ -244,12 +237,14 @@
# Apply selection criteria to the set.
if times is not None:
- my_outputs = self._get_outputs_by_time(times, tolerance=tolerance,
- outputs=my_all_outputs)
+ my_outputs = self._get_outputs_by_key("time", times,
+ tolerance=tolerance,
+ outputs=my_all_outputs)
elif redshifts is not None:
- my_outputs = self._get_outputs_by_redshift(redshifts, tolerance=tolerance,
- outputs=my_all_outputs)
+ my_outputs = self._get_outputs_by_key("redshift", redshifts,
+ tolerance=tolerance,
+ outputs=my_all_outputs)
elif initial_cycle is not None or final_cycle is not None:
if initial_cycle is None:
@@ -272,9 +267,11 @@
elif isinstance(initial_time, tuple) and len(initial_time) == 2:
initial_time = self.quan(*initial_time)
elif not isinstance(initial_time, YTArray):
- raise RuntimeError("Error: initial_time must be given as a float or tuple of (value, units).")
+ raise RuntimeError(
+ "Error: initial_time must be given as a float or " +
+ "tuple of (value, units).")
elif initial_redshift is not None:
- my_initial_time = self.enzo_cosmology.t_from_z(initial_redshift)
+ my_initial_time = self.cosmology.t_from_z(initial_redshift)
else:
my_initial_time = self.initial_time
@@ -284,10 +281,12 @@
elif isinstance(final_time, tuple) and len(final_time) == 2:
final_time = self.quan(*final_time)
elif not isinstance(final_time, YTArray):
- raise RuntimeError("Error: final_time must be given as a float or tuple of (value, units).")
+ raise RuntimeError(
+ "Error: final_time must be given as a float or " +
+ "tuple of (value, units).")
my_final_time = final_time.in_units("s")
elif final_redshift is not None:
- my_final_time = self.enzo_cosmology.t_from_z(final_redshift)
+ my_final_time = self.cosmology.t_from_z(final_redshift)
else:
my_final_time = self.final_time
@@ -390,8 +389,9 @@
raise MissingParameter(self.parameter_filename, v)
setattr(self, a, self.parameters[v])
else:
+ self.cosmological_simulation = 0
self.omega_lambda = self.omega_matter = \
- self.hubble_constant = self.cosmological_simulation = 0.0
+ self.hubble_constant = 0.0
# make list of redshift outputs
self.all_redshift_outputs = []
@@ -405,16 +405,10 @@
del output['index']
self.all_redshift_outputs = redshift_outputs
- def _calculate_redshift_dump_times(self):
- "Calculates time from redshift of redshift outputs."
-
- if not self.cosmological_simulation: return
- for output in self.all_redshift_outputs:
- output['time'] = self.enzo_cosmology.t_from_z(output['redshift'])
- self.all_redshift_outputs.sort(key=lambda obj:obj['time'])
-
def _calculate_time_outputs(self):
- "Calculate time outputs and their redshifts if cosmological."
+ """
+ Calculate time outputs and their redshifts if cosmological.
+ """
self.all_time_outputs = []
if self.final_time is None or \
@@ -432,7 +426,7 @@
output = {'index': index, 'filename': filename, 'time': current_time.copy()}
output['time'] = min(output['time'], self.final_time)
if self.cosmological_simulation:
- output['redshift'] = self.enzo_cosmology.z_from_t(current_time)
+ output['redshift'] = self.cosmology.z_from_t(current_time)
self.all_time_outputs.append(output)
if np.abs(self.final_time - current_time) / self.final_time < 1e-4: break
@@ -440,7 +434,9 @@
index += 1
def _calculate_cycle_outputs(self):
- "Calculate cycle outputs."
+ """
+ Calculate cycle outputs.
+ """
mylog.warn('Calculating cycle outputs. Dataset times will be unavailable.')
@@ -460,7 +456,9 @@
index += 1
def _get_all_outputs(self, find_outputs=False):
- "Get all potential datasets and combine into a time-sorted list."
+ """
+ Get all potential datasets and combine into a time-sorted list.
+ """
# Create the set of outputs from which further selection will be done.
if find_outputs:
@@ -468,8 +466,12 @@
elif self.parameters['dtDataDump'] > 0 and \
self.parameters['CycleSkipDataDump'] > 0:
- mylog.info("Simulation %s has both dtDataDump and CycleSkipDataDump set.", self.parameter_filename )
- mylog.info(" Unable to calculate datasets. Attempting to search in the current directory")
+ mylog.info(
+ "Simulation %s has both dtDataDump and CycleSkipDataDump set.",
+ self.parameter_filename )
+ mylog.info(
+ " Unable to calculate datasets. " +
+ "Attempting to search in the current directory")
self._find_outputs()
else:
@@ -480,7 +482,10 @@
self._calculate_time_outputs()
# Calculate times for redshift outputs.
- self._calculate_redshift_dump_times()
+ if self.cosmological_simulation:
+ for output in self.all_redshift_outputs:
+ output["time"] = self.cosmology.t_from_z(output["redshift"])
+ self.all_redshift_outputs.sort(key=lambda obj:obj["time"])
self.all_outputs = self.all_time_outputs + self.all_redshift_outputs
if self.parameters['CycleSkipDataDump'] <= 0:
@@ -496,9 +501,9 @@
# Convert initial/final redshifts to times.
if self.cosmological_simulation:
- self.initial_time = self.enzo_cosmology.t_from_z(self.initial_redshift)
+ self.initial_time = self.cosmology.t_from_z(self.initial_redshift)
self.initial_time.units.registry = self.unit_registry
- self.final_time = self.enzo_cosmology.t_from_z(self.final_redshift)
+ self.final_time = self.cosmology.t_from_z(self.final_redshift)
self.final_time.units.registry = self.unit_registry
# If not a cosmology simulation, figure out the stopping criteria.
@@ -516,11 +521,15 @@
'StopCycle' in self.parameters):
raise NoStoppingCondition(self.parameter_filename)
if self.final_time is None:
- mylog.warn('Simulation %s has no stop time set, stopping condition will be based only on cycles.',
- self.parameter_filename)
+ mylog.warn(
+ "Simulation %s has no stop time set, stopping condition " +
+ "will be based only on cycles.",
+ self.parameter_filename)
def _set_parameter_defaults(self):
- "Set some default parameters to avoid problems if they are not in the parameter file."
+ """
+ Set some default parameters to avoid problems if they are not in the parameter file.
+ """
self.parameters['GlobalDir'] = self.directory
self.parameters['DataDumpName'] = "data"
@@ -570,7 +579,9 @@
self.final_redshift = self.all_outputs[-1]['redshift']
def _check_for_outputs(self, potential_outputs):
- r"""Check a list of files to see if they are valid datasets."""
+ """
+ Check a list of files to see if they are valid datasets.
+ """
only_on_root(mylog.info, "Checking %d potential outputs.",
len(potential_outputs))
@@ -603,112 +614,10 @@
return my_outputs
- def _get_outputs_by_key(self, key, values, tolerance=None, outputs=None):
- r"""Get datasets at or near to given values.
-
- Parameters
- ----------
- key: str
- The key by which to retrieve outputs, usually 'time' or
- 'redshift'.
- values: array_like
- A list of values, given as floats.
- tolerance : float
- If not None, do not return a dataset unless the value is
- within the tolerance value. If None, simply return the
- nearest dataset.
- Default: None.
- outputs : list
- The list of outputs from which to choose. If None,
- self.all_outputs is used.
- Default: None.
-
- Examples
- --------
- >>> datasets = es.get_outputs_by_key('redshift', [0, 1, 2], tolerance=0.1)
-
- """
-
- if not isinstance(values, np.ndarray):
- values = ensure_list(values)
- if outputs is None:
- outputs = self.all_outputs
- my_outputs = []
- if not outputs:
- return my_outputs
- for value in values:
- outputs.sort(key=lambda obj:np.abs(value - obj[key]))
- if (tolerance is None or np.abs(value - outputs[0][key]) <= tolerance) \
- and outputs[0] not in my_outputs:
- my_outputs.append(outputs[0])
- else:
- mylog.error("No dataset added for %s = %f.", key, value)
-
- outputs.sort(key=lambda obj: obj['time'])
- return my_outputs
-
- def _get_outputs_by_redshift(self, redshifts, tolerance=None, outputs=None):
- r"""Get datasets at or near to given redshifts.
-
- Parameters
- ----------
- redshifts: array_like
- A list of redshifts, given as floats.
- tolerance : float
- If not None, do not return a dataset unless the value is
- within the tolerance value. If None, simply return the
- nearest dataset.
- Default: None.
- outputs : list
- The list of outputs from which to choose. If None,
- self.all_outputs is used.
- Default: None.
-
- Examples
- --------
- >>> datasets = es.get_outputs_by_redshift([0, 1, 2], tolerance=0.1)
-
- """
-
- return self._get_outputs_by_key('redshift', redshifts, tolerance=tolerance,
- outputs=outputs)
-
- def _get_outputs_by_time(self, times, tolerance=None, outputs=None):
- r"""Get datasets at or near to given times.
-
- Parameters
- ----------
- times: tuple of type (float array, str)
- A list of times for which outputs will be found and the units
- of those values. For example, ([0, 1, 2, 3], "s").
- tolerance : float
- If not None, do not return a dataset unless the time is
- within the tolerance value. If None, simply return the
- nearest dataset.
- Default = None.
- outputs : list
- The list of outputs from which to choose. If None,
- self.all_outputs is used.
- Default: None.
-
- Examples
- --------
- >>> datasets = es.get_outputs_by_time([600, 500, 400], tolerance=10.)
-
- """
-
- if not isinstance(times, YTArray):
- if isinstance(times, tuple) and len(times) == 2:
- times = self.arr(*times)
- else:
- times = self.arr(times, "code_time")
- times = times.in_units("s")
- return self._get_outputs_by_key('time', times, tolerance=tolerance,
- outputs=outputs)
-
def _write_cosmology_outputs(self, filename, outputs, start_index,
decimals=3):
- r"""Write cosmology output parameters for a cosmology splice.
+ """
+ Write cosmology output parameters for a cosmology splice.
"""
mylog.info("Writing redshift output list to %s.", filename)
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/frontends/fits/data_structures.py
--- a/yt/frontends/fits/data_structures.py
+++ b/yt/frontends/fits/data_structures.py
@@ -165,7 +165,7 @@
units = self._determine_image_units(hdu.header, known_units)
try:
# Grab field name from btype
- fname = hdu.header["btype"].lower()
+ fname = hdu.header["btype"]
except KeyError:
# Try to guess the name from the units
fname = self._guess_name_from_units(units)
@@ -205,18 +205,6 @@
"the same dimensions as the primary and will not be " +
"available as a field.")
- # For line fields, we still read the primary field. Not sure how to extend this
- # For now, we pick off the first field from the field list.
- line_db = self.dataset.line_database
- primary_fname = self.field_list[0][1]
- for k, v in iteritems(line_db):
- mylog.info("Adding line field: %s at frequency %g GHz" % (k, v))
- self.field_list.append((self.dataset_type, k))
- self._ext_map[k] = self._ext_map[primary_fname]
- self._axis_map[k] = self._axis_map[primary_fname]
- self._file_map[k] = self._file_map[primary_fname]
- self.dataset.field_units[k] = self.dataset.field_units[primary_fname]
-
def _count_grids(self):
self.num_grids = self.ds.parameters["nprocs"]
@@ -242,19 +230,11 @@
bbox = np.array([[le,re] for le, re in zip(ds.domain_left_edge,
ds.domain_right_edge)])
dims = np.array(ds.domain_dimensions)
- # If we are creating a dataset of lines, only decompose along the position axes
- if len(ds.line_database) > 0:
- dims[ds.spec_axis] = 1
psize = get_psize(dims, self.num_grids)
gle, gre, shapes, slices = decompose_array(dims, psize, bbox)
self.grid_left_edge = self.ds.arr(gle, "code_length")
self.grid_right_edge = self.ds.arr(gre, "code_length")
self.grid_dimensions = np.array([shape for shape in shapes], dtype="int32")
- # If we are creating a dataset of lines, only decompose along the position axes
- if len(ds.line_database) > 0:
- self.grid_left_edge[:,ds.spec_axis] = ds.domain_left_edge[ds.spec_axis]
- self.grid_right_edge[:,ds.spec_axis] = ds.domain_right_edge[ds.spec_axis]
- self.grid_dimensions[:,ds.spec_axis] = ds.domain_dimensions[ds.spec_axis]
else:
self.grid_left_edge[0,:] = ds.domain_left_edge
self.grid_right_edge[0,:] = ds.domain_right_edge
@@ -322,8 +302,6 @@
nan_mask=None,
spectral_factor=1.0,
z_axis_decomp=False,
- line_database=None,
- line_width=None,
suppress_astropy_warnings=True,
parameters=None,
units_override=None):
@@ -336,19 +314,6 @@
self.z_axis_decomp = z_axis_decomp
self.spectral_factor = spectral_factor
- if line_width is not None:
- self.line_width = YTQuantity(line_width[0], line_width[1])
- self.line_units = line_width[1]
- mylog.info("For line folding, spectral_factor = 1.0")
- self.spectral_factor = 1.0
- else:
- self.line_width = None
-
- self.line_database = {}
- if line_database is not None:
- for k in line_database:
- self.line_database[k] = YTQuantity(line_database[k], self.line_units)
-
if suppress_astropy_warnings:
warnings.filterwarnings('ignore', module="astropy", append=True)
auxiliary_files = ensure_list(auxiliary_files)
@@ -361,13 +326,13 @@
self.nan_mask = {"all":nan_mask}
elif isinstance(nan_mask, dict):
self.nan_mask = nan_mask
- if isinstance(self.filenames[0], _astropy.pyfits.hdu.image._ImageBaseHDU):
- self._handle = FITSFileHandler(self.filenames[0])
- fn = "InMemoryFITSImage_%s" % (uuid.uuid4().hex)
+ self._handle = FITSFileHandler(self.filenames[0])
+ if (isinstance(self.filenames[0], _astropy.pyfits.hdu.image._ImageBaseHDU) or
+ isinstance(self.filenames[0], _astropy.pyfits.HDUList)):
+ fn = "InMemoryFITSFile_%s" % uuid.uuid4().hex
else:
- self._handle = FITSFileHandler(self.filenames[0])
fn = self.filenames[0]
- self._handle._fits_files = [self._handle]
+ self._handle._fits_files.append(self._handle)
if self.num_files > 1:
for fits_file in auxiliary_files:
if isinstance(fits_file, _astropy.pyfits.hdu.image._ImageBaseHDU):
@@ -540,20 +505,14 @@
# If nprocs is None, do some automatic decomposition of the domain
if self.specified_parameters["nprocs"] is None:
- if len(self.line_database) > 0:
- dims = 2
- else:
- dims = self.dimensionality
if self.z_axis_decomp:
nprocs = np.around(self.domain_dimensions[2]/8).astype("int")
else:
- nprocs = np.around(np.prod(self.domain_dimensions)/32**dims).astype("int")
+ nprocs = np.around(np.prod(self.domain_dimensions)/32**self.dimensionality).astype("int")
self.parameters["nprocs"] = max(min(nprocs, 512), 1)
else:
self.parameters["nprocs"] = self.specified_parameters["nprocs"]
- self.reversed = False
-
# Check to see if this data is in some kind of (Lat,Lon,Vel) format
self.spec_cube = False
x = 0
@@ -618,41 +577,23 @@
self._z0 = self.wcs.wcs.crval[self.spec_axis]
self.spec_unit = str(self.wcs.wcs.cunit[self.spec_axis])
- if self.line_width is not None:
- if self._dz < 0.0:
- self.reversed = True
- le = self.dims[self.spec_axis]+0.5
- else:
- le = 0.5
- self.line_width = self.line_width.in_units(self.spec_unit)
- self.freq_begin = (le-self._p0)*self._dz + self._z0
- # We now reset these so that they are consistent
- # with the new setup
- self._dz = np.abs(self._dz)
- self._p0 = 0.0
- self._z0 = 0.0
- nz = np.rint(self.line_width.value/self._dz).astype("int")
- self.line_width = self._dz*nz
- self.domain_left_edge[self.spec_axis] = -0.5*float(nz)
- self.domain_right_edge[self.spec_axis] = 0.5*float(nz)
- self.domain_dimensions[self.spec_axis] = nz
- else:
- if self.spectral_factor == "auto":
- self.spectral_factor = float(max(self.domain_dimensions[[self.lon_axis,
- self.lat_axis]]))
- self.spectral_factor /= self.domain_dimensions[self.spec_axis]
- mylog.info("Setting the spectral factor to %f" % (self.spectral_factor))
- Dz = self.domain_right_edge[self.spec_axis]-self.domain_left_edge[self.spec_axis]
- self.domain_right_edge[self.spec_axis] = self.domain_left_edge[self.spec_axis] + \
- self.spectral_factor*Dz
- self._dz /= self.spectral_factor
- self._p0 = (self._p0-0.5)*self.spectral_factor + 0.5
+ if self.spectral_factor == "auto":
+ self.spectral_factor = float(max(self.domain_dimensions[[self.lon_axis,
+ self.lat_axis]]))
+ self.spectral_factor /= self.domain_dimensions[self.spec_axis]
+ mylog.info("Setting the spectral factor to %f" % (self.spectral_factor))
+ Dz = self.domain_right_edge[self.spec_axis]-self.domain_left_edge[self.spec_axis]
+ self.domain_right_edge[self.spec_axis] = self.domain_left_edge[self.spec_axis] + \
+ self.spectral_factor*Dz
+ self._dz /= self.spectral_factor
+ self._p0 = (self._p0-0.5)*self.spectral_factor + 0.5
+
else:
self.wcs_2d = self.wcs
self.spec_axis = 2
self.spec_name = "z"
- self.spec_unit = "code length"
+ self.spec_unit = "code_length"
def spec2pixel(self, spec_value):
sv = self.arr(spec_value).in_units(self.spec_unit)
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/frontends/fits/io.py
--- a/yt/frontends/fits/io.py
+++ b/yt/frontends/fits/io.py
@@ -24,12 +24,6 @@
super(IOHandlerFITS, self).__init__(ds)
self.ds = ds
self._handle = ds._handle
- if self.ds.line_width is not None:
- self.line_db = self.ds.line_database
- self.dz = self.ds.line_width/self.domain_dimensions[self.ds.spec_axis]
- else:
- self.line_db = None
- self.dz = 1.
def _read_particles(self, fields_to_read, type, args, grid_list,
count_list, conv_factors):
@@ -79,32 +73,15 @@
dx = self.ds.domain_width/self.ds.domain_dimensions
for field in fields:
ftype, fname = field
- tmp_fname = fname
- if fname in self.ds.line_database:
- fname = self.ds.field_list[0][1]
f = self.ds.index._file_map[fname]
ds = f[self.ds.index._ext_map[fname]]
bzero, bscale = self.ds.index._scale_map[fname]
- fname = tmp_fname
ind = 0
for chunk in chunks:
for g in chunk.objs:
start = ((g.LeftEdge-self.ds.domain_left_edge)/dx).to_ndarray().astype("int")
end = start + g.ActiveDimensions
- if self.line_db is not None and fname in self.line_db:
- my_off = self.line_db.get(fname).in_units(self.ds.spec_unit).value
- my_off = my_off - 0.5*self.ds.line_width
- my_off = int((my_off-self.ds.freq_begin)/self.dz)
- my_off = max(my_off, 0)
- my_off = min(my_off, self.ds.dims[self.ds.spec_axis]-1)
- start[self.ds.spec_axis] += my_off
- end[self.ds.spec_axis] += my_off
- mylog.debug("Reading from " + str(start) + str(end))
slices = [slice(start[i],end[i]) for i in range(3)]
- if self.ds.reversed:
- new_start = self.ds.dims[self.ds.spec_axis]-1-start[self.ds.spec_axis]
- new_end = max(self.ds.dims[self.ds.spec_axis]-1-end[self.ds.spec_axis],0)
- slices[self.ds.spec_axis] = slice(new_start,new_end,-1)
if self.ds.dimensionality == 2:
nx, ny = g.ActiveDimensions[:2]
nz = 1
@@ -115,13 +92,6 @@
data = ds.data[idx,slices[2],slices[1],slices[0]].transpose()
else:
data = ds.data[slices[2],slices[1],slices[0]].transpose()
- if self.line_db is not None:
- nz1 = data.shape[self.ds.spec_axis]
- nz2 = g.ActiveDimensions[self.ds.spec_axis]
- if nz1 != nz2:
- old_data = data.copy()
- data = np.zeros(g.ActiveDimensions)
- data[:,:,nz2-nz1:] = old_data
if fname in self.ds.nan_mask:
data[np.isnan(data)] = self.ds.nan_mask[fname]
elif "all" in self.ds.nan_mask:
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/frontends/fits/misc.py
--- a/yt/frontends/fits/misc.py
+++ b/yt/frontends/fits/misc.py
@@ -17,6 +17,8 @@
from yt.utilities.on_demand_imports import _astropy
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
import os
@@ -68,6 +70,70 @@
validators = [ValidateSpatial()],
display_name="Counts (%s-%s keV)" % (emin, emax))
+def create_spectral_slabs(filename, slab_centers, slab_width,
+ **kwargs):
+ r"""
+ Given a dictionary of spectral slab centers and a width in
+ spectral units, extract data from a spectral cube at these slab
+ centers and return a `FITSDataset` instance containing the different
+ slabs as separate yt fields. Useful for extracting individual
+ lines from a spectral cube and separating them out as different fields.
+
+ Requires the SpectralCube (http://spectral-cube.readthedocs.org)
+ library.
+
+ All keyword arguments will be passed on to the `FITSDataset` constructor.
+
+ Parameters
+ ----------
+ filename : string
+ The spectral cube FITS file to extract the data from.
+ slab_centers : dict of (float, string) tuples or YTQuantities
+ The centers of the slabs, where the keys are the names
+ of the new fields and the values are (float, string) tuples or
+ YTQuantities, specifying a value for each center and its unit.
+ slab_width : YTQuantity or (float, string) tuple
+ The width of the slab along the spectral axis.
+
+ Examples
+ --------
+ >>> slab_centers = {'13CN': (218.03117, 'GHz'),
+ ... 'CH3CH2CHO': (218.284256, 'GHz'),
+ ... 'CH3NH2': (218.40956, 'GHz')}
+ >>> slab_width = (0.05, "GHz")
+ >>> ds = create_spectral_slabs("intensity_cube.fits",
+ ... slab_centers, slab_width,
+ ... nan_mask=0.0)
+ """
+ from spectral_cube import SpectralCube
+ from yt.frontends.fits.api import FITSDataset
+ cube = SpectralCube.read(filename)
+ if not isinstance(slab_width, YTQuantity):
+ slab_width = YTQuantity(slab_width[0], slab_width[1])
+ slab_data = {}
+ field_units = cube.header.get("bunit", "dimensionless")
+ for k, v in slab_centers.items():
+ if not isinstance(v, YTQuantity):
+ slab_center = YTQuantity(v[0], v[1])
+ else:
+ slab_center = v
+ mylog.info("Adding slab field %s at %g %s" %
+ (k, slab_center.v, slab_center.units))
+ slab_lo = (slab_center-0.5*slab_width).to_astropy()
+ slab_hi = (slab_center+0.5*slab_width).to_astropy()
+ subcube = cube.spectral_slab(slab_lo, slab_hi)
+ slab_data[k] = YTArray(subcube.filled_data[:,:,:], field_units)
+ width = subcube.header["naxis3"]*cube.header["cdelt3"]
+ w = subcube.wcs.copy()
+ w.wcs.crpix[-1] = 0.5
+ w.wcs.crval[-1] = -0.5*width
+ fid = FITSImageBuffer(slab_data, wcs=w)
+ for hdu in fid:
+ hdu.header.pop("RESTFREQ", None)
+ hdu.header.pop("RESTFRQ", None)
+ ds = FITSDataset(fid, **kwargs)
+ return ds
+
def ds9_region(ds, reg, obj=None, field_parameters=None):
r"""
Create a data container from a ds9 region file. Requires the pyregion
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/frontends/gadget/api.py
--- a/yt/frontends/gadget/api.py
+++ b/yt/frontends/gadget/api.py
@@ -7,7 +7,7 @@
"""
#-----------------------------------------------------------------------------
-# Copyright (c) 2014, yt Development Team.
+# Copyright (c) 2014-2015, yt Development Team.
#
# Distributed under the terms of the Modified BSD License.
#
@@ -23,4 +23,7 @@
IOHandlerGadgetBinary, \
IOHandlerGadgetHDF5
+from .simulation_handling import \
+ GadgetSimulation
+
from . import tests
diff -r 696a6134f3f070d4b55fcfc4c2554db563aab11e -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 yt/frontends/gadget/data_structures.py
--- a/yt/frontends/gadget/data_structures.py
+++ b/yt/frontends/gadget/data_structures.py
@@ -18,6 +18,7 @@
import h5py
import numpy as np
import stat
+import struct
import os
import types
@@ -242,10 +243,59 @@
self.mass_unit = self.quan(mass_unit[0], mass_unit[1])
self.time_unit = self.length_unit / self.velocity_unit
+ @staticmethod
+ def _validate_header(filename):
+ '''
+ This method automatically detects whether the Gadget file is big/little endian
+ and is not corrupt/invalid using the first 4 bytes in the file. It returns a
+ tuple of (Valid, endianswap) where Valid is a boolean that is true if the file
+ is a Gadget binary file, and endianswap is the endianness character '>' or '<'.
+ '''
+ try:
+ f = open(filename,'rb')
+ except IOError:
+ try:
+ f = open(filename+".0")
+ except IOError:
+ return False, 1
+
+ # First int32 is 256 for a Gadget2 binary file with SnapFormat=1,
+ # 8 for a Gadget2 binary file with SnapFormat=2 file,
+ # or the byte swapped equivalents (65536 and 134217728).
+ # The int32 following the header (first 4+256 bytes) must equal this
+ # number.
+ (rhead,) = struct.unpack('<I',f.read(4))
+ # Use value to check endianess
+ if rhead == 256:
+ endianswap = '<'
+ elif rhead == 65536:
+ endianswap = '>'
+ # Disabled for now (does any one still use SnapFormat=2?)
+ # If so, alternate read would be needed based on header.
+ # elif rhead == 8:
+ # return True, '<'
+ # elif rhead == 134217728:
+ # return True, '>'
+ else:
+ f.close()
+ return False, 1
+ # Read in particle number from header
+ np0 = sum(struct.unpack(endianswap+'IIIIII',f.read(6*4)))
+ # Read in size of position block. It should be 4 bytes per float,
+ # with 3 coordinates (x,y,z) per particle. (12 bytes per particle)
+ f.seek(4+256+4,0)
+ np1 = struct.unpack(endianswap+'I',f.read(4))[0]/(4*3)
+ f.close()
+ # Compare
+ if np0 == np1:
+ return True, endianswap
+ else:
+ return False, 1
+
@classmethod
def _is_valid(self, *args, **kwargs):
- # We do not allow load() of these files.
- return False
+ # First 4 bytes used to check load
+ return GadgetDataset._validate_header(args[0])[0]
class GadgetHDF5Dataset(GadgetDataset):
_file_class = ParticleFile
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/yt_analysis/yt/commits/45467416f055/
Changeset: 45467416f055
Branch: yt
User: jzuhone
Date: 2015-06-02 00:36:26+00:00
Summary: Merge
Affected #: 16 files
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe distribute_setup.py
--- a/distribute_setup.py
+++ /dev/null
@@ -1,541 +0,0 @@
-#!python
-"""Bootstrap distribute installation
-
-If you want to use setuptools in your package's setup.py, just include this
-file in the same directory with it, and add this to the top of your setup.py::
-
- from distribute_setup import use_setuptools
- use_setuptools()
-
-If you want to require a specific version of setuptools, set a download
-mirror, or use an alternate download directory, you can do so by supplying
-the appropriate options to ``use_setuptools()``.
-
-This file can also be run as a script to install or upgrade setuptools.
-"""
-import os
-import shutil
-import sys
-import time
-import fnmatch
-import tempfile
-import tarfile
-import optparse
-
-from distutils import log
-
-try:
- from site import USER_SITE
-except ImportError:
- USER_SITE = None
-
-try:
- import subprocess
-
- def _python_cmd(*args):
- args = (sys.executable,) + args
- return subprocess.call(args) == 0
-
-except ImportError:
- # will be used for python 2.3
- def _python_cmd(*args):
- args = (sys.executable,) + args
- # quoting arguments if windows
- if sys.platform == 'win32':
- def quote(arg):
- if ' ' in arg:
- return '"%s"' % arg
- return arg
- args = [quote(arg) for arg in args]
- return os.spawnl(os.P_WAIT, sys.executable, *args) == 0
-
-DEFAULT_VERSION = "0.6.32"
-DEFAULT_URL = "http://pypi.python.org/packages/source/d/distribute/"
-SETUPTOOLS_FAKED_VERSION = "0.6c11"
-
-SETUPTOOLS_PKG_INFO = """\
-Metadata-Version: 1.0
-Name: setuptools
-Version: %s
-Summary: xxxx
-Home-page: xxx
-Author: xxx
-Author-email: xxx
-License: xxx
-Description: xxx
-""" % SETUPTOOLS_FAKED_VERSION
-
-
-def _install(tarball, install_args=()):
- # extracting the tarball
- tmpdir = tempfile.mkdtemp()
- log.warn('Extracting in %s', tmpdir)
- old_wd = os.getcwd()
- try:
- os.chdir(tmpdir)
- tar = tarfile.open(tarball)
- _extractall(tar)
- tar.close()
-
- # going in the directory
- subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0])
- os.chdir(subdir)
- log.warn('Now working in %s', subdir)
-
- # installing
- log.warn('Installing Distribute')
- if not _python_cmd('setup.py', 'install', *install_args):
- log.warn('Something went wrong during the installation.')
- log.warn('See the error message above.')
- # exitcode will be 2
- return 2
- finally:
- os.chdir(old_wd)
- shutil.rmtree(tmpdir)
-
-
-def _build_egg(egg, tarball, to_dir):
- # extracting the tarball
- tmpdir = tempfile.mkdtemp()
- log.warn('Extracting in %s', tmpdir)
- old_wd = os.getcwd()
- try:
- os.chdir(tmpdir)
- tar = tarfile.open(tarball)
- _extractall(tar)
- tar.close()
-
- # going in the directory
- subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0])
- os.chdir(subdir)
- log.warn('Now working in %s', subdir)
-
- # building an egg
- log.warn('Building a Distribute egg in %s', to_dir)
- _python_cmd('setup.py', '-q', 'bdist_egg', '--dist-dir', to_dir)
-
- finally:
- os.chdir(old_wd)
- shutil.rmtree(tmpdir)
- # returning the result
- log.warn(egg)
- if not os.path.exists(egg):
- raise IOError('Could not build the egg.')
-
-
-def _do_download(version, download_base, to_dir, download_delay):
- egg = os.path.join(to_dir, 'distribute-%s-py%d.%d.egg'
- % (version, sys.version_info[0], sys.version_info[1]))
- if not os.path.exists(egg):
- tarball = download_setuptools(version, download_base,
- to_dir, download_delay)
- _build_egg(egg, tarball, to_dir)
- sys.path.insert(0, egg)
- import setuptools
- setuptools.bootstrap_install_from = egg
-
-
-def use_setuptools(version=DEFAULT_VERSION, download_base=DEFAULT_URL,
- to_dir=os.curdir, download_delay=15, no_fake=True):
- # making sure we use the absolute path
- to_dir = os.path.abspath(to_dir)
- was_imported = 'pkg_resources' in sys.modules or \
- 'setuptools' in sys.modules
- try:
- try:
- import pkg_resources
- if not hasattr(pkg_resources, '_distribute'):
- if not no_fake:
- _fake_setuptools()
- raise ImportError
- except ImportError:
- return _do_download(version, download_base, to_dir, download_delay)
- try:
- pkg_resources.require("distribute>=" + version)
- return
- except pkg_resources.VersionConflict:
- e = sys.exc_info()[1]
- if was_imported:
- sys.stderr.write(
- "The required version of distribute (>=%s) is not available,\n"
- "and can't be installed while this script is running. Please\n"
- "install a more recent version first, using\n"
- "'easy_install -U distribute'."
- "\n\n(Currently using %r)\n" % (version, e.args[0]))
- sys.exit(2)
- else:
- del pkg_resources, sys.modules['pkg_resources'] # reload ok
- return _do_download(version, download_base, to_dir,
- download_delay)
- except pkg_resources.DistributionNotFound:
- return _do_download(version, download_base, to_dir,
- download_delay)
- finally:
- if not no_fake:
- _create_fake_setuptools_pkg_info(to_dir)
-
-
-def download_setuptools(version=DEFAULT_VERSION, download_base=DEFAULT_URL,
- to_dir=os.curdir, delay=15):
- """Download distribute from a specified location and return its filename
-
- `version` should be a valid distribute version number that is available
- as an egg for download under the `download_base` URL (which should end
- with a '/'). `to_dir` is the directory where the egg will be downloaded.
- `delay` is the number of seconds to pause before an actual download
- attempt.
- """
- # making sure we use the absolute path
- to_dir = os.path.abspath(to_dir)
- try:
- from urllib.request import urlopen
- except ImportError:
- from urllib2 import urlopen
- tgz_name = "distribute-%s.tar.gz" % version
- url = download_base + tgz_name
- saveto = os.path.join(to_dir, tgz_name)
- src = dst = None
- if not os.path.exists(saveto): # Avoid repeated downloads
- try:
- log.warn("Downloading %s", url)
- src = urlopen(url)
- # Read/write all in one block, so we don't create a corrupt file
- # if the download is interrupted.
- data = src.read()
- dst = open(saveto, "wb")
- dst.write(data)
- finally:
- if src:
- src.close()
- if dst:
- dst.close()
- return os.path.realpath(saveto)
-
-
-def _no_sandbox(function):
- def __no_sandbox(*args, **kw):
- try:
- from setuptools.sandbox import DirectorySandbox
- if not hasattr(DirectorySandbox, '_old'):
- def violation(*args):
- pass
- DirectorySandbox._old = DirectorySandbox._violation
- DirectorySandbox._violation = violation
- patched = True
- else:
- patched = False
- except ImportError:
- patched = False
-
- try:
- return function(*args, **kw)
- finally:
- if patched:
- DirectorySandbox._violation = DirectorySandbox._old
- del DirectorySandbox._old
-
- return __no_sandbox
-
-
-def _patch_file(path, content):
- """Will backup the file then patch it"""
- existing_content = open(path).read()
- if existing_content == content:
- # already patched
- log.warn('Already patched.')
- return False
- log.warn('Patching...')
- _rename_path(path)
- f = open(path, 'w')
- try:
- f.write(content)
- finally:
- f.close()
- return True
-
-_patch_file = _no_sandbox(_patch_file)
-
-
-def _same_content(path, content):
- return open(path).read() == content
-
-
-def _rename_path(path):
- new_name = path + '.OLD.%s' % time.time()
- log.warn('Renaming %s to %s', path, new_name)
- os.rename(path, new_name)
- return new_name
-
-
-def _remove_flat_installation(placeholder):
- if not os.path.isdir(placeholder):
- log.warn('Unknown installation at %s', placeholder)
- return False
- found = False
- for file in os.listdir(placeholder):
- if fnmatch.fnmatch(file, 'setuptools*.egg-info'):
- found = True
- break
- if not found:
- log.warn('Could not locate setuptools*.egg-info')
- return
-
- log.warn('Moving elements out of the way...')
- pkg_info = os.path.join(placeholder, file)
- if os.path.isdir(pkg_info):
- patched = _patch_egg_dir(pkg_info)
- else:
- patched = _patch_file(pkg_info, SETUPTOOLS_PKG_INFO)
-
- if not patched:
- log.warn('%s already patched.', pkg_info)
- return False
- # now let's move the files out of the way
- for element in ('setuptools', 'pkg_resources.py', 'site.py'):
- element = os.path.join(placeholder, element)
- if os.path.exists(element):
- _rename_path(element)
- else:
- log.warn('Could not find the %s element of the '
- 'Setuptools distribution', element)
- return True
-
-_remove_flat_installation = _no_sandbox(_remove_flat_installation)
-
-
-def _after_install(dist):
- log.warn('After install bootstrap.')
- placeholder = dist.get_command_obj('install').install_purelib
- _create_fake_setuptools_pkg_info(placeholder)
-
-
-def _create_fake_setuptools_pkg_info(placeholder):
- if not placeholder or not os.path.exists(placeholder):
- log.warn('Could not find the install location')
- return
- pyver = '%s.%s' % (sys.version_info[0], sys.version_info[1])
- setuptools_file = 'setuptools-%s-py%s.egg-info' % \
- (SETUPTOOLS_FAKED_VERSION, pyver)
- pkg_info = os.path.join(placeholder, setuptools_file)
- if os.path.exists(pkg_info):
- log.warn('%s already exists', pkg_info)
- return
-
- log.warn('Creating %s', pkg_info)
- try:
- f = open(pkg_info, 'w')
- except EnvironmentError:
- log.warn("Don't have permissions to write %s, skipping", pkg_info)
- return
- try:
- f.write(SETUPTOOLS_PKG_INFO)
- finally:
- f.close()
-
- pth_file = os.path.join(placeholder, 'setuptools.pth')
- log.warn('Creating %s', pth_file)
- f = open(pth_file, 'w')
- try:
- f.write(os.path.join(os.curdir, setuptools_file))
- finally:
- f.close()
-
-_create_fake_setuptools_pkg_info = _no_sandbox(
- _create_fake_setuptools_pkg_info
-)
-
-
-def _patch_egg_dir(path):
- # let's check if it's already patched
- pkg_info = os.path.join(path, 'EGG-INFO', 'PKG-INFO')
- if os.path.exists(pkg_info):
- if _same_content(pkg_info, SETUPTOOLS_PKG_INFO):
- log.warn('%s already patched.', pkg_info)
- return False
- _rename_path(path)
- os.mkdir(path)
- os.mkdir(os.path.join(path, 'EGG-INFO'))
- pkg_info = os.path.join(path, 'EGG-INFO', 'PKG-INFO')
- f = open(pkg_info, 'w')
- try:
- f.write(SETUPTOOLS_PKG_INFO)
- finally:
- f.close()
- return True
-
-_patch_egg_dir = _no_sandbox(_patch_egg_dir)
-
-
-def _before_install():
- log.warn('Before install bootstrap.')
- _fake_setuptools()
-
-
-def _under_prefix(location):
- if 'install' not in sys.argv:
- return True
- args = sys.argv[sys.argv.index('install') + 1:]
- for index, arg in enumerate(args):
- for option in ('--root', '--prefix'):
- if arg.startswith('%s=' % option):
- top_dir = arg.split('root=')[-1]
- return location.startswith(top_dir)
- elif arg == option:
- if len(args) > index:
- top_dir = args[index + 1]
- return location.startswith(top_dir)
- if arg == '--user' and USER_SITE is not None:
- return location.startswith(USER_SITE)
- return True
-
-
-def _fake_setuptools():
- log.warn('Scanning installed packages')
- try:
- import pkg_resources
- except ImportError:
- # we're cool
- log.warn('Setuptools or Distribute does not seem to be installed.')
- return
- ws = pkg_resources.working_set
- try:
- setuptools_dist = ws.find(
- pkg_resources.Requirement.parse('setuptools', replacement=False)
- )
- except TypeError:
- # old distribute API
- setuptools_dist = ws.find(
- pkg_resources.Requirement.parse('setuptools')
- )
-
- if setuptools_dist is None:
- log.warn('No setuptools distribution found')
- return
- # detecting if it was already faked
- setuptools_location = setuptools_dist.location
- log.warn('Setuptools installation detected at %s', setuptools_location)
-
- # if --root or --preix was provided, and if
- # setuptools is not located in them, we don't patch it
- if not _under_prefix(setuptools_location):
- log.warn('Not patching, --root or --prefix is installing Distribute'
- ' in another location')
- return
-
- # let's see if its an egg
- if not setuptools_location.endswith('.egg'):
- log.warn('Non-egg installation')
- res = _remove_flat_installation(setuptools_location)
- if not res:
- return
- else:
- log.warn('Egg installation')
- pkg_info = os.path.join(setuptools_location, 'EGG-INFO', 'PKG-INFO')
- if (os.path.exists(pkg_info) and
- _same_content(pkg_info, SETUPTOOLS_PKG_INFO)):
- log.warn('Already patched.')
- return
- log.warn('Patching...')
- # let's create a fake egg replacing setuptools one
- res = _patch_egg_dir(setuptools_location)
- if not res:
- return
- log.warn('Patching complete.')
- _relaunch()
-
-
-def _relaunch():
- log.warn('Relaunching...')
- # we have to relaunch the process
- # pip marker to avoid a relaunch bug
- _cmd1 = ['-c', 'install', '--single-version-externally-managed']
- _cmd2 = ['-c', 'install', '--record']
- if sys.argv[:3] == _cmd1 or sys.argv[:3] == _cmd2:
- sys.argv[0] = 'setup.py'
- args = [sys.executable] + sys.argv
- sys.exit(subprocess.call(args))
-
-
-def _extractall(self, path=".", members=None):
- """Extract all members from the archive to the current working
- directory and set owner, modification time and permissions on
- directories afterwards. `path' specifies a different directory
- to extract to. `members' is optional and must be a subset of the
- list returned by getmembers().
- """
- import copy
- import operator
- from tarfile import ExtractError
- directories = []
-
- if members is None:
- members = self
-
- for tarinfo in members:
- if tarinfo.isdir():
- # Extract directories with a safe mode.
- directories.append(tarinfo)
- tarinfo = copy.copy(tarinfo)
- tarinfo.mode = 448 # decimal for oct 0700
- self.extract(tarinfo, path)
-
- # Reverse sort directories.
- if sys.version_info < (2, 4):
- def sorter(dir1, dir2):
- return cmp(dir1.name, dir2.name)
- directories.sort(sorter)
- directories.reverse()
- else:
- directories.sort(key=operator.attrgetter('name'), reverse=True)
-
- # Set correct owner, mtime and filemode on directories.
- for tarinfo in directories:
- dirpath = os.path.join(path, tarinfo.name)
- try:
- self.chown(tarinfo, dirpath)
- self.utime(tarinfo, dirpath)
- self.chmod(tarinfo, dirpath)
- except ExtractError:
- e = sys.exc_info()[1]
- if self.errorlevel > 1:
- raise
- else:
- self._dbg(1, "tarfile: %s" % e)
-
-
-def _build_install_args(options):
- """
- Build the arguments to 'python setup.py install' on the distribute package
- """
- install_args = []
- if options.user_install:
- if sys.version_info < (2, 6):
- log.warn("--user requires Python 2.6 or later")
- raise SystemExit(1)
- install_args.append('--user')
- return install_args
-
-def _parse_args():
- """
- Parse the command line for options
- """
- parser = optparse.OptionParser()
- parser.add_option(
- '--user', dest='user_install', action='store_true', default=False,
- help='install in user site package (requires Python 2.6 or later)')
- parser.add_option(
- '--download-base', dest='download_base', metavar="URL",
- default=DEFAULT_URL,
- help='alternative URL from where to download the distribute package')
- options, args = parser.parse_args()
- # positional arguments are ignored
- return options
-
-def main(version=DEFAULT_VERSION):
- """Install or upgrade setuptools and EasyInstall"""
- options = _parse_args()
- tarball = download_setuptools(download_base=options.download_base)
- return _install(tarball, _build_install_args(options))
-
-if __name__ == '__main__':
- sys.exit(main())
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe doc/install_script.sh
--- a/doc/install_script.sh
+++ b/doc/install_script.sh
@@ -1,18 +1,14 @@
#
# Hi there! Welcome to the yt installation script.
#
+# First things first, if you experience problems, please visit the Help
+# section at http://yt-project.org.
+#
# This script is designed to create a fully isolated Python installation
# with the dependencies you need to run yt.
#
-# There are a few options, but you only need to set *one* of them. And
-# that's the next one, DEST_DIR. But, if you want to use an existing HDF5
-# installation you can set HDF5_DIR, or if you want to use some other
-# subversion checkout of yt, you can set YT_DIR, too. (It'll already
-# check the current directory and one up.
-#
-# If you experience problems, please visit the Help section at
-# http://yt-project.org.
-#
+# There are a few options, but you only need to set *one* of them, which is
+# the next one, DEST_DIR:
DEST_SUFFIX="yt-`uname -m`"
DEST_DIR="`pwd`/${DEST_SUFFIX/ /}" # Installation location
@@ -23,16 +19,25 @@
DEST_DIR=${YT_DEST}
fi
+# What follows are some other options that you may or may not need to change.
+
# Here's where you put the HDF5 path if you like; otherwise it'll download it
# and install it on its own
#HDF5_DIR=
+# If you've got yt some other place, set this to point to it. The script will
+# already check the current directory and the one above it in the tree.
+YT_DIR=""
+
# If you need to supply arguments to the NumPy or SciPy build, supply them here
# This one turns on gfortran manually:
#NUMPY_ARGS="--fcompiler=gnu95"
# If you absolutely can't get the fortran to work, try this:
#NUMPY_ARGS="--fcompiler=fake"
+INST_PY3=0 # Install Python 3 along with Python 2. If this is turned
+ # on, all Python packages (including yt) will be installed
+ # in Python 3 (except Mercurial, which requires Python 2).
INST_HG=1 # Install Mercurial or not? If hg is not already
# installed, yt cannot be installed.
INST_ZLIB=1 # On some systems (Kraken) matplotlib has issues with
@@ -50,9 +55,6 @@
INST_ROCKSTAR=0 # Install the Rockstar halo finder?
INST_SCIPY=0 # Install scipy?
-# If you've got yt some other place, set this to point to it.
-YT_DIR=""
-
# If you need to pass anything to matplotlib, do so here.
MPL_SUPP_LDFLAGS=""
MPL_SUPP_CFLAGS=""
@@ -111,6 +113,7 @@
echo INST_SQLITE3=${INST_SQLITE3} >> ${CONFIG_FILE}
echo INST_PYX=${INST_PYX} >> ${CONFIG_FILE}
echo INST_0MQ=${INST_0MQ} >> ${CONFIG_FILE}
+ echo INST_PY3=${INST_PY3} >> ${CONFIG_FILE}
echo INST_ROCKSTAR=${INST_ROCKSTAR} >> ${CONFIG_FILE}
echo INST_SCIPY=${INST_SCIPY} >> ${CONFIG_FILE}
echo YT_DIR=${YT_DIR} >> ${CONFIG_FILE}
@@ -415,6 +418,10 @@
get_willwont ${INST_SQLITE3}
echo "be installing SQLite3"
+printf "%-15s = %s so I " "INST_PY3" "${INST_PY3}"
+get_willwont ${INST_PY3}
+echo "be installing Python 3"
+
printf "%-15s = %s so I " "INST_HG" "${INST_HG}"
get_willwont ${INST_HG}
echo "be installing Mercurial"
@@ -487,6 +494,13 @@
exit 1
}
+if [ $INST_PY3 -eq 1 ]
+then
+ PYTHON_EXEC='python3.4'
+else
+ PYTHON_EXEC='python2.7'
+fi
+
function do_setup_py
{
[ -e $1/done ] && return
@@ -501,21 +515,27 @@
[ ! -e $LIB/extracted ] && tar xfz $LIB.tar.gz
touch $LIB/extracted
BUILD_ARGS=""
+ if [[ $LIB =~ .*mercurial.* ]]
+ then
+ PYEXE="python2.7"
+ else
+ PYEXE=${PYTHON_EXEC}
+ fi
case $LIB in
*h5py*)
pushd $LIB &> /dev/null
- ( ${DEST_DIR}/bin/python2.7 setup.py configure --hdf5=${HDF5_DIR} 2>&1 ) 1>> ${LOG_FILE} || do_exit
+ ( ${DEST_DIR}/bin/${PYTHON_EXEC} setup.py configure --hdf5=${HDF5_DIR} 2>&1 ) 1>> ${LOG_FILE} || do_exit
popd &> /dev/null
;;
*numpy*)
- if [ -e ${DEST_DIR}/lib/python2.7/site-packages/numpy/__init__.py ]
+ if [ -e ${DEST_DIR}/lib/${PYTHON_EXEC}/site-packages/numpy/__init__.py ]
then
- VER=$(${DEST_DIR}/bin/python -c 'from distutils.version import StrictVersion as SV; \
+ VER=$(${DEST_DIR}/bin/${PYTHON_EXEC} -c 'from distutils.version import StrictVersion as SV; \
import numpy; print SV(numpy.__version__) < SV("1.8.0")')
if [ $VER == "True" ]
then
echo "Removing previous NumPy instance (see issue #889)"
- rm -rf ${DEST_DIR}/lib/python2.7/site-packages/{numpy*,*.pth}
+ rm -rf ${DEST_DIR}/lib/${PYTHON_EXEC}/site-packages/{numpy*,*.pth}
fi
fi
;;
@@ -523,8 +543,8 @@
;;
esac
cd $LIB
- ( ${DEST_DIR}/bin/python2.7 setup.py build ${BUILD_ARGS} $* 2>&1 ) 1>> ${LOG_FILE} || do_exit
- ( ${DEST_DIR}/bin/python2.7 setup.py install 2>&1 ) 1>> ${LOG_FILE} || do_exit
+ ( ${DEST_DIR}/bin/${PYEXE} setup.py build ${BUILD_ARGS} $* 2>&1 ) 1>> ${LOG_FILE} || do_exit
+ ( ${DEST_DIR}/bin/${PYEXE} setup.py install 2>&1 ) 1>> ${LOG_FILE} || do_exit
touch done
cd ..
}
@@ -592,14 +612,15 @@
# Set paths to what they should be when yt is activated.
export PATH=${DEST_DIR}/bin:$PATH
export LD_LIBRARY_PATH=${DEST_DIR}/lib:$LD_LIBRARY_PATH
-export PYTHONPATH=${DEST_DIR}/lib/python2.7/site-packages
+export PYTHONPATH=${DEST_DIR}/lib/${PYTHON_EXEC}/site-packages
mkdir -p ${DEST_DIR}/src
cd ${DEST_DIR}/src
+PYTHON2='Python-2.7.9'
+PYTHON3='Python-3.4.3'
CYTHON='Cython-0.22'
PYX='PyX-0.12.1'
-PYTHON='Python-2.7.9'
BZLIB='bzip2-1.0.6'
FREETYPE_VER='freetype-2.4.12'
H5PY='h5py-2.5.0'
@@ -620,11 +641,13 @@
TORNADO='tornado-4.0.2'
ZEROMQ='zeromq-4.0.5'
ZLIB='zlib-1.2.8'
+SETUPTOOLS='setuptools-16.0'
# Now we dump all our SHA512 files out.
echo '856220fa579e272ac38dcef091760f527431ff3b98df9af6e68416fcf77d9659ac5abe5c7dee41331f359614637a4ff452033085335ee499830ed126ab584267 Cython-0.22.tar.gz' > Cython-0.22.tar.gz.sha512
echo '4941f5aa21aff3743546495fb073c10d2657ff42b2aff401903498638093d0e31e344cce778980f28a7170c6d29eab72ac074277b9d4088376e8692dc71e55c1 PyX-0.12.1.tar.gz' > PyX-0.12.1.tar.gz.sha512
echo 'a42f28ed8e49f04cf89e2ea7434c5ecbc264e7188dcb79ab97f745adf664dd9ab57f9a913543731635f90859536244ac37dca9adf0fc2aa1b215ba884839d160 Python-2.7.9.tgz' > Python-2.7.9.tgz.sha512
+echo '609cc82586fabecb25f25ecb410f2938e01d21cde85dd3f8824fe55c6edde9ecf3b7609195473d3fa05a16b9b121464f5414db1a0187103b78ea6edfa71684a7 Python-3.4.3.tgz' > Python-3.4.3.tgz.sha512
echo '276bd9c061ec9a27d478b33078a86f93164ee2da72210e12e2c9da71dcffeb64767e4460b93f257302b09328eda8655e93c4b9ae85e74472869afbeae35ca71e blas.tar.gz' > blas.tar.gz.sha512
echo '00ace5438cfa0c577e5f578d8a808613187eff5217c35164ffe044fbafdfec9e98f4192c02a7d67e01e5a5ccced630583ad1003c37697219b0f147343a3fdd12 bzip2-1.0.6.tar.gz' > bzip2-1.0.6.tar.gz.sha512
echo 'a296dfcaef7e853e58eed4e24b37c4fa29cfc6ac688def048480f4bb384b9e37ca447faf96eec7b378fd764ba291713f03ac464581d62275e28eb2ec99110ab6 reason-js-20120623.zip' > reason-js-20120623.zip.sha512
@@ -646,6 +669,7 @@
echo '93591068dc63af8d50a7925d528bc0cccdd705232c529b6162619fe28dddaf115e8a460b1842877d35160bd7ed480c1bd0bdbec57d1f359085bd1814e0c1c242 tornado-4.0.2.tar.gz' > tornado-4.0.2.tar.gz.sha512
echo '0d928ed688ed940d460fa8f8d574a9819dccc4e030d735a8c7db71b59287ee50fa741a08249e356c78356b03c2174f2f2699f05aa7dc3d380ed47d8d7bab5408 zeromq-4.0.5.tar.gz' > zeromq-4.0.5.tar.gz.sha512
echo 'ece209d4c7ec0cb58ede791444dc754e0d10811cbbdebe3df61c0fd9f9f9867c1c3ccd5f1827f847c005e24eef34fb5bf87b5d3f894d75da04f1797538290e4a zlib-1.2.8.tar.gz' > zlib-1.2.8.tar.gz.sha512
+echo '38a89aad89dc9aa682dbfbca623e2f69511f5e20d4a3526c01aabbc7e93ae78f20aac566676b431e111540b41540a1c4f644ce4174e7ecf052318612075e02dc setuptools-16.0.tar.gz' > setuptools-16.0.tar.gz.sha512
# Individual processes
[ -z "$HDF5_DIR" ] && get_ytproject $HDF5.tar.gz
[ $INST_ZLIB -eq 1 ] && get_ytproject $ZLIB.tar.gz
@@ -660,10 +684,11 @@
[ $INST_SCIPY -eq 1 ] && get_ytproject $SCIPY.tar.gz
[ $INST_SCIPY -eq 1 ] && get_ytproject blas.tar.gz
[ $INST_SCIPY -eq 1 ] && get_ytproject $LAPACK.tar.gz
-get_ytproject $PYTHON.tgz
+[ $INST_HG -eq 1 ] && get_ytproject $MERCURIAL.tar.gz
+[ $INST_PY3 -eq 1 ] && get_ytproject $PYTHON3.tgz
+get_ytproject $PYTHON2.tgz
get_ytproject $NUMPY.tar.gz
get_ytproject $MATPLOTLIB.tar.gz
-get_ytproject $MERCURIAL.tar.gz
get_ytproject $IPYTHON.tar.gz
get_ytproject $H5PY.tar.gz
get_ytproject $CYTHON.tar.gz
@@ -671,6 +696,7 @@
get_ytproject $NOSE.tar.gz
get_ytproject $PYTHON_HGLIB.tar.gz
get_ytproject $SYMPY.tar.gz
+get_ytproject $SETUPTOOLS.tar.gz
if [ $INST_BZLIB -eq 1 ]
then
if [ ! -e $BZLIB/done ]
@@ -787,11 +813,11 @@
fi
fi
-if [ ! -e $PYTHON/done ]
+if [ ! -e $PYTHON2/done ]
then
- echo "Installing Python. This may take a while, but don't worry. yt loves you."
- [ ! -e $PYTHON ] && tar xfz $PYTHON.tgz
- cd $PYTHON
+ echo "Installing Python 2. This may take a while, but don't worry. yt loves you."
+ [ ! -e $PYTHON2 ] && tar xfz $PYTHON2.tgz
+ cd $PYTHON2
( ./configure --prefix=${DEST_DIR}/ ${PYCONF_ARGS} 2>&1 ) 1>> ${LOG_FILE} || do_exit
( make ${MAKE_PROCS} 2>&1 ) 1>> ${LOG_FILE} || do_exit
@@ -802,7 +828,30 @@
cd ..
fi
-export PYTHONPATH=${DEST_DIR}/lib/python2.7/site-packages/
+if [ $INST_PY3 -eq 1 ]
+then
+ if [ ! -e $PYTHON3/done ]
+ then
+ echo "Installing Python 3. Because two Pythons are better than one."
+ [ ! -e $PYTHON3 ] && tar xfz $PYTHON3.tgz
+ cd $PYTHON3
+ ( ./configure --prefix=${DEST_DIR}/ ${PYCONF_ARGS} 2>&1 ) 1>> ${LOG_FILE} || do_exit
+
+ ( make ${MAKE_PROCS} 2>&1 ) 1>> ${LOG_FILE} || do_exit
+ ( make install 2>&1 ) 1>> ${LOG_FILE} || do_exit
+ ( ln -sf ${DEST_DIR}/bin/python3.4 ${DEST_DIR}/bin/pyyt 2>&1 ) 1>> ${LOG_FILE}
+ ( ln -sf ${DEST_DIR}/bin/python3.4 ${DEST_DIR}/bin/python 2>&1 ) 1>> ${LOG_FILE}
+ ( ln -sf ${DEST_DIR}/bin/python3-config ${DEST_DIR}/bin/python-config 2>&1 ) 1>> ${LOG_FILE}
+ ( make clean 2>&1) 1>> ${LOG_FILE} || do_exit
+ touch done
+ cd ..
+ fi
+fi
+
+export PYTHONPATH=${DEST_DIR}/lib/${PYTHON_EXEC}/site-packages/
+
+# Install setuptools
+do_setup_py $SETUPTOOLS
if [ $INST_HG -eq 1 ]
then
@@ -847,12 +896,10 @@
# This fixes problems with gfortran linking.
unset LDFLAGS
-
-echo "Installing distribute"
-( ${DEST_DIR}/bin/python2.7 ${YT_DIR}/distribute_setup.py 2>&1 ) 1>> ${LOG_FILE} || do_exit
-
+
echo "Installing pip"
-( ${DEST_DIR}/bin/easy_install-2.7 pip 2>&1 ) 1>> ${LOG_FILE} || do_exit
+( ${GETFILE} https://bootstrap.pypa.io/get-pip.py 2>&1 ) 1>> ${LOG_FILE} || do_exit
+( ${DEST_DIR}/bin/${PYTHON_EXEC} get-pip.py 2>&1 ) 1>> ${LOG_FILE} || do_exit
if [ $INST_SCIPY -eq 0 ]
then
@@ -986,13 +1033,14 @@
echo "Installing yt"
[ $INST_PNG -eq 1 ] && echo $PNG_DIR > png.cfg
-( export PATH=$DEST_DIR/bin:$PATH ; ${DEST_DIR}/bin/python2.7 setup.py develop 2>&1 ) 1>> ${LOG_FILE} || do_exit
+( export PATH=$DEST_DIR/bin:$PATH ; ${DEST_DIR}/bin/${PYTHON_EXEC} setup.py develop 2>&1 ) 1>> ${LOG_FILE} || do_exit
touch done
cd $MY_PWD
-if !( ( ${DEST_DIR}/bin/python2.7 -c "import readline" 2>&1 )>> ${LOG_FILE})
+if !( ( ${DEST_DIR}/bin/${PYTHON_EXEC} -c "import readline" 2>&1 )>> ${LOG_FILE}) || \
+ [[ "${MYOS##Darwin}" != "${MYOS}" && $INST_PY3 -eq 1 ]]
then
- if !( ( ${DEST_DIR}/bin/python2.7 -c "import gnureadline" 2>&1 )>> ${LOG_FILE})
+ if !( ( ${DEST_DIR}/bin/${PYTHON_EXEC} -c "import gnureadline" 2>&1 )>> ${LOG_FILE})
then
echo "Installing pure-python readline"
( ${DEST_DIR}/bin/pip install gnureadline 2>&1 ) 1>> ${LOG_FILE}
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe doc/source/analyzing/generating_processed_data.rst
--- a/doc/source/analyzing/generating_processed_data.rst
+++ b/doc/source/analyzing/generating_processed_data.rst
@@ -47,10 +47,30 @@
frb = FixedResolutionBuffer(sl, (0.3, 0.5, 0.6, 0.8), (512, 512))
my_image = frb["density"]
-This resultant array can be saved out to disk or visualized using a
-hand-constructed Matplotlib image, for instance using
+This image may then be used in a hand-constructed Matplotlib image, for instance using
:func:`~matplotlib.pyplot.imshow`.
+The buffer arrays can be saved out to disk in either HDF5 or FITS format:
+
+.. code-block:: python
+
+ frb.export_hdf5("my_images.h5", fields=["density","temperature"])
+ frb.export_fits("my_images.fits", fields=["density","temperature"],
+ clobber=True, units="kpc")
+
+In the FITS case, there is an option for setting the ``units`` of the coordinate system in
+the file. If you want to overwrite a file with the same name, set ``clobber=True``.
+
+The :class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer` can even be exported
+as a 2D dataset itself, which may be operated on in the same way as any other dataset in yt:
+
+.. code-block:: python
+
+ ds_frb = frb.export_dataset(fields=["density","temperature"], nprocs=8)
+ sp = ds_frb.sphere("c", (100.,"kpc"))
+
+where the ``nprocs`` parameter can be used to decompose the image into ``nprocs`` number of grids.
+
.. _generating-profiles-and-histograms:
Profiles and Histograms
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe setup.py
--- a/setup.py
+++ b/setup.py
@@ -13,11 +13,6 @@
sys.exit(1)
import setuptools
-from distutils.version import StrictVersion
-if StrictVersion(setuptools.__version__) < StrictVersion('0.7.0'):
- import distribute_setup
- distribute_setup.use_setuptools()
-
from distutils.command.build_py import build_py
from numpy.distutils.misc_util import appendpath
from numpy.distutils.command import install_data as np_install_data
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe yt/analysis_modules/absorption_spectrum/absorption_line.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_line.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_line.py
@@ -16,12 +16,9 @@
import numpy as np
from yt.utilities.physical_constants import \
charge_proton_cgs, \
- cm_per_km, \
- km_per_cm, \
mass_electron_cgs, \
speed_of_light_cgs
-
def voigt(a,u):
"""
NAME:
@@ -140,67 +137,75 @@
k1 = k1.astype(np.float64).clip(0)
return k1
-def tau_profile(lam0, fval, gamma, vkms, column_density,
- deltav=None, delta_lambda=None,
+def tau_profile(lamba_0, f_value, gamma, v_doppler, column_density,
+ delta_v=None, delta_lambda=None,
lambda_bins=None, n_lambda=12000, dlambda=0.01):
- """
+ r"""
Create an optical depth vs. wavelength profile for an
absorption line using a voigt profile.
- :param lam0 (float): central wavelength (angstroms).
- :param fval (float): f-value.
- :param gamma (float): gamma value.
- :param vkms (float): doppler b-parameter.
- :param column_density (float): column density (cm^-2).
- :param deltav (float): velocity offset from lam0 (km/s).
- Default: None (no shift).
- :param delta_lambda (float): wavelength offset in angstroms.
- Default: None (no shift).
- :param lambda_bins (array): array of wavelengths in angstroms.
- Default: None
- :param n_lambda (float): size of lambda bins to create
- array if lambda_bins is None. Default: 12000
- :param dlambda (float): lambda bin width if lambda_bins is
- None. Default: 0.01
+
+ Parameters
+ ----------
+
+ lamba_0 : float YTQuantity in length units
+ central wavelength.
+ f_value : float
+ absorption line f-value.
+ gamma : float
+ absorption line gamma value.
+ v_doppler : float YTQuantity in velocity units
+ doppler b-parameter.
+ column_density : float YTQuantity in (length units)^-2
+ column density.
+ delta_v : float YTQuantity in velocity units
+ velocity offset from lamba_0.
+ Default: None (no shift).
+ delta_lambda : float YTQuantity in length units
+ wavelength offset.
+ Default: None (no shift).
+ lambda_bins : YTArray in length units
+ wavelength array for line deposition. If None, one will be
+ created using n_lambda and dlambda.
+ Default: None.
+ n_lambda : int
+ size of lambda bins to create if lambda_bins is None.
+ Default: 12000.
+ dlambda : float
+ lambda bin width in angstroms if lambda_bins is None.
+ Default: 0.01.
+
"""
- ## constants
- me = mass_electron_cgs # grams mass electron
- e = charge_proton_cgs # esu
- c = speed_of_light_cgs * km_per_cm # km/s
- ccgs = speed_of_light_cgs # cm/s
-
- ## shift lam0 by deltav
- if deltav is not None:
- lam1 = lam0 * (1 + deltav / c)
+ ## shift lamba_0 by delta_v
+ if delta_v is not None:
+ lam1 = lamba_0 * (1 + delta_v / speed_of_light_cgs)
elif delta_lambda is not None:
- lam1 = lam0 + delta_lambda
+ lam1 = lamba_0 + delta_lambda
else:
- lam1 = lam0
+ lam1 = lamba_0
## conversions
- vdop = vkms * cm_per_km # in cm/s
- lam0cgs = lam0 / 1.e8 # rest wavelength in cm
- lam1cgs = lam1 / 1.e8 # line wavelength in cm
- nu1 = ccgs / lam1cgs # line freq in Hz
- nudop = vdop / ccgs * nu1 # doppler width in Hz
- lamdop = vdop / ccgs * lam1 # doppler width in Ang
+ nu1 = speed_of_light_cgs / lam1 # line freq in Hz
+ nudop = v_doppler / speed_of_light_cgs * nu1 # doppler width in Hz
+ lamdop = v_doppler / speed_of_light_cgs * lam1 # doppler width in Ang
## create wavelength
if lambda_bins is None:
lambda_bins = lam1 + \
np.arange(n_lambda, dtype=np.float) * dlambda - \
- n_lambda * dlambda / 2 # wavelength vector (angstroms)
- nua = ccgs / (lambda_bins / 1.e8) # frequency vector (Hz)
+ n_lambda * dlambda / 2 # wavelength vector (angstroms)
+ nua = (speed_of_light_cgs / lambda_bins) # frequency vector (Hz)
## tau_0
- tau_X = np.sqrt(np.pi) * e**2 / (me * ccgs) * \
- column_density * fval / vdop
- tau0 = tau_X * lam0cgs
+ tau_X = np.sqrt(np.pi) * charge_proton_cgs**2 / \
+ (mass_electron_cgs * speed_of_light_cgs) * \
+ column_density * f_value / v_doppler
+ tau0 = tau_X * lamba_0
# dimensionless frequency offset in units of doppler freq
- x = (nua - nu1) / nudop
- a = gamma / (4 * np.pi * nudop) # damping parameter
- phi = voigt(a, x) # profile
- tauphi = tau0 * phi # profile scaled with tau0
+ x = ((nua - nu1) / nudop).in_units("")
+ a = (gamma / (4 * np.pi * nudop)).in_units("s") # damping parameter
+ phi = voigt(a, x) # line profile
+ tauphi = (tau0 * phi).in_units("") # profile scaled with tau0
return (lambda_bins, tauphi)
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum.py
@@ -4,7 +4,7 @@
"""
-from __future__ import print_function
+
from __future__ import absolute_import
#-----------------------------------------------------------------------------
@@ -21,13 +21,12 @@
from .absorption_line import tau_profile
from yt.funcs import get_pbar, mylog
-from yt.units.yt_array import YTArray
+from yt.units.yt_array import YTArray, YTQuantity
from yt.utilities.physical_constants import \
- amu_cgs, boltzmann_constant_cgs, \
- speed_of_light_cgs, km_per_cm
+ boltzmann_constant_cgs, \
+ speed_of_light_cgs
from yt.utilities.on_demand_imports import _astropy
-speed_of_light_kms = speed_of_light_cgs * km_per_cm
pyfits = _astropy.pyfits
class AbsorptionSpectrum(object):
@@ -49,8 +48,10 @@
self.tau_field = None
self.flux_field = None
self.spectrum_line_list = None
- self.lambda_bins = np.linspace(lambda_min, lambda_max, n_lambda)
- self.bin_width = (lambda_max - lambda_min) / float(n_lambda - 1)
+ self.lambda_bins = YTArray(np.linspace(lambda_min, lambda_max, n_lambda),
+ "angstrom")
+ self.bin_width = YTQuantity((lambda_max - lambda_min) /
+ float(n_lambda - 1), "angstrom")
self.line_list = []
self.continuum_list = []
@@ -76,8 +77,10 @@
mass of atom in amu.
"""
self.line_list.append({'label': label, 'field_name': field_name,
- 'wavelength': wavelength, 'f_value': f_value,
- 'gamma': gamma, 'atomic_mass': atomic_mass,
+ 'wavelength': YTQuantity(wavelength, "angstrom"),
+ 'f_value': f_value,
+ 'gamma': gamma,
+ 'atomic_mass': YTQuantity(atomic_mass, "amu"),
'label_threshold': label_threshold})
def add_continuum(self, label, field_name, wavelength,
@@ -209,16 +212,15 @@
# include factor of (1 + z) because our velocity is in proper frame.
delta_lambda += line['wavelength'] * (1 + field_data['redshift']) * \
field_data['velocity_los'] / speed_of_light_cgs
- thermal_b = km_per_cm * np.sqrt((2 * boltzmann_constant_cgs *
- field_data['temperature']) /
- (amu_cgs * line['atomic_mass']))
- thermal_b.convert_to_cgs()
+ thermal_b = np.sqrt((2 * boltzmann_constant_cgs *
+ field_data['temperature']) /
+ line['atomic_mass'])
center_bins = np.digitize((delta_lambda + line['wavelength']),
self.lambda_bins)
# ratio of line width to bin width
width_ratio = ((line['wavelength'] + delta_lambda) * \
- thermal_b / speed_of_light_kms / self.bin_width).value
+ thermal_b / speed_of_light_cgs / self.bin_width).in_units("").d
if (width_ratio < 1.0).any():
mylog.warn(("%d out of %d line components are unresolved, " +
@@ -240,11 +242,13 @@
my_bin_ratio = spectrum_bin_ratio
while True:
lambda_bins, line_tau = \
- tau_profile(line['wavelength'], line['f_value'],
- line['gamma'], thermal_b[lixel],
- column_density[lixel],
- delta_lambda=delta_lambda[lixel],
- lambda_bins=self.lambda_bins[left_index[lixel]:right_index[lixel]])
+ tau_profile(
+ line['wavelength'], line['f_value'],
+ line['gamma'], thermal_b[lixel].in_units("km/s"),
+ column_density[lixel],
+ delta_lambda=delta_lambda[lixel],
+ lambda_bins=self.lambda_bins[left_index[lixel]:right_index[lixel]])
+
# Widen wavelength window until optical depth reaches a max value at the ends.
if (line_tau[0] < max_tau and line_tau[-1] < max_tau) or \
(left_index[lixel] <= 0 and right_index[lixel] >= self.n_lambda):
@@ -260,7 +264,7 @@
if line['label_threshold'] is not None and \
column_density[lixel] >= line['label_threshold']:
if use_peculiar_velocity:
- peculiar_velocity = km_per_cm * field_data['velocity_los'][lixel]
+ peculiar_velocity = field_data['velocity_los'][lixel].in_units("km/s")
else:
peculiar_velocity = 0.0
self.spectrum_line_list.append({'label': line['label'],
@@ -271,7 +275,7 @@
'redshift': field_data['redshift'][lixel],
'v_pec': peculiar_velocity})
pbar.update(i)
- pbar.finish()
+ pbar.finish()
del column_density, delta_lambda, thermal_b, \
center_bins, width_ratio, left_index, right_index
@@ -280,7 +284,7 @@
"""
Write out list of spectral lines.
"""
- print("Writing spectral line list: %s." % filename)
+ mylog.info("Writing spectral line list: %s." % filename)
self.spectrum_line_list.sort(key=lambda obj: obj['wavelength'])
f = open(filename, 'w')
f.write('#%-14s %-14s %-12s %-12s %-12s %-12s\n' %
@@ -295,7 +299,7 @@
"""
Write spectrum to an ascii file.
"""
- print("Writing spectrum to ascii file: %s." % filename)
+ mylog.info("Writing spectrum to ascii file: %s." % filename)
f = open(filename, 'w')
f.write("# wavelength[A] tau flux\n")
for i in range(self.lambda_bins.size):
@@ -307,7 +311,7 @@
"""
Write spectrum to a fits file.
"""
- print("Writing spectrum to fits file: %s." % filename)
+ mylog.info("Writing spectrum to fits file: %s." % filename)
col1 = pyfits.Column(name='wavelength', format='E', array=self.lambda_bins)
col2 = pyfits.Column(name='flux', format='E', array=self.flux_field)
cols = pyfits.ColDefs([col1, col2])
@@ -319,7 +323,7 @@
Write spectrum to an hdf5 file.
"""
- print("Writing spectrum to hdf5 file: %s." % filename)
+ mylog.info("Writing spectrum to hdf5 file: %s." % filename)
output = h5py.File(filename, 'w')
output.create_dataset('wavelength', data=self.lambda_bins)
output.create_dataset('tau', data=self.tau_field)
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe yt/analysis_modules/absorption_spectrum/absorption_spectrum_fit.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum_fit.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum_fit.py
@@ -1,9 +1,14 @@
-from __future__ import print_function
+import h5py
import numpy as np
-import h5py
-from yt.analysis_modules.absorption_spectrum.absorption_line \
- import voigt
-from yt.utilities.on_demand_imports import _scipy
+
+from yt.analysis_modules.absorption_spectrum.absorption_line import \
+ voigt
+from yt.funcs import \
+ mylog
+from yt.units.yt_array import \
+ YTArray
+from yt.utilities.on_demand_imports import \
+ _scipy
optimize = _scipy.optimize
@@ -79,6 +84,10 @@
absorption profiles. Same size as x.
"""
+ # convert to NumPy array if we have a YTArray
+ if isinstance(x, YTArray):
+ x = x.d
+
#Empty dictionary for fitted lines
allSpeciesLines = {}
@@ -1007,6 +1016,5 @@
f.create_dataset("{0}/b".format(ion),data=params['b'])
f.create_dataset("{0}/z".format(ion),data=params['z'])
f.create_dataset("{0}/complex".format(ion),data=params['group#'])
- print('Writing spectrum fit to {0}'.format(file_name))
+ mylog.info('Writing spectrum fit to {0}'.format(file_name))
f.close()
-
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -709,7 +709,7 @@
pdata = pdata_ftype
# This will update the stream handler too
assign_particle_data(sds, pdata)
-
+
return sds
def load_amr_grids(grid_data, domain_dimensions,
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe yt/visualization/base_plot_types.py
--- a/yt/visualization/base_plot_types.py
+++ b/yt/visualization/base_plot_types.py
@@ -54,6 +54,7 @@
self._type_name = "CuttingPlane"
else:
self._type_name = viewer._plot_type
+ self.aspect = window_plot._aspect
self.font_properties = font_properties
self.font_color = font_color
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe yt/visualization/fixed_resolution.py
--- a/yt/visualization/fixed_resolution.py
+++ b/yt/visualization/fixed_resolution.py
@@ -20,6 +20,7 @@
from yt.utilities.lib.pixelization_routines import \
pixelize_cylinder
from yt.utilities.lib.api import add_points_to_greyscale_image
+from yt.frontends.stream.api import load_uniform_grid
from . import _MPL
import numpy as np
@@ -73,13 +74,13 @@
To make a projection and then several images, you can generate a
single FRB and then access multiple fields:
- >>> proj = ds.proj(0, "Density")
+ >>> proj = ds.proj(0, "density")
>>> frb1 = FixedResolutionBuffer(proj, (0.2, 0.3, 0.4, 0.5),
- (1024, 1024))
- >>> print frb1["Density"].max()
- 1.0914e-9
- >>> print frb1["Temperature"].max()
- 104923.1
+ ... (1024, 1024))
+ >>> print frb1["density"].max()
+ 1.0914e-9 g/cm**3
+ >>> print frb1["temperature"].max()
+ 104923.1 K
"""
_exclude_fields = ('pz','pdz','dx','x','y','z',
'r', 'dr', 'phi', 'dphi', 'theta', 'dtheta',
@@ -289,7 +290,7 @@
These fields will be pixelized and output.
"""
import h5py
- if fields is None: fields = self.data.keys()
+ if fields is None: fields = list(self.data.keys())
output = h5py.File(filename, "a")
for field in fields:
output.create_dataset(field,data=self[field])
@@ -307,30 +308,68 @@
filename : string
The name of the FITS file to be written.
fields : list of strings
- These fields will be pixelized and output.
+ These fields will be pixelized and output. If "None", the keys of the
+ FRB will be used.
clobber : boolean
If the file exists, this governs whether we will overwrite.
other_keys : dictionary, optional
A set of header keys and values to write into the FITS header.
units : string, optional
- the length units that the coordinates are written in, default 'cm'
- If units are set to "deg" then assume that sky coordinates are
- requested.
+ the length units that the coordinates are written in, default 'cm'.
"""
from yt.utilities.fits_image import FITSImageData
- extra_fields = ['x','y','z','px','py','pz','pdx','pdy','pdz','weight_field']
- if fields is None:
- fields = [field[-1] for field in self.data_source.field_data
- if field not in extra_fields]
+ if fields is None: fields = list(self.data.keys())
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)
fib.writeto(filename, clobber=clobber)
-
+
+ def export_dataset(self, fields=None, nprocs=1):
+ r"""Export a set of pixelized fields to an in-memory dataset that can be
+ analyzed as any other in yt. Unit information and other parameters (e.g.,
+ geometry, current_time, etc.) will be taken from the parent dataset.
+
+ Parameters
+ ----------
+ fields : list of strings, optional
+ These fields will be pixelized and output. If "None", the keys of the
+ FRB will be used.
+ nprocs: integer, optional
+ If greater than 1, will create this number of subarrays out of data
+
+ Examples
+ --------
+ >>> import yt
+ >>> ds = yt.load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150")
+ >>> slc = ds.slice(2, 0.0)
+ >>> frb = slc.to_frb((500.,"kpc"), 500)
+ >>> ds2 = frb.export_dataset(fields=["density","temperature"], nprocs=32)
+ """
+ nx, ny = self.buff_size
+ data = {}
+ if fields is None:
+ fields = list(self.keys())
+ for field in fields:
+ arr = self[field]
+ data[field] = (arr.d.T.reshape(nx,ny,1), str(arr.units))
+ bounds = [b.in_units("code_length").v for b in self.bounds]
+ bbox = np.array([[bounds[0],bounds[1]],[bounds[2],bounds[3]],[0.,1.]])
+ return load_uniform_grid(data, [nx,ny,1],
+ length_unit=self.ds.length_unit,
+ bbox=bbox,
+ sim_time=self.ds.current_time.in_units("s").v,
+ mass_unit=self.ds.mass_unit,
+ time_unit=self.ds.time_unit,
+ velocity_unit=self.ds.velocity_unit,
+ magnetic_unit=self.ds.magnetic_unit,
+ periodicity=(False,False,False),
+ geometry=self.ds.geometry,
+ nprocs=nprocs)
+
@property
def limits(self):
rv = dict(x = None, y = None, z = None)
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe yt/visualization/plot_container.py
--- a/yt/visualization/plot_container.py
+++ b/yt/visualization/plot_container.py
@@ -14,7 +14,9 @@
#-----------------------------------------------------------------------------
from yt.extern.six.moves import builtins
from yt.extern.six import iteritems
+
import base64
+import errno
import numpy as np
import matplotlib
import os
@@ -536,7 +538,13 @@
name = str(self.ds)
name = os.path.expanduser(name)
if name[-1] == os.sep and not os.path.isdir(name):
- os.mkdir(name)
+ try:
+ os.mkdir(name)
+ except OSError as e:
+ if e.errno == errno.EEXIST:
+ pass
+ else:
+ raise
if os.path.isdir(name) and name != str(self.ds):
name = name + (os.sep if name[-1] != os.sep else '') + str(self.ds)
if suffix is None:
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -23,6 +23,7 @@
from matplotlib.patches import Circle
from matplotlib.colors import colorConverter
+from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
from yt.funcs import *
from yt.extern.six import add_metaclass
@@ -1665,11 +1666,11 @@
# Setting pos overrides corner argument
if self.pos[0] is None or self.pos[1] is None:
if self.corner == 'upper_left':
- self.pos = (0.03, 0.97)
+ self.pos = (0.03, 0.96)
self.text_args['horizontalalignment'] = 'left'
self.text_args['verticalalignment'] = 'top'
elif self.corner == 'upper_right':
- self.pos = (0.97, 0.97)
+ self.pos = (0.97, 0.96)
self.text_args['horizontalalignment'] = 'right'
self.text_args['verticalalignment'] = 'top'
elif self.corner == 'lower_left':
@@ -1726,47 +1727,49 @@
class ScaleCallback(PlotCallback):
"""
annotate_scale(corner='lower_right', coeff=None, unit=None, pos=None,
- max_frac=0.2, min_frac=0.018, coord_system='axis',
- text_args=None, plot_args=None)
+ max_frac=0.16, min_frac=0.015, coord_system='axis',
+ size_bar_args=None, draw_inset_box=False,
+ inset_box_args=None)
Annotates the scale of the plot at a specified location in the image
(either in a preset corner, or by specifying (x,y) image coordinates with
- the pos argument. Coeff and units (e.g. 1 Mpc or 100 kpc) refer to the
- distance scale you desire to show on the plot. If no coeff and units are
- specified, an appropriate pair will be determined such that your scale bar
- is never smaller than min_frac or greater than max_frac of your plottable
- axis length. For additional text and plot arguments for the text and line,
- include them as dictionaries to pass to text_args and plot_args.
-
+ the pos argument. Coeff and units (e.g. 1 Mpc or 100 kpc) refer to the
+ distance scale you desire to show on the plot. If no coeff and units are
+ specified, an appropriate pair will be determined such that your scale bar
+ is never smaller than min_frac or greater than max_frac of your plottable
+ axis length. Additional customization of the scale bar is possible by
+ adjusting the size_bar_args dictionary. This accepts keyword arguments
+ for the AnchoredSizeBar class in matplotlib's axes_grid toolkit.
+
Parameters
----------
corner : string, optional
- Corner sets up one of 4 predeterimined locations for the timestamp
+ Corner sets up one of 4 predeterimined locations for the scale bar
to be displayed in the image: 'upper_left', 'upper_right', 'lower_left',
- 'lower_right' (also allows None). This value will be overridden by the
+ 'lower_right' (also allows None). This value will be overridden by the
optional 'pos' keyword.
coeff : float, optional
The coefficient of the unit defining the distance scale (e.g. 10 kpc or
- 100 Mpc) for overplotting. If set to None along with unit keyword,
+ 100 Mpc) for overplotting. If set to None along with unit keyword,
coeff will be automatically determined to be a power of 10
relative to the best-fit unit.
unit : string, optional
- unit must be a valid yt distance unit (e.g. 'm', 'km', 'AU', 'pc',
+ unit must be a valid yt distance unit (e.g. 'm', 'km', 'AU', 'pc',
'kpc', etc.) or set to None. If set to None, will be automatically
determined to be the best-fit to the data.
pos : 2- or 3-element tuples, lists, or arrays, optional
- The image location of the timestamp in the coord system defined by the
- coord_system kwarg. Setting pos overrides the corner parameter.
+ The image location of the scale bar in the plot coordinate system.
+ Setting pos overrides the corner parameter.
min_frac, max_frac: float, optional
- The minimum/maximum fraction of the axis width for the scale bar to
+ The minimum/maximum fraction of the axis width for the scale bar to
extend. A value of 1 would allow the scale bar to extend across the
- entire axis width. Only used for automatically calculating
- best-fit coeff and unit when neither is specified, otherwise
+ entire axis width. Only used for automatically calculating
+ best-fit coeff and unit when neither is specified, otherwise
disregarded.
coord_system : string, optional
@@ -1783,17 +1786,23 @@
"figure" -- the MPL figure coordinates: (0,0) is lower left, (1,1)
is upper right
- text_args : dictionary, optional
- A dictionary of any arbitrary parameters to be passed to the Matplotlib
- text object. Defaults: {'color':'white',
- 'horizontalalignment':'center', 'verticalalignment':'top'}.
+ size_bar_args : dictionary, optional
+ A dictionary of parameters to be passed to the Matplotlib
+ AnchoredSizeBar initializer.
+ Defaults: {'pad': 0.25, 'sep': 5, 'borderpad': 1, 'color': 'w'}
- plot_args : dictionary, optional
- A dictionary of any arbitrary parameters to be passed to the Matplotlib
- line object. Defaults: {'color':'white', 'linewidth':3}.
+ draw_inset_box : boolean, optional
+ Whether or not an inset box should be included around the scale bar.
+
+ inset_box_args : dictionary, optional
+ A dictionary of keyword arguments to be passed to the matplotlib Patch
+ object that represents the inset box.
+ Defaults: {'facecolor': 'black', 'linewidth': 3, 'edgecolor', 'white',
+ 'alpha': 0.5, 'boxstyle': 'square'}
+
Example
- -------
+ -------
>>> import yt
>>> ds = yt.load('Enzo_64/DD0020/data0020')
@@ -1802,11 +1811,23 @@
"""
_type_name = "scale"
def __init__(self, corner='lower_right', coeff=None, unit=None, pos=None,
- max_frac=0.20, min_frac=0.018, coord_system='axis',
- text_args=None, plot_args=None):
+ max_frac=0.16, min_frac=0.015, coord_system='axis',
+ size_bar_args=None, draw_inset_box=False, inset_box_args=None):
- def_text_args = {'color':'white'}
- def_plot_args = {'color':'white', 'linewidth':3}
+ def_size_bar_args = {
+ 'pad': 0.05,
+ 'sep': 5,
+ 'borderpad': 1,
+ 'color': 'w'
+ }
+
+ inset_box_args = {
+ 'facecolor': 'black',
+ 'linewidth': 3,
+ 'edgecolor': 'white',
+ 'alpha': 0.5,
+ 'boxstyle': 'square',
+ }
# Set position based on corner argument.
self.corner = corner
@@ -1816,33 +1837,34 @@
self.max_frac = max_frac
self.min_frac = min_frac
self.coord_system = coord_system
- if text_args is None: text_args = def_text_args
- self.text_args = text_args
- # This assures the line and the text are aligned
- self.text_args['horizontalalignment'] = 'center'
- self.text_args['verticalalignment'] = 'top'
- if plot_args is None: plot_args = def_plot_args
- self.plot_args = plot_args
+ if size_bar_args is None:
+ self.size_bar_args = def_size_bar_args
+ else:
+ self.size_bar_args = size_bar_args
+ if inset_box_args is None:
+ self.inset_box_args = def_inset_box_args
+ else:
+ self.inset_box_args = inset_box_args
+ self.draw_inset_box = draw_inset_box
def __call__(self, plot):
# Callback only works for plots with axis ratios of 1
xsize = plot.xlim[1] - plot.xlim[0]
- ysize = plot.ylim[1] - plot.ylim[0]
- if xsize != ysize:
- raise RuntimeError("Scale callback only works for plots with "
- "axis ratios of 1. Here: xsize = %s, ysize "
- " = %s." % (xsize, ysize))
+ if plot.aspect != 1.0:
+ raise NotImplementedError(
+ "Scale callback has only been implemented for plots with no "
+ "aspect ratio scaling. (aspect = {%s})".format(plot._aspect))
# Setting pos overrides corner argument
if self.pos is None:
if self.corner == 'upper_left':
- self.pos = (0.12, 0.971)
+ self.pos = (0.11, 0.952)
elif self.corner == 'upper_right':
- self.pos = (0.88, 0.971)
+ self.pos = (0.89, 0.952)
elif self.corner == 'lower_left':
- self.pos = (0.12, 0.062)
+ self.pos = (0.11, 0.052)
elif self.corner == 'lower_right':
- self.pos = (0.88, 0.062)
+ self.pos = (0.89, 0.052)
elif self.corner is None:
self.pos = (0.5, 0.5)
else:
@@ -1851,8 +1873,8 @@
"'lower_right', or None")
# When identifying a best fit distance unit, do not allow scale marker
- # to be greater than max_frac fraction of xaxis or under min_frac
- # fraction of xaxis
+ # to be greater than max_frac fraction of xaxis or under min_frac
+ # fraction of xaxis
max_scale = self.max_frac * xsize
min_scale = self.min_frac * xsize
@@ -1861,30 +1883,38 @@
# If no units are set, then identify a best fit distance unit
if self.unit is None:
- min_scale = plot.ds.get_smallest_appropriate_unit(min_scale,
- return_quantity=True)
- max_scale = plot.ds.get_smallest_appropriate_unit(max_scale,
- return_quantity=True)
+ min_scale = plot.ds.get_smallest_appropriate_unit(
+ min_scale, return_quantity=True)
+ max_scale = plot.ds.get_smallest_appropriate_unit(
+ max_scale, return_quantity=True)
self.coeff = max_scale.v
self.unit = max_scale.units
self.scale = YTQuantity(self.coeff, self.unit)
- self.text = "{scale} {units}".format(scale=int(self.coeff),
- units=self.unit)
- image_scale = (plot.frb.convert_distance_x(self.scale) / \
+ text = "{scale} {units}".format(scale=int(self.coeff), units=self.unit)
+ image_scale = (plot.frb.convert_distance_x(self.scale) /
plot.frb.convert_distance_x(xsize)).v
- # This is just a fancy wrapper around the TextLabelCallback and the
- # ImageLineCallback
- pos_line_start = (self.pos[0]-image_scale/2, self.pos[1]+0.01)
- pos_line_end = (self.pos[0]+image_scale/2, self.pos[1]+0.01)
- icb = LinePlotCallback(pos_line_start, pos_line_end,
- coord_system=self.coord_system,
- plot_args=self.plot_args)
- icb(plot)
- tcb = TextLabelCallback(self.pos, self.text,
- coord_system=self.coord_system,
- text_args=self.text_args)
- return tcb(plot)
+ size_vertical = self.size_bar_args.pop('size_vertical', .005)
+ fontproperties = self.size_bar_args.pop(
+ 'fontproperties', plot.font_properties)
+ frameon = self.size_bar_args.pop('frameon', self.draw_inset_box)
+
+ # this "anchors" the size bar to a box centered on self.pos in axis
+ # coordinates
+ self.size_bar_args['bbox_to_anchor'] = self.pos
+ self.size_bar_args['bbox_transform'] = plot._axes.transAxes
+
+ bar = AnchoredSizeBar(plot._axes.transAxes, image_scale, text, 10,
+ size_vertical=size_vertical,
+ fontproperties=fontproperties,
+ frameon=frameon,
+ **self.size_bar_args)
+
+ bar.patch.set(**self.inset_box_args)
+
+ plot._axes.add_artist(bar)
+
+ return plot
class RayCallback(PlotCallback):
"""
diff -r 16985637841f11d9bb814ad0f2edcc9391cbcf60 -r 45467416f055ff41bac76b5d90f392e5f0f5cafe yt/visualization/tests/test_export_frb.py
--- /dev/null
+++ b/yt/visualization/tests/test_export_frb.py
@@ -0,0 +1,39 @@
+"""
+Tests for exporting an FRB as a dataset
+
+
+
+"""
+from __future__ import absolute_import
+
+#-----------------------------------------------------------------------------
+# 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 numpy as np
+from yt.testing import \
+ fake_random_ds, assert_equal, \
+ assert_allclose
+
+def setup():
+ """Test specific setup."""
+ from yt.config import ytcfg
+ ytcfg["yt", "__withintesting"] = "True"
+
+
+def test_export_frb():
+ test_ds = fake_random_ds(128)
+ slc = test_ds.slice(0,0.5)
+ frb = slc.to_frb((0.5,"unitary"), 64)
+ frb_ds = frb.export_dataset(fields=["density"], nprocs=8)
+ dd_frb = frb_ds.all_data()
+
+ yield assert_equal, frb_ds.domain_left_edge.v, np.array([0.25,0.25,0.0])
+ yield assert_equal, frb_ds.domain_right_edge.v, np.array([0.75,0.75,1.0])
+ yield assert_equal, frb_ds.domain_width.v, np.array([0.5,0.5,1.0])
+ yield assert_equal, frb_ds.domain_dimensions, np.array([64,64,1], dtype="int64")
+ yield assert_allclose, frb["density"].sum(), dd_frb.quantities.total_quantity("density")
+ yield assert_equal, frb_ds.index.num_grids, 8
https://bitbucket.org/yt_analysis/yt/commits/a8be9a38fde1/
Changeset: a8be9a38fde1
Branch: yt
User: jzuhone
Date: 2015-06-03 16:47:01+00:00
Summary: Replacing this here as well
Affected #: 1 file
diff -r 45467416f055ff41bac76b5d90f392e5f0f5cafe -r a8be9a38fde15f42a2aec4d2bac7015b03a85ec8 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)
https://bitbucket.org/yt_analysis/yt/commits/2ab401515cc1/
Changeset: 2ab401515cc1
Branch: yt
User: jzuhone
Date: 2015-06-11 21:39:58+00:00
Summary: Fixing a few things noticed by Matt
Affected #: 2 files
diff -r a8be9a38fde15f42a2aec4d2bac7015b03a85ec8 -r 2ab401515cc14a74d47d79ff9f0a788aa15241a7 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
@@ -312,8 +312,8 @@
def _create_intensity(self):
def _intensity(field, data):
- v = self.current_v-data["v_los"].v
- T = data["temperature"].v
+ 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
diff -r a8be9a38fde15f42a2aec4d2bac7015b03a85ec8 -r 2ab401515cc14a74d47d79ff9f0a788aa15241a7 yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -163,7 +163,7 @@
if isinstance(width[0], YTQuantity):
cdelt = [wh.in_units(units).v/n for wh, n in zip(width, self.shape)]
else:
- cdelt = [wh/n for wh, n in zip(width, self.shape)]
+ 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
@@ -429,7 +429,7 @@
width = ds.coordinates.sanitize_width(axis, width, None)
unit = str(width[0].units)
if image_res is None:
- ddims = ds.domain_dimensions*2**ds.index.max_level
+ ddims = ds.domain_dimensions*ds.refine_by**ds.index.max_level
if iterable(axis):
nx = ddims.max()
ny = ddims.max()
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