[yt-svn] commit/yt: ngoldbaum: Merged in jzuhone/yt (pull request #2455)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Dec 5 09:48:39 PST 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/b9c09a2edd4e/
Changeset: b9c09a2edd4e
Branch: yt
User: ngoldbaum
Date: 2016-12-05 17:48:12+00:00
Summary: Merged in jzuhone/yt (pull request #2455)
Allow multiple WCS coordinate systems in FITSImageData, a couple of new units, some simplifications
Affected #: 12 files
diff -r 49a03a144a9648a461207eb0942534f6762bea72 -r b9c09a2edd4e61caf2678093fcfdec8dceb0b20c doc/source/visualizing/FITSImageData.ipynb
--- a/doc/source/visualizing/FITSImageData.ipynb
+++ b/doc/source/visualizing/FITSImageData.ipynb
@@ -27,9 +27,9 @@
},
"outputs": [],
"source": [
- "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\")})"
+ "ds = yt.load(\"MHDSloshing/virgo_low_res.0054.vtk\", units_override={\"length_unit\":(1.0,\"Mpc\"),\n",
+ " \"mass_unit\":(1.0e14,\"Msun\"),\n",
+ " \"time_unit\":(1.0,\"Myr\")})"
]
},
{
@@ -348,7 +348,27 @@
"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:"
+ "So far, the FITS images we have shown have linear spatial coordinates. We can see this by looking at the header for one of the fields, and examining the `CTYPE1` and `CTYPE2` keywords:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prj_fits[\"temperature\"].header"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `WCSNAME` keyword is set to `\"yt\"` by default. \n",
+ "\n",
+ "However, 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:"
]
},
{
@@ -368,7 +388,7 @@
"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:"
+ "By 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:"
]
},
{
@@ -386,6 +406,49 @@
"cell_type": "markdown",
"metadata": {},
"source": [
+ "and now the `WCSNAME` has been set to `\"celestial\"`. If you don't want to override the default WCS but to add another one, then you can make the call to `create_sky_wcs` and set `replace_old_wcs=False`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "prj_fits3.create_sky_wcs(sky_center, sky_scale, ctype=[\"RA---TAN\",\"DEC--TAN\"], replace_old_wcs=False)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We now can see that there are two WCSes in the header, with the celestial WCS keywords having the \"A\" designation:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prj_fits3[\"temperature\"].header"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Any further WCSes that are added will have \"B\", \"C\", etc."
+ ]
+ },
+ {
+ "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`:"
]
},
@@ -415,22 +478,11 @@
}
],
"metadata": {
+ "anaconda-cloud": {},
"kernelspec": {
- "display_name": "Python 3",
+ "display_name": "Python [default]",
"language": "python",
"name": "python3"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.5.1"
}
},
"nbformat": 4,
diff -r 49a03a144a9648a461207eb0942534f6762bea72 -r b9c09a2edd4e61caf2678093fcfdec8dceb0b20c yt/analysis_modules/particle_trajectories/particle_trajectories.py
--- a/yt/analysis_modules/particle_trajectories/particle_trajectories.py
+++ b/yt/analysis_modules/particle_trajectories/particle_trajectories.py
@@ -99,7 +99,7 @@
pbar = get_pbar("Constructing trajectory information", len(self.data_series))
for i, (sto, ds) in enumerate(self.data_series.piter(storage=my_storage)):
dd = ds.all_data()
- newtags = dd[idx_field].ndarray_view().astype("int64")
+ newtags = dd[idx_field].d.astype("int64")
mask = np.in1d(newtags, indices, assume_unique=True)
sort = np.argsort(newtags[mask])
array_indices = np.where(np.in1d(indices, newtags, assume_unique=True))[0]
@@ -197,7 +197,6 @@
Examples
________
- >>> from yt.mods import *
>>> trajs = ParticleTrajectories(my_fns, indices)
>>> trajs.add_fields(["particle_mass", "particle_gpot"])
"""
@@ -247,15 +246,15 @@
dd = ds.all_data()
for field in new_particle_fields:
# This is easy... just get the particle fields
- pfield[field] = dd[fds[field]].ndarray_view()[mask][sort]
+ pfield[field] = dd[fds[field]].d[mask][sort]
if grid_fields:
# This is hard... must loop over grids
for field in grid_fields:
- pfield[field] = np.zeros((self.num_indices))
- x = self["particle_position_x"][:,step].ndarray_view()
- y = self["particle_position_y"][:,step].ndarray_view()
- z = self["particle_position_z"][:,step].ndarray_view()
+ pfield[field] = np.zeros(self.num_indices)
+ x = self["particle_position_x"][:,step].d
+ y = self["particle_position_y"][:,step].d
+ z = self["particle_position_z"][:,step].d
particle_grids, particle_grid_inds = ds.index._find_points(x,y,z)
# This will fail for non-grid index objects
@@ -375,10 +374,10 @@
>>> trajs.write_out_h5("orbit_trajectories")
"""
fid = h5py.File(filename, "w")
- fields = [field for field in sorted(self.field_data.keys())]
fid.create_dataset("particle_indices", dtype=np.int64,
data=self.indices)
- fid.create_dataset("particle_time", data=self.times)
+ fid.close()
+ self.times.write_hdf5(filename, dataset_name="particle_times")
+ fields = [field for field in sorted(self.field_data.keys())]
for field in fields:
- fid.create_dataset("%s" % field, data=self[field])
- fid.close()
+ self[field].write_hdf5(filename, dataset_name="%s" % field)
diff -r 49a03a144a9648a461207eb0942534f6762bea72 -r b9c09a2edd4e61caf2678093fcfdec8dceb0b20c yt/frontends/fits/misc.py
--- a/yt/frontends/fits/misc.py
+++ b/yt/frontends/fits/misc.py
@@ -171,14 +171,18 @@
nx = ds.domain_dimensions[ds.lon_axis]
ny = ds.domain_dimensions[ds.lat_axis]
mask = filter.mask((ny,nx)).transpose()
+ if ds.events_data:
+ prefix = "event_"
+ else:
+ prefix = ""
def _reg_field(field, data):
- i = data["xyz"[ds.lon_axis]].ndarray_view().astype("int")-1
- j = data["xyz"[ds.lat_axis]].ndarray_view().astype("int")-1
+ i = data[prefix+"xyz"[ds.lon_axis]].d.astype("int")-1
+ j = data[prefix+"xyz"[ds.lat_axis]].d.astype("int")-1
new_mask = mask[i,j]
- ret = data["zeros"].copy()
+ ret = np.zeros(data[prefix+"x"].shape)
ret[new_mask] = 1.
return ret
- ds.add_field(("gas",reg_name), sampling_type="cell", function=_reg_field)
+ ds.add_field(("gas", reg_name), sampling_type="cell", function=_reg_field)
if obj is None:
obj = ds.all_data()
if field_parameters is not None:
diff -r 49a03a144a9648a461207eb0942534f6762bea72 -r b9c09a2edd4e61caf2678093fcfdec8dceb0b20c yt/units/dimensions.py
--- a/yt/units/dimensions.py
+++ b/yt/units/dimensions.py
@@ -32,6 +32,8 @@
rate = 1 / time
frequency = rate
+solid_angle = angle * angle
+
velocity = length / time
acceleration = length / time**2
jerk = length / time**3
@@ -53,6 +55,8 @@
angular_momentum = mass*length*velocity
specific_angular_momentum = angular_momentum / mass
specific_energy = energy / mass
+count_flux = 1 / (area*time)
+count_intensity = count_flux / solid_angle
# Gaussian electromagnetic units
charge_cgs = (energy * length)**Rational(1, 2) # proper 1/2 power
@@ -77,8 +81,6 @@
resistance = resistance_cgs
current = current_cgs
-solid_angle = angle * angle
-
derived_dimensions = [rate, velocity, acceleration, jerk, snap, crackle, pop,
momentum, force, energy, power, charge_cgs, electric_field_cgs,
magnetic_field_cgs, solid_angle, flux, specific_flux, volume,
diff -r 49a03a144a9648a461207eb0942534f6762bea72 -r b9c09a2edd4e61caf2678093fcfdec8dceb0b20c yt/units/unit_lookup_table.py
--- a/yt/units/unit_lookup_table.py
+++ b/yt/units/unit_lookup_table.py
@@ -124,8 +124,8 @@
"me": (mass_electron_grams, dimensions.mass, 0.0, r"m_e"),
"mp": (mass_hydrogen_grams, dimensions.mass, 0.0, r"m_p"),
"mol": (1.0/amu_grams, dimensions.dimensionless, 0.0, r"\rm{mol}"),
- 'Sv': (cm_per_m**2, dimensions.specific_energy, 0.0,
- r"\rm{Sv}"),
+ 'Sv': (cm_per_m**2, dimensions.specific_energy, 0.0, r"\rm{Sv}"),
+ "rayleigh": (0.25e6/np.pi, dimensions.count_intensity, 0.0, r"\rm{R}"),
# for AstroPy compatibility
"solMass": (mass_sun_grams, dimensions.mass, 0.0, r"M_\odot"),
diff -r 49a03a144a9648a461207eb0942534f6762bea72 -r b9c09a2edd4e61caf2678093fcfdec8dceb0b20c yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -166,16 +166,26 @@
else:
self.set_wcs(wcs)
- def set_wcs(self, wcs):
+ def set_wcs(self, wcs, wcsname=None, suffix=None):
"""
Set the WCS coordinate information for all images
with a WCS object *wcs*.
"""
- self.wcs = wcs
- h = self.wcs.to_header()
+ if wcsname is None:
+ wcs.wcs.name = 'yt'
+ else:
+ wcs.wcs.name = wcsname
+ h = wcs.to_header()
+ if suffix is None:
+ self.wcs = wcs
+ else:
+ setattr(self, "wcs"+suffix, wcs)
for img in self.hdulist:
for k, v in h.items():
- img.header[k] = v
+ kk = k
+ if suffix is not None:
+ kk += suffix
+ img.header[kk] = v
def update_header(self, field, key, value):
"""
@@ -348,7 +358,9 @@
def create_sky_wcs(self, sky_center, sky_scale,
ctype=["RA---TAN","DEC--TAN"],
- crota=None, cd=None, pc=None):
+ crota=None, cd=None, pc=None,
+ wcsname="celestial",
+ replace_old_wcs=True):
"""
Takes a Cartesian WCS and converts it to one in a
celestial coordinate system.
@@ -369,6 +381,10 @@
Dimensioned coordinate transformation matrix.
pc : 2x2-element ndarray, optional
Coordinate transformation matrix.
+ replace_old_wcs : boolean, optional
+ Whether or not to overwrite the default WCS of the
+ FITSImageData instance. If false, a second WCS will
+ be added to the header. Default: True.
"""
old_wcs = self.wcs
naxis = old_wcs.naxis
@@ -402,7 +418,10 @@
new_wcs.wcs.cd = cd
if pc is not None:
new_wcs.wcs.cd = pc
- self.set_wcs(new_wcs)
+ if replace_old_wcs:
+ self.set_wcs(new_wcs, wcsname=wcsname)
+ else:
+ self.set_wcs(new_wcs, wcsname=wcsname, suffix="a")
class FITSImageBuffer(FITSImageData):
pass
diff -r 49a03a144a9648a461207eb0942534f6762bea72 -r b9c09a2edd4e61caf2678093fcfdec8dceb0b20c yt/utilities/tests/test_fits_image.py
--- a/yt/utilities/tests/test_fits_image.py
+++ b/yt/utilities/tests/test_fits_image.py
@@ -113,6 +113,16 @@
yield assert_equal, fid4.get_data("density"), fits_oap.get_data("density")
+ fid4.create_sky_wcs([30., 45.], (1.0, "arcsec/kpc"), replace_old_wcs=False)
+ assert fid4.wcs.wcs.cunit[0] == "cm"
+ assert fid4.wcs.wcs.cunit[1] == "cm"
+ assert fid4.wcs.wcs.ctype[0] == "linear"
+ assert fid4.wcs.wcs.ctype[1] == "linear"
+ assert fid4.wcsa.wcs.cunit[0] == "deg"
+ assert fid4.wcsa.wcs.cunit[1] == "deg"
+ assert fid4.wcsa.wcs.ctype[0] == "RA---TAN"
+ assert fid4.wcsa.wcs.ctype[1] == "DEC--TAN"
+
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"])
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