Hi,
I have a question about volume rendering. I'd like to run camera.allsky_projection on a data object with a particular shape, in this case a sphere with a cylindrical hole in the center. I was able to define such an object using the boolean operator "NOT". But feeding it to camera.allsky_projection leads to an error, with the traceback below.
Traceback (most recent call last): File "sky_maps.py", line 70, in <module> inner_radius = rr, rotation = rotation) File "/Users/mkjoung/work/yt/yt-i386/src/yt-hg/yt/visualization/volume_rendering/camera.py", line 1526, in allsky_projection dx = min(g.dds.min() for g in pf.h.find_point(center)[0])
AttributeError: 'AMRBooleanRegion' object has no attribute 'h'
(allsky_projection runs successfully when I input 'pf' as an argument.) Is there a workaround so that allsky_projection will accept the boolean object as pf-like?
Thanks, -Ryan
Hi Ryan,
I have a question about volume rendering. I'd like to run camera.allsky_projection on a data object with a particular shape, in this case a sphere with a cylindrical hole in the center. I was able to define such an object using the boolean operator "NOT". But feeding it to camera.allsky_projection leads to an error, with the traceback below.
I think I might have an idea on how to address this, but since I'd like to A) test the fix myself before pushing it and B) have a better understanding of what you're doing in the first place, would you mind sharing the script you're using? That would be most helpful, thanks!
-- Stephen Skory s@skory.us http://stephenskory.com/... 510.621.3687 (google voice)
Hi Stephen,
Sure, here it is. (I should give credit to Matt for a large chunk of the script.)
Thanks! -Ryan
from yt.mods import * from yt.utilities.math_utils import ortho_find import yt.visualization.volume_rendering.camera as camera from yt.utilities.lib import healpix_aitoff_proj, pixelize_healpix from itertools import izip, chain, repeat import h5py
pf = load("~/enzo/data/cosmo/r0058/redshift0058")
com = na.array([ 0.40329299, 0.47176085, 0.46132502]) L = na.array([ 0.25527702, -0.79298414, -0.55318153])
# We now have our angular momentum vector.
v1, v2, v3 = ortho_find(L) # v1 = L, v2 & v3 = are two vectors orthogonal to L
c = v2 * 8.0/pf['kpc'] + com
print "L =", L print "v1, v2, v3 =", v1, v2, v3
sp = pf.h.sphere(c, (1.0, 'kpc')) stars = sp["particle_type"] == 2 xv, yv, zv = [sp["particlevelocity%s" % ax] for ax in 'xyz'] bv = na.array([na.average(xv), na.average(yv), na.average(zv)]) rotation = na.array([v2,v3,v1])
sphvir = pf.h.sphere(c, (100.0, 'kpc'))
rdisk = 18./pf['kpc']
hdisk = 4./pf['kpc']
disk = pf.h.disk(c, v1, rdisk, hdisk)
halo = pf.h.boolean([sphvir, "NOT", disk]) # exclude the disk from the sphere
for g in pf.h.grids: g.set_field_parameter("center", c) g.set_field_parameter("bulk_velocity", bv) g.set_field_parameter("normal",np.array([0,0,1]))
NP = 256 resolution = 512 rr = 10 ntheta = 720 nphi = 360
field_names = []
for vi in na.mgrid[-250:250:10]:
fname = "VelocityTransmittedHI%+03i%+03i" % (vi, vi+10)
print fname
def _get(vel):
@derived_field(name=fname)
def VelocityTransmittedHI(field, data):
return data["HI_Density"] (
(((vel 1e5) < data["RadialVelocity"])
&(data["RadialVelocity"] < (vel + 10) *
1e5))).astype("float64")
field_names.append(fname)
_get(vi)
for field in field_names:
f = h5py.File("hi_maps_pixel.h5", "a")
if field in f["/"]: continue
allsky = camera.allsky_projection(halo, c, #pf, c,
100.0/pf['kpc'], NP, field,
inner_radius = rr, rotation = rotation)
print "Computed", field, allsky.min(), allsky.max()
img, count = pixelize_healpix(NP, allsky, ntheta, nphi, na.eye(3))
print "Saving", field, img.min(), img.max()
f.create_dataset("%sPixel" % field, data=img)
f.flush()
f.close()
On Oct 23, 2012, at 9:28 AM, Stephen Skory wrote:
Hi Ryan,
I have a question about volume rendering. I'd like to run camera.allsky_projection on a data object with a particular shape, in this case a sphere with a cylindrical hole in the center. I was able to define such an object using the boolean operator "NOT". But feeding it to camera.allsky_projection leads to an error, with the traceback below.
I think I might have an idea on how to address this, but since I'd like to A) test the fix myself before pushing it and B) have a better understanding of what you're doing in the first place, would you mind sharing the script you're using? That would be most helpful, thanks!
-- Stephen Skory s@skory.us http://stephenskory.com/... 510.621.3687 (google voice)
yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org...
Hi Ryan, Stephen,
The main issue here is that the cameras are explicitly fed a parameter file, not a data object. For some of the cameras, this was done because we don't currently (though it is in the works) have a brick homogenization routine that acts on a data object instead of a parameter file. For others (off-axis and allsky in particular), there would be a way to get around this, but it will require modifying the source.
Around line 1530 of yt/visualization/volume_rendering/camera.py, there are the following lines:
grids = pf.h.sphere(center, radius)._grids
sampler = ProjectionSampler(positions, vs, center, (0.0, 0.0, 0.0, 0.0),
image, uv, uv, np.zeros(3, dtype='float64'))
pb = get_pbar("Sampling ", len(grids))
for i,grid in enumerate(grids):
data = [grid[field] * grid.child_mask.astype('float64')
for field in fields]
I think what you would want to do is add a source= keyword argument to the allsky projection call, and then use the object to get the grids in question. How to get the boolean regions to work correctly, which might need modification of the above loop, I'll leave up to Stephen.
Ryan, If you get this to work, please feel free to fork the repository and issue a pull request with the new feature.
Cheers, Sam
On Tue, Oct 23, 2012 at 7:28 AM, Stephen Skory s@skory.us... wrote:
Hi Ryan,
I have a question about volume rendering. I'd like to run camera.allsky_projection on a data object with a particular shape, in this case a sphere with a cylindrical hole in the center. I was able to define such an object using the boolean operator "NOT". But feeding it to camera.allsky_projection leads to an error, with the traceback below.
I think I might have an idea on how to address this, but since I'd like to A) test the fix myself before pushing it and B) have a better understanding of what you're doing in the first place, would you mind sharing the script you're using? That would be most helpful, thanks!
-- Stephen Skory s@skory.us http://stephenskory.com/... 510.621.3687 (google voice)
yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org...
Hi Ryan and Sam,
thanks for the suggestion, Sam. I'll take a look at this in a bit (in and out today; guests in town).
-- Stephen Skory s@skory.us http://stephenskory.com/... 510.621.3687 (google voice)
Hi Ryan,
I think I've made a change that should allow you to now do what you want. You'll need to update your yt so you can use this latest changeset:
https://bitbucket.org/yt_analysis/yt/pull-request/314/modifying-allsky-proje...
All you should need to do is to change the line in your script to look like:
allsky = camera.allsky_projection(pf, c, 100.0/pf['kpc'], NP, field, inner_radius = rr, rotation = rotation, source = halo)
Good luck, and let us know if something doesn't work.
-- Stephen Skory s@skory.us http://stephenskory.com/... 510.621.3687 (google voice)
-
M Ryan Joung -
Sam Skillman -
Stephen Skory