[yt-dev] Issues with 3D convolution and reshaping data?

[Nathan Goldbaum]

On Mon, Aug 8, 2016 at 2:49 PM, Rasmi Elasmar <re2300 at columbia.edu> wrote:

> Hi all,
>
> I'm trying to gaussian smooth a 3D field of gas. I'm convolving the field
> with a gaussian kernel, but my results look odd. I think this might be
> because of how I'm reshaping my data. You can see my scripts here
> <https://github.com/rasmi/deviant-cosmology/blob/master/scripts/convolve.py>
> and here
> <https://github.com/rasmi/deviant-cosmology/blob/4f6b4886c01c92cfe3987d4a66e25ba84ecebf86/scripts/convolve.py>,
> and a more detailed explanation with plots here
> <https://github.com/astropy/astropy/issues/5169#issuecomment-238304706>.
>
> My suspicion (and someone from astropy's) is that I'm reshaping data
> improperly. I access the fields like this:
> ad = ds.all_data()
> density = ad[fields['density']]
> Which gives me a 1D array of length (1024^3) = 1073741824. To convolve
> this with a 3D gaussian kernel, I attempt to reshape the 1D array into its
> original 3D shape:
> density = density.reshape(ds.domain_dimensions)
> Which gives me a 3D 1024^3 matrix, but I'm not sure if this is correct.
> Does this correctly transform the 1D data from ds.all_data()[field] back
> into its 3D grid, or are things being misaligned somehow? Is there a better
> way to access a certain field's data in 3D without getting it in 1D first
> and reshaping?
>

HI Ramsi,

As you've discovered, yt data objects return data as a flattened, unordered
1D array. Simply reshaping the data as you've done will not produce what
you're expect them to produce, since the data are unordered.

It sounds like you're dealing with a uniform resolution dataset, and what
I'm saying to you probably doesn't make much sense in that context, but
keep in mind that yt can deal with multiresolution data natively, where
this flattened, unordered representation is the most convenient way to deal
with the data.

For this particular use case, I'd bet that what you really want is either a
covering_grid data object or an arbitrary_grid data object:

http://yt-project.org/doc/examining/low_level_inspection.html#examining-grid-data-in-a-fixed-resolution-array
http://yt-project.org/doc/analyzing/objects.html#arbitrary-grid

http://yt-project.org/doc/reference/api/generated/yt.data_objects.construction_data_containers.YTArbitraryGrid.html
http://yt-project.org/doc/reference/api/generated/yt.data_objects.construction_data_containers.YTCoveringGrid.html

A covering grid is a uniform resolution 3D version of the data object. For
multiresolution data, this will simply sample the data at a single
resolution level, using nearest neighbor interpolation where the data are
too coarse, and sampling when the data are too finely resolved. Covering
grids are somewhat restrictive, as the resolution of the grid much be a
power-of-two multiple of the coarsest resolution of your data.

Arbitrary_grid lets you create uniform resolution data at an arbitrary
position or resolution, so it's a bit more flexible, although it's not
necessarily as "honest" a representation of the original dataset.

Hope that helped clear things up, please let us know if you have additional
questions.

-Nathan



>
>
> Thanks,
>
> Rasmi
>
> _______________________________________________
> yt-dev mailing list
> yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-dev-spacepope.org/attachments/20160808/9811f6fc/attachment.html>