[yt-svn] commit/yt: MatthewTurk: Merged in jzuhone/yt-3.x/yt-3.0 (pull request #879)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon May 5 19:13:14 PDT 2014
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/907cb91838fe/
Changeset: 907cb91838fe
Branch: yt-3.0
User: MatthewTurk
Date: 2014-05-06 04:13:06
Summary: Merged in jzuhone/yt-3.x/yt-3.0 (pull request #879)
Round-trip unit conversions to AstroPy, on-demand imports module, Dataset-specific add_field
Affected #: 13 files
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f doc/source/analyzing/units/2)_Data_Selection_and_fields.ipynb
--- a/doc/source/analyzing/units/2)_Data_Selection_and_fields.ipynb
+++ b/doc/source/analyzing/units/2)_Data_Selection_and_fields.ipynb
@@ -1,7 +1,7 @@
{
"metadata": {
"name": "",
- "signature": "sha256:9e7ac626b3609cf5f3fb2d4ebc6e027ed923ab1c22f0acc212e42fc7535e3205"
+ "signature": "sha256:b7541e0167001c6dd74306c8490385ace7bdb0533a829286f0505c0b24c67f16"
},
"nbformat": 3,
"nbformat_minor": 0,
@@ -296,6 +296,166 @@
"language": "python",
"metadata": {},
"outputs": []
+ },
+ {
+ "cell_type": "heading",
+ "level": 3,
+ "metadata": {},
+ "source": [
+ "Round-Trip Conversions to and from AstroPy's Units System"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Finally, a `YTArray` or `YTQuantity` may be converted to an [AstroPy quantity](http://astropy.readthedocs.org/en/latest/units/), which is a NumPy array or a scalar associated with units from AstroPy's units system. You may use this facility if you have AstroPy installed. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Some examples of converting from AstroPy units to yt:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "from astropy import units as u\n",
+ "x = 42.0 * u.meter\n",
+ "y = YTQuantity(x)\n",
+ "y2 = YTQuantity.from_astropy(x) # Another way to create the quantity"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "print x, type(x)\n",
+ "print y, type(y)\n",
+ "print y2, type(y2)"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "a = np.random.random(size=10) * u.km/u.s\n",
+ "b = YTArray(a)\n",
+ "b2 = YTArray.from_astropy(a) # Another way to create the quantity"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "print a, type(a)\n",
+ "print b, type(b)\n",
+ "print b2, type(b2)"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "It also works the other way around, converting a `YTArray` or `YTQuantity` to an AstroPy quantity via the method `to_astropy`. For arrays:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "temp = dd[\"temperature\"]\n",
+ "atemp = temp.to_astropy()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "print temp, type(temp)\n",
+ "print atemp, type(atemp)"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and quantities:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "from yt.utilities.physical_constants import kboltz\n",
+ "kb = kboltz.to_astropy()"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "print kboltz, type(kboltz)\n",
+ "print kb, type(kb)"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As a sanity check, you can show that it works round-trip:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "k1 = kboltz.to_astropy()\n",
+ "k2 = YTQuantity(kb)\n",
+ "print k1 == k2"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "c = YTArray(a)\n",
+ "d = c.to_astropy()\n",
+ "print a == d"
+ ],
+ "language": "python",
+ "metadata": {},
+ "outputs": []
}
],
"metadata": {}
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/analysis_modules/photon_simulator/photon_simulator.py
--- a/yt/analysis_modules/photon_simulator/photon_simulator.py
+++ b/yt/analysis_modules/photon_simulator/photon_simulator.py
@@ -21,22 +21,19 @@
#-----------------------------------------------------------------------------
import numpy as np
-from numpy.testing import assert_allclose
from yt.funcs import *
-from yt.utilities.physical_constants import clight, \
- cm_per_km, erg_per_keV
+from yt.utilities.physical_constants import clight
from yt.utilities.cosmology import Cosmology
from yt.utilities.orientation import Orientation
-from yt.utilities.definitions import mpc_conversion
from yt.utilities.parallel_tools.parallel_analysis_interface import \
communication_system, parallel_root_only, get_mpi_type, \
- op_names, parallel_capable
+ parallel_capable
from yt import units
from yt.units.yt_array import YTQuantity
import h5py
-from yt.frontends.fits.data_structures import ap
-pyfits = ap.pyfits
-pywcs = ap.pywcs
+from yt.utilities.on_demand_imports import _astropy
+pyfits = _astropy.pyfits
+pywcs = _astropy.pywcs
comm = communication_system.communicators[-1]
@@ -471,7 +468,9 @@
"but not both!")
if sky_center is None:
- sky_center = YTArray([30.,45.], "degree")
+ sky_center = YTArray([30.,45.], "degree")
+ else:
+ sky_center = YTArray(sky_center, "degree")
dx = self.photons["dx"].ndarray_view()
nx = self.parameters["Dimension"]
@@ -897,39 +896,42 @@
tbhdu = pyfits.new_table(coldefs)
tbhdu.update_ext_name("EVENTS")
- tbhdu.header.update("MTYPE1", "sky")
- tbhdu.header.update("MFORM1", "x,y")
- tbhdu.header.update("MTYPE2", "EQPOS")
- tbhdu.header.update("MFORM2", "RA,DEC")
- tbhdu.header.update("TCTYP2", "RA---TAN")
- tbhdu.header.update("TCTYP3", "DEC--TAN")
- tbhdu.header.update("TCRVL2", float(self.parameters["sky_center"][0]))
- tbhdu.header.update("TCRVL3", float(self.parameters["sky_center"][1]))
- tbhdu.header.update("TCDLT2", -float(self.parameters["dtheta"]))
- tbhdu.header.update("TCDLT3", float(self.parameters["dtheta"]))
- tbhdu.header.update("TCRPX2", self.parameters["pix_center"][0])
- tbhdu.header.update("TCRPX3", self.parameters["pix_center"][1])
- tbhdu.header.update("TLMIN2", 0.5)
- tbhdu.header.update("TLMIN3", 0.5)
- tbhdu.header.update("TLMAX2", 2.*self.parameters["pix_center"][0]-0.5)
- tbhdu.header.update("TLMAX3", 2.*self.parameters["pix_center"][1]-0.5)
- tbhdu.header.update("EXPOSURE", float(self.parameters["ExposureTime"]))
- tbhdu.header.update("AREA", float(self.parameters["Area"]))
- tbhdu.header.update("D_A", float(self.parameters["AngularDiameterDistance"]))
- tbhdu.header.update("REDSHIFT", self.parameters["Redshift"])
- tbhdu.header.update("HDUVERS", "1.1.0")
- tbhdu.header.update("RADECSYS", "FK5")
- tbhdu.header.update("EQUINOX", 2000.0)
+ tbhdu.header["MTYPE1"] = "sky"
+ tbhdu.header["MFORM1"] = "x,y"
+ tbhdu.header["MTYPE2"] = "EQPOS"
+ tbhdu.header["MFORM2"] = "RA,DEC"
+ tbhdu.header["TCTYP2"] = "RA---TAN"
+ tbhdu.header["TCTYP3"] = "DEC--TAN"
+ tbhdu.header["TCRVL2"] = float(self.parameters["sky_center"][0])
+ tbhdu.header["TCRVL3"] = float(self.parameters["sky_center"][1])
+ tbhdu.header["TCDLT2"] = -float(self.parameters["dtheta"])
+ tbhdu.header["TCDLT3"] = float(self.parameters["dtheta"])
+ tbhdu.header["TCRPX2"] = self.parameters["pix_center"][0]
+ tbhdu.header["TCRPX3"] = self.parameters["pix_center"][1]
+ tbhdu.header["TLMIN2"] = 0.5
+ tbhdu.header["TLMIN3"] = 0.5
+ tbhdu.header["TLMAX2"] = 2.*self.parameters["pix_center"][0]-0.5
+ tbhdu.header["TLMAX3"] = 2.*self.parameters["pix_center"][1]-0.5
+ tbhdu.header["EXPOSURE"] = float(self.parameters["ExposureTime"])
+ if isinstance(self.parameters["Area"], basestring):
+ tbhdu.header["AREA"] = self.parameters["Area"]
+ else:
+ tbhdu.header["AREA"] = float(self.parameters["Area"])
+ tbhdu.header["D_A"] = float(self.parameters["AngularDiameterDistance"])
+ tbhdu.header["REDSHIFT"] = self.parameters["Redshift"]
+ tbhdu.header["HDUVERS"] = "1.1.0"
+ tbhdu.header["RADECSYS"] = "FK5"
+ tbhdu.header["EQUINOX"] = 2000.0
if "RMF" in self.parameters:
- tbhdu.header.update("RMF", self.parameters["RMF"])
+ tbhdu.header["RMF"] = self.parameters["RMF"]
if "ARF" in self.parameters:
- tbhdu.header.update("ARF", self.parameters["ARF"])
+ tbhdu.header["ARF"] = self.parameters["ARF"]
if "ChannelType" in self.parameters:
- tbhdu.header.update("CHANTYPE", self.parameters["ChannelType"])
+ tbhdu.header["CHANTYPE"] = self.parameters["ChannelType"]
if "Telescope" in self.parameters:
- tbhdu.header.update("TELESCOP", self.parameters["Telescope"])
+ tbhdu.header["TELESCOP"] = self.parameters["Telescope"]
if "Instrument" in self.parameters:
- tbhdu.header.update("INSTRUME", self.parameters["Instrument"])
+ tbhdu.header["INSTRUME"] = self.parameters["Instrument"]
tbhdu.writeto(fitsfile, clobber=clobber)
@@ -978,15 +980,15 @@
tbhdu = pyfits.new_table(coldefs)
tbhdu.update_ext_name("PHLIST")
- tbhdu.header.update("HDUCLASS", "HEASARC/SIMPUT")
- tbhdu.header.update("HDUCLAS1", "PHOTONS")
- tbhdu.header.update("HDUVERS", "1.1.0")
- tbhdu.header.update("EXTVER", 1)
- tbhdu.header.update("REFRA", 0.0)
- tbhdu.header.update("REFDEC", 0.0)
- tbhdu.header.update("TUNIT1", "keV")
- tbhdu.header.update("TUNIT2", "deg")
- tbhdu.header.update("TUNIT3", "deg")
+ tbhdu.header["HDUCLASS"] = "HEASARC/SIMPUT"
+ tbhdu.header["HDUCLAS1"] = "PHOTONS"
+ tbhdu.header["HDUVERS"] = "1.1.0"
+ tbhdu.header["EXTVER"] = 1
+ tbhdu.header["REFRA"] = 0.0
+ tbhdu.header["REFDEC"] = 0.0
+ tbhdu.header["TUNIT1"] = "keV"
+ tbhdu.header["TUNIT2"] = "deg"
+ tbhdu.header["TUNIT3"] = "deg"
phfile = prefix+"_phlist.fits"
@@ -1006,17 +1008,17 @@
wrhdu = pyfits.new_table(coldefs)
wrhdu.update_ext_name("SRC_CAT")
- wrhdu.header.update("HDUCLASS", "HEASARC")
- wrhdu.header.update("HDUCLAS1", "SIMPUT")
- wrhdu.header.update("HDUCLAS2", "SRC_CAT")
- wrhdu.header.update("HDUVERS", "1.1.0")
- wrhdu.header.update("RADECSYS", "FK5")
- wrhdu.header.update("EQUINOX", 2000.0)
- wrhdu.header.update("TUNIT2", "deg")
- wrhdu.header.update("TUNIT3", "deg")
- wrhdu.header.update("TUNIT4", "keV")
- wrhdu.header.update("TUNIT5", "keV")
- wrhdu.header.update("TUNIT6", "erg/s/cm**2")
+ wrhdu.header["HDUCLASS"] = "HEASARC"
+ wrhdu.header["HDUCLAS1"] = "SIMPUT"
+ wrhdu.header["HDUCLAS2"] = "SRC_CAT"
+ wrhdu.header["HDUVERS"] = "1.1.0"
+ wrhdu.header["RADECSYS"] = "FK5"
+ wrhdu.header["EQUINOX"] = 2000.0
+ wrhdu.header["TUNIT2"] = "deg"
+ wrhdu.header["TUNIT3"] = "deg"
+ wrhdu.header["TUNIT4"] = "keV"
+ wrhdu.header["TUNIT5"] = "keV"
+ wrhdu.header["TUNIT6"] = "erg/s/cm**2"
simputfile = prefix+"_simput.fits"
@@ -1100,19 +1102,19 @@
hdu = pyfits.PrimaryHDU(H.T)
- hdu.header.update("MTYPE1", "EQPOS")
- hdu.header.update("MFORM1", "RA,DEC")
- hdu.header.update("CTYPE1", "RA---TAN")
- hdu.header.update("CTYPE2", "DEC--TAN")
- hdu.header.update("CRPIX1", 0.5*(nx+1))
- hdu.header.update("CRPIX2", 0.5*(nx+1))
- hdu.header.update("CRVAL1", float(self.parameters["sky_center"][0]))
- hdu.header.update("CRVAL2", float(self.parameters["sky_center"][1]))
- hdu.header.update("CUNIT1", "deg")
- hdu.header.update("CUNIT2", "deg")
- hdu.header.update("CDELT1", -float(self.parameters["dtheta"]))
- hdu.header.update("CDELT2", float(self.parameters["dtheta"]))
- hdu.header.update("EXPOSURE", float(self.parameters["ExposureTime"]))
+ hdu.header["MTYPE1"] = "EQPOS"
+ hdu.header["MFORM1"] = "RA,DEC"
+ hdu.header["CTYPE1"] = "RA---TAN"
+ hdu.header["CTYPE2"] = "DEC--TAN"
+ hdu.header["CRPIX1"] = 0.5*(nx+1)
+ hdu.header["CRPIX2"] = 0.5*(nx+1)
+ hdu.header["CRVAL1"] = float(self.parameters["sky_center"][0])
+ hdu.header["CRVAL2"] = float(self.parameters["sky_center"][1])
+ hdu.header["CUNIT1"] = "deg"
+ hdu.header["CUNIT2"] = "deg"
+ hdu.header["CDELT1"] = -float(self.parameters["dtheta"])
+ hdu.header["CDELT2"] = float(self.parameters["dtheta"])
+ hdu.header["EXPOSURE"] = float(self.parameters["ExposureTime"])
hdu.writeto(imagefile, clobber=clobber)
@@ -1183,41 +1185,41 @@
tbhdu.update_ext_name("SPECTRUM")
if not energy_bins:
- tbhdu.header.update("DETCHANS", spec.shape[0])
- tbhdu.header.update("TOTCTS", spec.sum())
- tbhdu.header.update("EXPOSURE", float(self.parameters["ExposureTime"]))
- tbhdu.header.update("LIVETIME", float(self.parameters["ExposureTime"]))
- tbhdu.header.update("CONTENT", spectype)
- tbhdu.header.update("HDUCLASS", "OGIP")
- tbhdu.header.update("HDUCLAS1", "SPECTRUM")
- tbhdu.header.update("HDUCLAS2", "TOTAL")
- tbhdu.header.update("HDUCLAS3", "TYPE:I")
- tbhdu.header.update("HDUCLAS4", "COUNT")
- tbhdu.header.update("HDUVERS", "1.1.0")
- tbhdu.header.update("HDUVERS1", "1.1.0")
- tbhdu.header.update("CHANTYPE", spectype)
- tbhdu.header.update("BACKFILE", "none")
- tbhdu.header.update("CORRFILE", "none")
- tbhdu.header.update("POISSERR", True)
+ tbhdu.header["DETCHANS"] = spec.shape[0]
+ tbhdu.header["TOTCTS"] = spec.sum()
+ tbhdu.header["EXPOSURE"] = float(self.parameters["ExposureTime"])
+ tbhdu.header["LIVETIME"] = float(self.parameters["ExposureTime"])
+ tbhdu.header["CONTENT"] = spectype
+ tbhdu.header["HDUCLASS"] = "OGIP"
+ tbhdu.header["HDUCLAS1"] = "SPECTRUM"
+ tbhdu.header["HDUCLAS2"] = "TOTAL"
+ tbhdu.header["HDUCLAS3"] = "TYPE:I"
+ tbhdu.header["HDUCLAS4"] = "COUNT"
+ tbhdu.header["HDUVERS"] = "1.1.0"
+ tbhdu.header["HDUVERS1"] = "1.1.0"
+ tbhdu.header["CHANTYPE"] = spectype
+ tbhdu.header["BACKFILE"] = "none"
+ tbhdu.header["CORRFILE"] = "none"
+ tbhdu.header["POISSERR"] = True
if self.parameters.has_key("RMF"):
- tbhdu.header.update("RESPFILE", self.parameters["RMF"])
+ tbhdu.header["RESPFILE"] = self.parameters["RMF"]
else:
- tbhdu.header.update("RESPFILE", "none")
+ tbhdu.header["RESPFILE"] = "none"
if self.parameters.has_key("ARF"):
- tbhdu.header.update("ANCRFILE", self.parameters["ARF"])
+ tbhdu.header["ANCRFILE"] = self.parameters["ARF"]
else:
- tbhdu.header.update("ANCRFILE", "none")
+ tbhdu.header["ANCRFILE"] = "none"
if self.parameters.has_key("Telescope"):
- tbhdu.header.update("TELESCOP", self.parameters["Telescope"])
+ tbhdu.header["TELESCOP"] = self.parameters["Telescope"]
else:
- tbhdu.header.update("TELESCOP", "none")
+ tbhdu.header["TELESCOP"] = "none"
if self.parameters.has_key("Instrument"):
- tbhdu.header.update("INSTRUME", self.parameters["Instrument"])
+ tbhdu.header["INSTRUME"] = self.parameters["Instrument"]
else:
- tbhdu.header.update("INSTRUME", "none")
- tbhdu.header.update("AREASCAL", 1.0)
- tbhdu.header.update("CORRSCAL", 0.0)
- tbhdu.header.update("BACKSCAL", 1.0)
+ tbhdu.header["INSTRUME"] = "none"
+ tbhdu.header["AREASCAL"] = 1.0
+ tbhdu.header["CORRSCAL"] = 0.0
+ tbhdu.header["BACKSCAL"] = 1.0
hdulist = pyfits.HDUList([pyfits.PrimaryHDU(), tbhdu])
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/analysis_modules/photon_simulator/spectral_models.py
--- a/yt/analysis_modules/photon_simulator/spectral_models.py
+++ b/yt/analysis_modules/photon_simulator/spectral_models.py
@@ -28,8 +28,8 @@
from scipy import stats
except ImportError:
pass
-from yt.frontends.fits.data_structures import ap
-pyfits = ap.pyfits
+from yt.utilities.on_demand_imports import _astropy
+pyfits = _astropy.pyfits
from yt.utilities.physical_constants import hcgs, clight, erg_per_keV, amu_cgs
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f 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
@@ -11,7 +11,7 @@
#-----------------------------------------------------------------------------
import numpy as np
-from yt.frontends.fits.data_structures import ap
+from yt.utilities.on_demand_imports import _astropy
from yt.utilities.orientation import Orientation
from yt.utilities.fits_image import FITSImageBuffer
from yt.visualization.volume_rendering.camera import off_axis_projection
@@ -150,7 +150,7 @@
if length_unit[1] == "deg":
dx *= -1.
- w = ap.pywcs.WCS(naxis=3)
+ w = _astropy.pywcs.WCS(naxis=3)
w.wcs.crpix = [0.5*(self.nx+1), 0.5*(self.ny+1), 0.5*(self.nv+1)]
w.wcs.cdelt = [dx,dy,dv]
w.wcs.crval = [center[0], center[1], v_center]
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -653,6 +653,15 @@
registry = self.unit_registry)
return self._quan
+ def add_field(self, name, function=None, **kwargs):
+ """
+ Dataset-specific call to add_field
+ """
+ self.index
+ self.field_info.add_field(name, function=function, **kwargs)
+ deps, _ = self.field_info.check_derived_fields([name])
+ self.field_dependencies.update(deps)
+
def _reconstruct_pf(*args, **kwargs):
pfs = ParameterFileStore()
pf = pfs.get_pf_hash(*args)
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/frontends/fits/data_structures.py
--- a/yt/frontends/fits/data_structures.py
+++ b/yt/frontends/fits/data_structures.py
@@ -39,58 +39,7 @@
unit_prefixes
from yt.units import dimensions
from yt.units.yt_array import YTQuantity
-
-class astropy_imports:
- _pyfits = None
- @property
- def pyfits(self):
- if self._pyfits is None:
- try:
- import astropy.io.fits as pyfits
- self.log
- except ImportError:
- pyfits = None
- self._pyfits = pyfits
- return self._pyfits
-
- _pywcs = None
- @property
- def pywcs(self):
- if self._pywcs is None:
- try:
- import astropy.wcs as pywcs
- self.log
- except ImportError:
- pywcs = None
- self._pywcs = pywcs
- return self._pywcs
-
- _log = None
- @property
- def log(self):
- if self._log is None:
- try:
- from astropy import log
- if log.exception_logging_enabled():
- log.disable_exception_logging()
- except ImportError:
- log = None
- self._log = log
- return self._log
-
- _conv = None
- @property
- def conv(self):
- if self._conv is None:
- try:
- import astropy.convolution as conv
- self.log
- except ImportError:
- conv = None
- self._conv = conv
- return self._conv
-
-ap = astropy_imports()
+from yt.utilities.on_demand_imports import _astropy
lon_prefixes = ["X","RA","GLON"]
lat_prefixes = ["Y","DEC","GLAT"]
@@ -145,7 +94,8 @@
# FITS units always return upper-case, so we need to get
# the right case by comparing against known units. This
# only really works for common units.
- units = re.split(regex_pattern, field_units)
+ units = set(re.split(regex_pattern, field_units))
+ units.remove('')
n = int(0)
for unit in units:
if unit in known_units:
@@ -262,7 +212,7 @@
pf.domain_right_edge)])
dims = np.array(pf.domain_dimensions)
# If we are creating a dataset of lines, only decompose along the position axes
- if pf.line_database is not None:
+ if len(pf.line_database) > 0:
dims[pf.vel_axis] = 1
psize = get_psize(dims, pf.nprocs)
gle, gre, shapes, slices = decompose_array(dims, psize, bbox)
@@ -270,7 +220,7 @@
self.grid_right_edge = self.pf.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 pf.line_database is not None:
+ if len(pf.line_database) > 0:
self.grid_left_edge[:,pf.vel_axis] = pf.domain_left_edge[pf.vel_axis]
self.grid_right_edge[:,pf.vel_axis] = pf.domain_right_edge[pf.vel_axis]
self.grid_dimensions[:,pf.vel_axis] = pf.domain_dimensions[pf.vel_axis]
@@ -373,7 +323,7 @@
elif isinstance(nan_mask, dict):
self.nan_mask = nan_mask
self.nprocs = nprocs
- self._handle = ap.pyfits.open(self.filenames[0],
+ self._handle = _astropy.pyfits.open(self.filenames[0],
memmap=True,
do_not_scale_image_data=True,
ignore_blank=True)
@@ -384,7 +334,7 @@
fn = fits_file
else:
fn = os.path.join(ytcfg.get("yt","test_data_dir"),fits_file)
- f = ap.pyfits.open(fn, memmap=True,
+ f = _astropy.pyfits.open(fn, memmap=True,
do_not_scale_image_data=True,
ignore_blank=True)
self._fits_files.append(f)
@@ -394,7 +344,7 @@
self.first_image = 1
self.primary_header = self._handle[self.first_image].header
self.naxis = 2
- self.wcs = ap.pywcs.WCS(naxis=2)
+ self.wcs = _astropy.pywcs.WCS(naxis=2)
self.events_info = {}
for k,v in self.primary_header.items():
if k.startswith("TTYP"):
@@ -428,7 +378,7 @@
self.events_data = False
self.first_image = 0
self.primary_header = self._handle[self.first_image].header
- self.wcs = ap.pywcs.WCS(header=self.primary_header)
+ self.wcs = _astropy.pywcs.WCS(header=self.primary_header)
self.naxis = self.primary_header["naxis"]
self.axis_names = [self.primary_header["ctype%d" % (i+1)]
for i in xrange(self.naxis)]
@@ -531,6 +481,8 @@
32**self.dimensionality).astype("int")
self.nprocs = max(min(self.nprocs, 512), 1)
+ self.reversed = False
+
# Check to see if this data is in some kind of (Lat,Lon,Vel) format
self.ppv_data = False
x = 0
@@ -572,7 +524,7 @@
self.vel_axis = np.where(self.vel_axis)[0][0]
self.vel_name = ctypes[self.vel_axis].split("-")[0].lower()
- self.wcs_2d = ap.pywcs.WCS(naxis=2)
+ self.wcs_2d = _astropy.pywcs.WCS(naxis=2)
self.wcs_2d.wcs.crpix = self.wcs.wcs.crpix[[self.lon_axis, self.lat_axis]]
self.wcs_2d.wcs.cdelt = self.wcs.wcs.cdelt[[self.lon_axis, self.lat_axis]]
self.wcs_2d.wcs.crval = self.wcs.wcs.crval[[self.lon_axis, self.lat_axis]]
@@ -591,7 +543,6 @@
le = self.dims[self.vel_axis]+0.5
re = 0.5
else:
- self.reversed = False
le = 0.5
re = self.dims[self.vel_axis]+0.5
self.domain_left_edge[self.vel_axis] = (le-x0)*dz + z0
@@ -632,7 +583,7 @@
try:
with warnings.catch_warnings():
warnings.filterwarnings('ignore', category=UserWarning, append=True)
- fileh = ap.pyfits.open(args[0])
+ fileh = _astropy.pyfits.open(args[0])
valid = fileh[0].header["naxis"] >= 2
if len(fileh) > 1 and fileh[1].name == "EVENTS":
valid = fileh[1].header["naxis"] >= 2
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/frontends/fits/fields.py
--- a/yt/frontends/fits/fields.py
+++ b/yt/frontends/fits/fields.py
@@ -25,17 +25,17 @@
def _setup_ppv_fields(self):
- def _get_2d_wcs(self, data, axis):
+ def _get_2d_wcs(data, axis):
w_coords = data.pf.wcs_2d.wcs_pix2world(data["x"], data["y"], 1)
return w_coords[axis]
def world_f(axis, unit):
def _world_f(field, data):
- return data.pf.arr(self._get_2d_wcs(data, axis), unit)
+ return data.pf.arr(_get_2d_wcs(data, axis), unit)
return _world_f
for (i, axis), name in zip(enumerate([self.pf.lon_axis, self.pf.lat_axis]),
- [self.pf.lon_name, self.pf.lat_name]):
+ [self.pf.lon_name, self.pf.lat_name]):
unit = str(self.pf.wcs_2d.wcs.cunit[i])
if unit.lower() == "deg": unit = "degree"
if unit.lower() == "rad": unit = "radian"
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/frontends/fits/misc.py
--- a/yt/frontends/fits/misc.py
+++ b/yt/frontends/fits/misc.py
@@ -14,7 +14,7 @@
from yt.fields.api import add_field
from yt.fields.derived_field import ValidateSpatial
from yt.funcs import mylog
-from .data_structures import ap
+from yt.utilities.on_demand_imports import _astropy
def _make_counts(emin, emax):
def _counts(field, data):
@@ -30,31 +30,34 @@
else:
sigma = None
if sigma is not None and sigma > 0.0:
- kern = ap.conv.Gaussian2DKernel(stddev=sigma)
- img[:,:,0] = ap.conv.convolve(img[:,:,0], kern)
+ kern = _astropy.conv.Gaussian2DKernel(stddev=sigma)
+ img[:,:,0] = _astropy.conv.convolve(img[:,:,0], kern)
return data.pf.arr(img, "counts/pixel")
return _counts
-def setup_counts_fields(ebounds):
+def setup_counts_fields(ds, ebounds, ftype="gas"):
r"""
Create deposited image fields from X-ray count data in energy bands.
Parameters
----------
+ ds : Dataset
+ The FITS events file dataset to add the counts fields to.
ebounds : list of tuples
A list of tuples, one for each field, with (emin, emax) as the
energy bounds for the image.
Examples
--------
+ >>> ds = yt.load("evt.fits")
>>> ebounds = [(0.1,2.0),(2.0,3.0)]
- >>> setup_counts_fields(ebounds)
+ >>> setup_counts_fields(ds, ebounds)
"""
for (emin, emax) in ebounds:
cfunc = _make_counts(emin, emax)
fname = "counts_%s-%s" % (emin, emax)
mylog.info("Creating counts field %s." % fname)
- add_field(("gas",fname), function=cfunc,
- units="counts/pixel",
- validators = [ValidateSpatial()],
- display_name="Counts (%s-%s keV)" % (emin, emax))
\ No newline at end of file
+ ds.add_field((ftype,fname), function=cfunc,
+ units="counts/pixel",
+ validators = [ValidateSpatial()],
+ display_name="Counts (%s-%s keV)" % (emin, emax))
\ No newline at end of file
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/units/tests/test_ytarray.py
--- a/yt/units/tests/test_ytarray.py
+++ b/yt/units/tests/test_ytarray.py
@@ -26,7 +26,7 @@
unary_operators, binary_operators
from yt.utilities.exceptions import \
YTUnitOperationError, YTUfuncUnitError
-from yt.testing import fake_random_pf
+from yt.testing import fake_random_pf, requires_module
from yt.funcs import fix_length
import numpy as np
import copy
@@ -650,3 +650,28 @@
for op in [operator.abs, operator.neg, operator.pos]:
yield unary_op_registry_comparison, op
+
+ at requires_module("astropy")
+def test_astropy():
+ from yt.utilities.on_demand_imports import _astropy
+
+ ap_arr = np.arange(10)*_astropy.units.km/_astropy.units.hr
+ yt_arr = YTArray(np.arange(10), "km/hr")
+ yt_arr2 = YTArray.from_astropy(ap_arr)
+
+ ap_quan = 10.*_astropy.units.Msun**0.5/(_astropy.units.kpc**3)
+ yt_quan = YTQuantity(10.,"sqrt(Msun)/kpc**3")
+ yt_quan2 = YTQuantity.from_astropy(ap_quan)
+
+ yield assert_array_equal, ap_arr, yt_arr.to_astropy()
+ yield assert_array_equal, yt_arr, YTArray(ap_arr)
+ yield assert_array_equal, yt_arr, yt_arr2
+
+ yield assert_equal, ap_quan, yt_quan.to_astropy()
+ yield assert_equal, yt_quan, YTQuantity(ap_quan)
+ yield assert_equal, yt_quan, yt_quan2
+
+ yield assert_array_equal, yt_arr, YTArray(yt_arr.to_astropy())
+ yield assert_equal, yt_quan, YTQuantity(yt_quan.to_astropy())
+
+
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/units/unit_lookup_table.py
--- a/yt/units/unit_lookup_table.py
+++ b/yt/units/unit_lookup_table.py
@@ -59,12 +59,12 @@
"yr": (sec_per_year, dimensions.time),
# Solar units
- "Msun": ( mass_sun_grams, dimensions.mass),
- "msun": ( mass_sun_grams, dimensions.mass),
- "Rsun": ( cm_per_rsun, dimensions.length),
- "rsun": ( cm_per_rsun, dimensions.length),
- "Lsun": ( luminosity_sun_ergs_per_sec, dimensions.power),
- "Tsun": ( temp_sun_kelvin, dimensions.temperature),
+ "Msun": (mass_sun_grams, dimensions.mass),
+ "msun": (mass_sun_grams, dimensions.mass),
+ "Rsun": (cm_per_rsun, dimensions.length),
+ "rsun": (cm_per_rsun, dimensions.length),
+ "Lsun": (luminosity_sun_ergs_per_sec, dimensions.power),
+ "Tsun": (temp_sun_kelvin, dimensions.temperature),
"Zsun": (metallicity_sun, dimensions.dimensionless),
"Mjup": (mass_jupiter_grams, dimensions.mass),
"Mearth": (mass_earth_grams, dimensions.mass),
@@ -89,6 +89,20 @@
"angstrom": (cm_per_ang, dimensions.length),
"Jy": (jansky_cgs, dimensions.specific_flux),
"counts": (1.0, dimensions.dimensionless),
+
+ # for AstroPy compatibility
+ "solMass": (mass_sun_grams, dimensions.mass),
+ "solRad": (cm_per_rsun, dimensions.length),
+ "solLum": (luminosity_sun_ergs_per_sec, dimensions.power),
+ "dyn": (1.0, dimensions.force),
+ "sr": (1.0, dimensions.solid_angle),
+ "rad": (1.0, dimensions.solid_angle),
+ "deg": (np.pi/180., dimensions.angle),
+ "Fr": (1.0, dimensions.charge),
+ "G": (1.0, dimensions.magnetic_field),
+ "d": (1.0, dimensions.time),
+ "Angstrom": (cm_per_ang, dimensions.length),
+
}
# Add LaTeX representations for units with trivial representations.
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/units/yt_array.py
--- a/yt/units/yt_array.py
+++ b/yt/units/yt_array.py
@@ -35,6 +35,8 @@
YTUnitOperationError, YTUnitConversionError, \
YTUfuncUnitError
from numbers import Number as numeric_type
+from yt.utilities.on_demand_imports import _astropy
+from sympy import Rational
# redefine this here to avoid a circular import from yt.funcs
def iterable(obj):
@@ -249,6 +251,9 @@
"Perhaps you meant to do something like this instead: \n"
"ds.arr(%s, \"%s\")" % (input_array, input_units)
)
+ if _astropy.units is not None:
+ if isinstance(input_array, _astropy.units.quantity.Quantity):
+ return cls.from_astropy(input_array)
if isinstance(input_array, YTArray):
if input_units is None:
if registry is None:
@@ -423,6 +428,38 @@
"""
return np.array(self)
+
+ @classmethod
+ def from_astropy(cls, arr):
+ """
+ Creates a new YTArray with the same unit information from an
+ AstroPy quantity *arr*.
+ """
+ # Converting from AstroPy Quantity
+ u = arr.unit
+ ap_units = []
+ for base, power in zip(u.bases, u.powers):
+ unit_str = base.to_string()
+ # we have to do this because AstroPy is silly and defines
+ # hour as "h"
+ if unit_str == "h": unit_str = "hr"
+ ap_units.append("%s**(%s)" % (unit_str, Rational(power)))
+ ap_units = "*".join(ap_units)
+ if isinstance(arr.value, np.ndarray):
+ return YTArray(arr.value, ap_units)
+ else:
+ return YTQuantity(arr.value, ap_units)
+
+
+ def to_astropy(self, **kwargs):
+ """
+ Creates a new AstroPy quantity with the same unit information.
+ """
+ if _astropy.units is None:
+ raise ImportError("You don't have AstroPy installed, so you can't convert to " +
+ "an AstroPy quantity.")
+ return self.value*_astropy.units.Unit(str(self.units), **kwargs)
+
#
# End unit conversion methods
#
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -14,11 +14,11 @@
from yt.funcs import mylog, iterable, fix_axis, ensure_list
from yt.visualization.fixed_resolution import FixedResolutionBuffer
from yt.data_objects.construction_data_containers import YTCoveringGridBase
-from yt.frontends.fits.data_structures import ap
+from yt.utilities.on_demand_imports import _astropy
from yt.units.yt_array import YTQuantity
-pyfits = ap.pyfits
-pywcs = ap.pywcs
+pyfits = _astropy.pyfits
+pywcs = _astropy.pywcs
class FITSImageBuffer(pyfits.HDUList):
@@ -295,7 +295,7 @@
"""
def __init__(self, ds, axis, fields, coord, **kwargs):
fields = ensure_list(fields)
- axis = fix_axis(axis)
+ axis = fix_axis(axis, ds)
if isinstance(coord, tuple):
coord = ds.quan(coord[0], coord[1]).in_units("code_length").value
elif isinstance(coord, YTQuantity):
@@ -303,6 +303,8 @@
slc = ds.slice(axis, coord, **kwargs)
w, frb = construct_image(slc)
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):
r"""
@@ -321,10 +323,12 @@
"""
def __init__(self, ds, axis, fields, weight_field=None, **kwargs):
fields = ensure_list(fields)
- axis = fix_axis(axis)
+ axis = fix_axis(axis, ds)
prj = ds.proj(fields[0], axis, weight_field=weight_field, **kwargs)
w, frb = construct_image(prj)
super(FITSProjection, self).__init__(frb, fields=fields, wcs=w)
+ for i, field in enumerate(fields):
+ self[i].header["bunit"] = str(frb[field].units)
diff -r ffe8407ff931962371b1aed113f363fcaa412ff5 -r 907cb91838fe6a2ef52306d9a0c98dab2dc23c1f yt/utilities/on_demand_imports.py
--- /dev/null
+++ b/yt/utilities/on_demand_imports.py
@@ -0,0 +1,74 @@
+"""
+A set of convenient on-demand imports
+"""
+
+#-----------------------------------------------------------------------------
+# 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.
+#-----------------------------------------------------------------------------
+
+class astropy_imports:
+ _pyfits = None
+ @property
+ def pyfits(self):
+ if self._pyfits is None:
+ try:
+ import astropy.io.fits as pyfits
+ self.log
+ except ImportError:
+ pyfits = None
+ self._pyfits = pyfits
+ return self._pyfits
+
+ _pywcs = None
+ @property
+ def pywcs(self):
+ if self._pywcs is None:
+ try:
+ import astropy.wcs as pywcs
+ self.log
+ except ImportError:
+ pywcs = None
+ self._pywcs = pywcs
+ return self._pywcs
+
+ _log = None
+ @property
+ def log(self):
+ if self._log is None:
+ try:
+ from astropy import log
+ if log.exception_logging_enabled():
+ log.disable_exception_logging()
+ except ImportError:
+ log = None
+ self._log = log
+ return self._log
+
+ _units = None
+ @property
+ def units(self):
+ if self._units is None:
+ try:
+ from astropy import units
+ except ImportError:
+ units = None
+ self._units = units
+ return self._units
+
+ _conv = None
+ @property
+ def conv(self):
+ if self._conv is None:
+ try:
+ import astropy.convolution as conv
+ self.log
+ except ImportError:
+ conv = None
+ self._conv = conv
+ return self._conv
+
+_astropy = astropy_imports()
\ No newline at end of file
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