Mailman 3 October 2016 - yt-users

Gadget velocity question
by Jared Coughlin 17 Oct '16

17 Oct '16

Hello! Just a quick question about the gadget velocities. For cosmological simulations, the default gadget snapshot does not actually contain the peculiar velocities v_pec, but, rather, the quantity u, where v_pec = u * sqrt(a), and a is the scale factor. My question is this: When I do: ds = yt.load("snapshot") and then access the velocity field, is yt giving me u or v_pec? That is, is yt doing the conversion from u to v_pec behind the scenes or not? I tried looking in the source code to find some conversion, but I didn't really know where to look. Thanks! -Jared

2 1

How to calculate a field THEN plot a slice, not have SlicePlot recalculate the field
by Rajika Kuruwita 17 Oct '16

17 Oct '16

Hi yt-users, I have created a field that calculates the enclosed mass in the simulation domain from a centre. This isn't a conventional field which can be calculated from one-to-one matrix multiplication like most fields, but needs to be aware of all the cells interior to the cell being calculated for. The field returns the correct values when I retrieve it using: dd = ds.all_data() dd['Enclosed_Mass'], But when creating a Slice plot it gives the wrong values. I realise this is because when creating a plot with one of those functions, yt calculates the field only for the cells used in the slice plot, hence severely underestimates the actual enclosed mass because it is ignoring all other cells. I know this is to make yt much faster and works for most fields, but I'm wondering if anyone else has encountered this issue, and figured out work arounds? Regards, Rajika -- Rajika Kuruwita PhD Candidate at Australian National University

2 1

Derived fields for generic array data
by Will Barnes 15 Oct '16

15 Oct '16

Apologies if this message shows up twice. I sent it from the wrong account the first time ---- Hi yt users, I have a question about derived fields in a dataset loaded from generic array data. I’ve created a dataset from a three dimensional Numpy array and traced some streamlines (say N of them) through the volume (this is a 3D magnetic field) such that I now have a set of N fieldlines. External to the dataset, I then have N Numpy arrays corresponding to each of the N streamlines, each the same length of the corresponding streamline, that represent the density as a function of the position along the streamline. What I’d like to do is create a new derived field (“density”) where the values from each of my external numpy arrays are mapped to the corresponding coordinate (as determined by the streamlines) and then “density” is zero everywhere else. I’ve read the documentation about creating derived fields and adding fields to a dataset, but it is not clear to me how to use external information (e.g. my set of N arrays) to create a new or derived field as all the examples show only how to use already existing fields in the dataset to create a new field. This seems like a relatively simple thing so perhaps I am just missing something or not quite understanding how a new field is added to a dataset? Since my explanation above is a bit confusing, I’ve included a code snippet below that hopefully illustrates a bit better what I’d like to do. Thanks! Will Barnes ``` # create 3d numpy array arr = np.random.random(size=(64,64,64)) # load generic data into a dataset data = dict(density = (arr, "g/cm**3")) bbox = np.array([[-1.5, 1.5], [-1.5, 1.5], [-1.5, 1.5]]) ds = yt.load_uniform_grid(data, arr.shape, length_unit="Mpc", bbox=bbox, nprocs=64) # set seed points c = ds.domain_center N = 100 scale = ds.domain_width[0] pos_dx = np.random.random((N,3))*scale-scale/2. pos = c+pos_dx # Create streamlines streamlines = Streamlines(ds, pos, 'velocity_x', 'velocity_y', 'velocity_z',length=1.0*Mpc, get_magnitude=True) streamlines.integrate_through_volume() # create a density field in ds initialized as 0 everywhere # for each streamline in streamlines, find coordinates from streamline and add the appropriate density value to the coordinates ```

1 0

Related to Off axis plot and Dark matter density field
by Prateek Gupta 14 Oct '16

14 Oct '16

Dear Yt users, I want to the off axis slice plot of density. But I didn't understand the use of north_vector. We already gave the normal vector to the plane and the point and these two things are enough to write a plane equation. I think I am confused with the concept. Please help me in understanding it very well. I am facing one more problem. I want the dark matter density Sliceplot. But when I used the field "dark_matter_density", it's not giving the correct plot I think. Here I attach the plot also for your reference. (check2.png) Then I searched and got that it can also be done if I plot the particle_density so I define field like this as given in yt-doc. def particle_density(field, data): pos = data[ptype, coord_name] mass = data[ptype, mass_name] pos.convert_to_units("code_length") mass.convert_to_units("code_mass") d = data.deposit(pos, [data[ptype, mass_name]], method = "sum") d = data.ds.arr(d, "code_mass") d /= data["index", "cell_volume"] return d but this also give an error says : *NameError: global name 'ptype' is not defined* *but "all_density" field is working fine. *But it will be total density (dark + baryon *)* Please help me regarding this issue. Thanks in advance. Prateek Gupta.

2 1

VTK how to Import 3D image data via imageImport on c++ for Volume Rendering
by Honghe 13 Oct '16

13 Oct '16

down votefavorite Thanks for VTK first, it is a nice Tool. But I meet a problem as follow. Through VTK Python API, I can import/make 3d image data and Volume Rendering as follow: And the Python code is: # coding=utf-8 import vtk from numpy import * import random num = 74 # We begin by creating the data we want to render. # For this tutorial, we create a 3D-image containing three overlaping cubes. # This data can of course easily be replaced by data from a medical CT-scan or anything else three dimensional. # The only limit is that the data must be reduced to unsigned 8 bit or 16 bit integers. data_matrix = zeros([num, num, num], dtype=uint8) # data_matrix[0:35, 0:35, 0:35] = 50 # data_matrix[25:55, 25:55, 25:55] = 100 # data_matrix[45:74, 45:74, 45:74] = 150 for i in range(0, num - 1): for j in range(0, num - 1): for k in range(0, num - 1): data_matrix[i:i + 1, j:j + 1, k:k + 1] = 100 # For VTK to be able to use the data, it must be stored as a VTK-image. This can be done by the vtkImageImport-class which # imports raw data and stores it. dataImporter = vtk.vtkImageImport() # The preaviusly created array is converted to a string of chars and imported. data_string = data_matrix.tostring() # dataImporter.CopyImportVoidPointer(data_string, len(data_string)) dataImporter.SetImportVoidPointer(data_matrix) # The type of the newly imported data is set to unsigned char (uint8) dataImporter.SetDataScalarTypeToUnsignedChar() # Because the data that is imported only contains an intensity value (it isnt RGB-coded or someting similar), the importer # must be told this is the case. dataImporter.SetNumberOfScalarComponents(1) # The following two functions describe how the data is stored and the dimensions of the array it is stored in. For this # simple case, all axes are of length 75 and begins with the first element. For other data, this is probably not the case. # I have to admit however, that I honestly dont know the difference between SetDataExtent() and SetWholeExtent() although # VTK complains if not both are used. dataImporter.SetDataExtent(0, 74, 0, 74, 0, 74) dataImporter.SetWholeExtent(0, 74, 0, 74, 0, 74) # The following class is used to store transparencyv-values for later retrival. In our case, we want the value 0 to be # completly opaque whereas the three different cubes are given different transperancy-values to show how it works. alphaChannelFunc = vtk.vtkPiecewiseFunction() alphaChannelFunc.AddPoint(0, 0.0) # alphaChannelFunc.AddPoint(50, 0.001) # alphaChannelFunc.AddPoint(100, 0.02) alphaChannelFunc.AddPoint(150, 0.02) # This class stores color data and can create color tables from a few color points. For this demo, we want the three cubes # to be of the colors red green and blue. colorFunc = vtk.vtkColorTransferFunction() colorFunc.AddRGBPoint(50, 0.0, 0.2, 0.2) colorFunc.AddRGBPoint(100, 0.0, 1.0, 0.0) colorFunc.AddRGBPoint(150, 0.0, 0.0, 1.0) # The preavius two classes stored properties. Because we want to apply these properties to the volume we want to render, # we have to store them in a class that stores volume prpoperties. volumeProperty = vtk.vtkVolumeProperty() volumeProperty.SetColor(colorFunc) volumeProperty.SetScalarOpacity(alphaChannelFunc) # This class describes how the volume is rendered (through ray tracing). compositeFunction = vtk.vtkVolumeRayCastCompositeFunction() # We can finally create our volume. We also have to specify the data for it, as well as how the data will be rendered. volumeMapper = vtk.vtkVolumeRayCastMapper() volumeMapper.SetVolumeRayCastFunction(compositeFunction) volumeMapper.SetInputConnection(dataImporter.GetOutputPort()) # The class vtkVolume is used to pair the preaviusly declared volume as well as the properties to be used when rendering that volume. volume = vtk.vtkVolume() volume.SetMapper(volumeMapper) volume.SetProperty(volumeProperty) # With almost everything else ready, its time to initialize the renderer and window, as well as creating a method for exiting the application renderer = vtk.vtkRenderer() renderWin = vtk.vtkRenderWindow() renderWin.AddRenderer(renderer) renderInteractor = vtk.vtkRenderWindowInteractor() renderInteractor.SetRenderWindow(renderWin) # We add the volume to the renderer ... renderer.AddVolume(volume) # ... set background color to white ... renderer.SetBackground(1, 1, 1) # ... and set window size. renderWin.SetSize(400, 400) # A simple function to be called when the user decides to quit the application. def exitCheck(obj, event): if obj.GetEventPending() != 0: obj.SetAbortRender(1) # Tell the application to use the function as an exit check. # renderWin.AddObserver("AbortCheckEvent", exitCheck) renderInteractor.Initialize() # Because nothing will be rendered without any input, we order the first render manually before control is handed over to the main-loop. renderWin.Render() renderInteractor.Start() But for c++ API, I can not get the same Volume Rendering through ImageImport. The following code is the C++ code I write refer to the above Python code. But the result is as follow： the C++ code I write refer to the above Python code: #include <iostream> #include <vtkStructuredPointsReader.h> #include <vtkVolumeRayCastCompositeFunction.h> #include <vtkVolumeRayCastMapper.h> #include <vtkRendererCollection.h> #include <vtkSmartPointer.h> #include <vtkImageImport.h> #include <stdlib.h> // VTK includes #include "vtkBoxWidget.h" #include "vtkCamera.h" #include "vtkCommand.h" #include "vtkColorTransferFunction.h" #include "vtkDICOMImageReader.h" #include "vtkImageData.h" #include "vtkImageResample.h" #include "vtkMetaImageReader.h" #include "vtkPiecewiseFunction.h" #include "vtkPlanes.h" #include "vtkProperty.h" #include "vtkRenderer.h" #include "vtkRenderWindow.h" #include "vtkRenderWindowInteractor.h" #include "vtkVolume.h" #include "vtkVolumeProperty.h" #include "vtkXMLImageDataReader.h" #include "vtkFixedPointVolumeRayCastMapper.h" int main() { std::cout << "Hello, World!" << std::endl; // create data // Create a c-style image const int width = 4; const int height = 4; const int depth = 4; unsigned char cImage[width * height * depth]; unsigned char value = 0; for (unsigned int row = 0; row < height; ++row) { for (unsigned int col = 0; col < width; ++col) { for (unsigned int z = 0; z < depth; ++z) { cImage[row * width * depth + col * depth + z] = 100; } } } // Convert the c-style image to a vtkImageData vtkSmartPointer<vtkImageImport> imageImport = vtkSmartPointer<vtkImageImport>::New(); imageImport->SetImportVoidPointer(cImage); imageImport->SetDataScalarTypeToUnsignedChar(); imageImport->SetNumberOfScalarComponents(1); imageImport->SetDataSpacing(1, 1, 1); imageImport->SetDataOrigin(0, 0, 0); imageImport->SetDataExtent(0, width - 1, 0, height - 1, 0, depth - 1); imageImport->SetWholeExtent(0, width - 1, 0, height - 1, 0, depth - 1); // imageImport->SetDataExtentToWholeExtent(); // imageImport->SetDataScalarTypeToUnsignedChar(); // imageImport->SetNumberOfScalarComponents(1); imageImport->Update(); // Create the standard ren, render window and interactor vtkRenderer *ren = vtkRenderer::New(); vtkRenderWindow *renWin = vtkRenderWindow::New(); renWin->AddRenderer(ren); vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New(); iren->SetRenderWindow(renWin); // Create the reader for the data // vtkStructuredPointsReader *reader = vtkStructuredPointsReader::New(); // reader->SetFileName("/home/honhe/ClionProjects/Demo_vtk_SimpleRayCast/ironProt.vtk"); // Create transfer mapping scalar value to opacity vtkPiecewiseFunction *opacityTransferFunction = vtkPiecewiseFunction::New(); opacityTransferFunction->AddPoint(0, 0.0); opacityTransferFunction->AddPoint(150, 0.2); // Create transfer mapping scalar value to color vtkColorTransferFunction *colorTransferFunction = vtkColorTransferFunction::New(); colorTransferFunction->AddRGBPoint(0.0, 0.0, 0.0, 0.0); colorTransferFunction->AddRGBPoint(100.0, 0.0, 1.0, 0.0); // The property describes how the data will look vtkVolumeProperty *volumeProperty = vtkVolumeProperty::New(); volumeProperty->SetColor(colorTransferFunction); volumeProperty->SetScalarOpacity(opacityTransferFunction); volumeProperty->ShadeOn(); volumeProperty->SetInterpolationTypeToLinear(); // The mapper / ray cast function know how to render the data vtkVolumeRayCastCompositeFunction *compositeFunction = vtkVolumeRayCastCompositeFunction::New(); vtkVolumeRayCastMapper *volumeMapper = vtkVolumeRayCastMapper::New(); volumeMapper->SetVolumeRayCastFunction(compositeFunction); volumeMapper->SetInputConnection(imageImport->GetOutputPort()); // The volume holds the mapper and the property and // can be used to position/orient the volume vtkVolume *volume = vtkVolume::New(); volume->SetMapper(volumeMapper); volume->SetProperty(volumeProperty); ren->AddVolume(volume); ren->SetBackground(1, 1, 1); renWin->SetSize(400, 400); iren->Initialize(); renWin->Render(); iren->Start(); return 0; } And I search example code, but only find the following, they can not work with 3d image data. This example show how to import 2d Image VTK/Examples/Cxx/Images/ImageImport And this example show how to deal with 3d image, but the API is out of date This example demonstrates how to set and access locations in a 3D image. This problem also post in http://stackoverflow.com/questions/39976669/vtk-how-to-import-3d-image-data…, but no responce right now. Any help appreciate! ------------------ Honghe

2 1

pyXSIM version 1.1.0
by John Zuhone 11 Oct '16

11 Oct '16

Hi everyone, pyXSIM has been updated to version 1.1.0. This version contains a bugfix and some minor new features. Fixed a bug which did not use the correct file names for AtomDB tables when using TableApecModel. Refactored the absorption model handling into a new class. No user-facing changes have been made. Added special classes for the TBabs and wabs absorption models. De-emphasizing XSpec-based spectral models in favor of the table-based alternatives. For more information see the website: http://hea-www.cfa.harvard.edu/~jzuhone/pyxsim <http://hea-www.cfa.harvard.edu/~jzuhone/pyxsim> pyXSIM can be installed in a few different ways. The simplest way is via the conda package if you have the Anaconda Python Distribution: [~]$ conda install -c jzuhone pyxsim The second way to install pyXSIM is via pip. pip will attempt to download the dependencies and install them, if they are not already installed in your Python distribution: [~]$ pip install pyxsim Alternatively, to install into your Python distribution from source (http://github.com/jzuhone/pyxsim <http://github.com/jzuhone/pyxsim>): [~]$ python setup.py install Thanks for using pyXSIM, John ZuHone

1 0

Projection questions
by Scott Feister 04 Oct '16

04 Oct '16

Hi all, 1) Is there a way to do a "mean" ProjectionPlot? As far as I can tell, the options are "integrate", "mip", and "sum", so I am guessing the answer is no. I know you can do this with regions by reg.mean(), and make your own plot in matplotlib. 2) If I do a Z projection on a domain with high X,Y spatial resolution, will it average cells along X and Y (within the pixel limits) at each Z depth as it projects? Or, just pick a single cell at X,Y for each pixel? Also, does anyone know if the reg.mean() function weights the cells equally in the mean if there is varied spatial resolution? For example, I know the reg.integrate() function multiplies by path length, but reg.sum() simply adds up cell values. Thanks! Scott Scott Feister, Ph.D. Postdoctoral Researcher, Flash Center for Computational Science University of Chicago, Department of Astronomy and Astrophysics

2 3

Re: [yt-users] Cut Region Failure
by Nathan Goldbaum 01 Oct '16

01 Oct '16

Hi Jason, Unfortunately no. I think I might be able to come back to this issue soonish, we're going to be introducing a new fast periodic KDTree which could be exploited here. For now I think you're going to need to make this plot in a different way, without using a cut_region. Perhaps by overplotting contours? Nathan On Saturday, October 1, 2016, Jason Galyardt <jason.galyardt(a)gmail.com> wrote: > Hi Nathan, > > Is there a previous version of yt in which compound cut regions work? > > Thanks, > Jason > > ----- > > Hi Jason, > > This is a known issue: > > https://bitbucket.org/yt_analysis/yt/issues/941/cutregionsel > ector-doesnt-work-with-some > > Unfortunately fixing it will require a modification to the cut_region > implementation to do a better job of searching for cells that belong to the > region for chained selectors. > > -Nathan > > On Fri, Sep 30, 2016 at 1:25 PM, Jason Galyardt <jason.galyardt(a)gmail.com > <javascript:_e(%7B%7D,'cvml','jason.galyardt(a)gmail.com');>> wrote: > Dear all, > > I would like to generate a compound cut region for a FLASH simulation, but > the output I'm getting is odd. I've pasted my simple script here: > > http://paste.yt-project.org/show/6839/ > > It loads a simulation file, creates a compound cut_region (1e4 K < T < 1e6 > K), and applies it to a slice plot. The resulting image is banded, as you > can see by the attached plots. Am I missing something, or is this a problem > related to yt? > > FYI, I get the same behavior when I use a projection plot or only one of > the conditionals (e.g. T < 1e6 K). I've tried both the stock 3.3.1 version, > as well as the self-contained development version that I installed earlier > this afternoon using install_script.sh. Let me know if any other info is > useful. > > Cheers, > Jason > >

1 0