[yt-svn] commit/yt: atmyers: Merged in jzuhone/yt-3.x (pull request #1588)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Fri May 29 12:20:36 PDT 2015
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/d6983e51e584/
Changeset: d6983e51e584
Branch: yt
User: atmyers
Date: 2015-05-29 19:20:26+00:00
Summary: Merged in jzuhone/yt-3.x (pull request #1588)
Removing line database stuff from the FITS frontend
Affected #: 5 files
diff -r 718444346962f551ca51fec7a6ef877b6badb473 -r d6983e51e584eb0e996dfd2a4475ff2a8e0a6014 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
^^^^^^^^^^^^^^^^^^^^^^^^^^^
diff -r 718444346962f551ca51fec7a6ef877b6badb473 -r d6983e51e584eb0e996dfd2a4475ff2a8e0a6014 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 718444346962f551ca51fec7a6ef877b6badb473 -r d6983e51e584eb0e996dfd2a4475ff2a8e0a6014 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 718444346962f551ca51fec7a6ef877b6badb473 -r d6983e51e584eb0e996dfd2a4475ff2a8e0a6014 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 718444346962f551ca51fec7a6ef877b6badb473 -r d6983e51e584eb0e996dfd2a4475ff2a8e0a6014 yt/utilities/file_handler.py
--- a/yt/utilities/file_handler.py
+++ b/yt/utilities/file_handler.py
@@ -54,12 +54,15 @@
class FITSFileHandler(HDF5FileHandler):
def __init__(self, filename):
from yt.utilities.on_demand_imports import _astropy
- if isinstance(filename, _astropy.pyfits.PrimaryHDU):
+ if isinstance(filename, _astropy.pyfits.hdu.image._ImageBaseHDU):
self.handle = _astropy.pyfits.HDUList(filename)
+ elif isinstance(filename, _astropy.pyfits.HDUList):
+ self.handle = filename
else:
self.handle = _astropy.pyfits.open(
filename, memmap=True, do_not_scale_image_data=True,
ignore_blank=True)
+ self._fits_files = []
def __del__(self):
for f in self._fits_files:
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