How to find densest point in subvolume of Gadget data
Hey everyone, I am trying to find the densest point in Gadget data in the subvolume with left edge [700]*3 and right edge [800]*3 where 700 and 800 are in code units. This is what I do [1], which was inspired by this page [2], except I try and take the all_data confined to a region. The output [3] says however where you will note the center is not inside the region. It found the densest point in the whole data set, not the region. So how would I do this instead for a subvolume of Gadget data? Thanks. [1] http://pastebin.com/kyUgEi3a [2] http://yt-project.org/doc/cookbook/gadget_notebook.html [3] http://pastebin.com/E4aLb5VC -- ------------------------------------------------------------------------ Joseph Smidt <josephsmidt@gmail.com> Theoretical Division P.O. Box 1663, Mail Stop B283 Los Alamos, NM 87545 Office: 505-665-9752 Fax: 505-667-1931 _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Hey Joseph, I don't know if this will work, but maybe try (sorry about the weird indenting...trying to type while also watching the kids in bath) 1. ds = load('myfile.hdf5',over_refine_factor=1) 2. bbox_lim = 50 center = whatever_your_desired_center bbox1 = [[center[0]-bbox_lim,center[0]+bbox_lim], [center[1]-bbox_lim,center[1]+bbox_lim], [center[2]-bbox_lim,center[2]+bbox_lim]] 3. ds1 = load(fname,bounding_box=bbox1,n_ref = 64 ,over_refine_factor=1)
ad1 = ds1.all_data() and search within ad1?
you might need to force periodicity off
ds1.periodicity = (False,False,False)
but maybe not. Anyways, I hope this works! -d On Sat, May 2, 2015 at 6:27 PM, Joseph Smidt <josephsmidt@gmail.com> wrote:
Hey everyone,
I am trying to find the densest point in Gadget data in the subvolume with left edge [700]*3 and right edge [800]*3 where 700 and 800 are in code units.
This is what I do [1], which was inspired by this page [2], except I try and take the all_data confined to a region. The output [3] says however where you will note the center is not inside the region. It found the densest point in the whole data set, not the region.
So how would I do this instead for a subvolume of Gadget data? Thanks.
[1] http://pastebin.com/kyUgEi3a
[2] http://yt-project.org/doc/cookbook/gadget_notebook.html
[3] http://pastebin.com/E4aLb5VC
-- ------------------------------------------------------------------------ Joseph Smidt <josephsmidt@gmail.com>
Theoretical Division P.O. Box 1663, Mail Stop B283 Los Alamos, NM 87545 Office: 505-665-9752 Fax: 505-667-1931 _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Hi all, Joe decided (offlist) to directly access the region object itself. Because re.ds references the dataset that the region came from, instead of using 1. ds=load('myfile.hdf5',over_refine_factor=1) 2. re=ds.region(center=[750]*3,left_edge=[700]*3,right_edge=[800]*3) 3. ad=re.ds.all_data() 4. density=ad[("PartType0","density")] 5. wdens=np.where(density==np.max(density)) 6. coordinates=ad[("PartType0","Coordinates")] 7. cen=coordinates[wdens][0] the question can be resolved with 1. ds=load('myfile.hdf5',over_refine_factor=1) 2. re=ds.region(center=[750]*3,left_edge=[700]*3,right_edge=[800]*3) 3. density=re[("PartType0","density")] 4. wdens=np.where(density==np.max(density)) 5. coordinates=re[("PartType0","Coordinates")] 6. cen=coordinates[wdens][0] Cheers, Aaron -- Aaron Smith NSF Graduate Research Fellow Department of Astronomy University of Texas at Austin www.as.utexas.edu/~asmith On 5/2/15 5:27 PM, Joseph Smidt wrote:
Hey everyone,
I am trying to find the densest point in Gadget data in the subvolume with left edge [700]*3 and right edge [800]*3 where 700 and 800 are in code units.
This is what I do [1], which was inspired by this page [2], except I try and take the all_data confined to a region. The output [3] says however where you will note the center is not inside the region. It found the densest point in the whole data set, not the region.
So how would I do this instead for a subvolume of Gadget data? Thanks.
[1] http://pastebin.com/kyUgEi3a
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Ah yes, the "ds" attribute of a data object is a reference to the dataset itself. Sorry that wasn't clear initially. On Saturday, May 2, 2015, Aaron Smith <asmith@astro.as.utexas.edu> wrote:
Hi all,
Joe decided (offlist) to directly access the region object itself. Because re.ds references the dataset that the region came from, instead of using
1. ds = load('myfile.hdf5',over_refine_factor=1) 2. re = ds.region(center=[750]*3, left_edge=[700]*3, right_edge=[800]* 3) 3. ad = re.ds.all_data() 4. density = ad[("PartType0","density")] 5. wdens = np.where(density == np.max(density)) 6. coordinates = ad[("PartType0","Coordinates")] 7. cen = coordinates[wdens][0]
the question can be resolved with
1. ds = load('myfile.hdf5',over_refine_factor=1) 2. re = ds.region(center=[750]*3, left_edge=[700]*3, right_edge=[800]* 3) 3. density = re[("PartType0","density")] 4. wdens = np.where(density == np.max(density)) 5. coordinates = re[("PartType0","Coordinates")] 6. cen = coordinates[wdens][0]
Cheers, Aaron
-- Aaron Smith NSF Graduate Research Fellow Department of Astronomy University of Texas at Austinwww.as.utexas.edu/~asmith
On 5/2/15 5:27 PM, Joseph Smidt wrote:
Hey everyone,
I am trying to find the densest point in Gadget data in the subvolume with left edge [700]*3 and right edge [800]*3 where 700 and 800 are in code units.
This is what I do [1], which was inspired by this page [2], except I try and take the all_data confined to a region. The output [3] says however where you will note the center is not inside the region. It found the densest point in the whole data set, not the region.
So how would I do this instead for a subvolume of Gadget data? Thanks.
[1] http://pastebin.com/kyUgEi3a
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
participants (4)
-
Aaron Smith -
Desika Narayanan -
Joseph Smidt -
Nathan Goldbaum