Hi,
I'm trying to rebin some uniform gridded data to polar coordinates. The original data (Z) is noisy but on a fine grid. ### snip # define grid in cartesian coordinates x = arange(-500, 500) y = arange(-500, 500) X, Y = meshgrid(x, y) Z = X**2 + Y**2 Z += uniform(-5000, 5000, size=Z.shape)) # transform to polar coordinates def cart2pol(x, y): theta = arctan2(y, x) rho = sqrt(x**2 + y**2) return (theta, rho) THETA_xy, RHO_xy = cart2pol(X, Y) # define new polar grid theta = arange(-pi, pi, 0.1) rho = arange(0.1 ,500, 2.0) THETA, RHO = meshgrid(theta, rho) # Tried to use griddata, but it gives of course very jumpy results, because it # does not take the average over the pixels in a patch (dTHETA, dRHO) but an interpolation at the # exact new point. interpolate.griddata((THETA_xy.ravel(), RHO_xy.ravel()), Z.ravel(), (THETA, RHO), method='linear') ### snip Is there a method that rebins the data to a new grid, possibly taking the average of the Z-values inside a patch (dTHETA, dRHO)? Thanks for any hints! |
On Tue, 28 Aug 2012 03:33:11 -0700 (PDT)
arsbbr <[hidden email]> wrote: > > Hi, > > I'm trying to rebin some uniform gridded data to polar coordinates. The > original data (Z) is noisy but on a fine grid. > > ### snip > > # define grid in cartesian coordinates > x = arange(-500, 500) > y = arange(-500, 500) > X, Y = meshgrid(x, y) > Z = X**2 + Y**2 > Z += uniform(-5000, 5000, size=Z.shape)) > > # transform to polar coordinates > > def cart2pol(x, y): > theta = arctan2(y, x) > rho = sqrt(x**2 + y**2) > return (theta, rho) > > THETA_xy, RHO_xy = cart2pol(X, Y) > > # define new polar grid > theta = arange(-pi, pi, 0.1) > rho = arange(0.1 ,500, 2.0) > > THETA, RHO = meshgrid(theta, rho) > > # Tried to use griddata, but it gives of course very jumpy results, because > it > # does not take the average over the pixels in a patch (dTHETA, dRHO) but an > interpolation at the > # exact new point. > > interpolate.griddata((THETA_xy.ravel(), RHO_xy.ravel()), > Z.ravel(), (THETA, RHO), method='linear') > > ### snip > > Is there a method that rebins the data to a new grid, possibly taking the > average > of the Z-values inside a patch (dTHETA, dRHO)? Try to use 2D-weighted histogram (divided by unweighted histograms). This unfortunately does dot split pixels over various bins so it works badly when the number of output bins is large. w=histogram2d(THETA_xy.ravel(),RHO_xy.ravel(),(int(2*pi/0.1),500//2),weights=Z.ravel()) u=histogram2d(THETA_xy.ravel(),RHO_xy.ravel(),(int(2*pi/0.1),500//2)) imshow (w[0]/u[0]) If you want to split pixels, I am working on more sophisticated algorithm to do that in https://github.com/kif/pyFAI Tell me if you are interested. Cheers, -- Jérôme Kieffer On-Line Data analysis / Software Group ISDD / ESRF tel +33 476 882 445 _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
Dear Jerome,
thank you for your tip! Unfortunately I'm getting big gaps for my real data in the small rho region. The optimal algorithm should interpolate in the small rho region and bin in the high rho region. I guess your split pixel solution would be best. But this a whole API with a lot of corrections and methods for real experimental data. Is it possible to use your sophisticated algorithm only for the transformation from Cartesian to a evenly stepped polar grid? Speed is not the issue with my data, so a numpy solution would be fine. Thanks a lot!
|
On Fri, 31 Aug 2012 00:36:21 -0700 (PDT)
arsbbr <[hidden email]> wrote: > thank you for your tip! Unfortunately I'm getting big gaps for my real data > in the small rho region. The optimal algorithm should interpolate in the > small rho region and bin in the high rho region. > I guess your split pixel solution would be best. But this a whole API with a > lot of corrections and methods for real experimental data. Is it possible to > use your sophisticated algorithm only for the transformation from Cartesian > to a evenly stepped polar grid? Speed is not the issue with my data, so a > numpy solution would be fine. I wish this would be possible using only numpy but I failed. I had to write it in a low-level (Cython) way to get decent speed because the algorithm is iterating over all input pixels. What you are looking for is in: https://github.com/kif/pyFAI/blob/master/src/splitBBox.pyx Use the function histoBBox2d(weights, pos0, delta_pos0, pos1, delta_pos1, bins=(100, 36)) If you refer to the manual of numpy.histogram2d: -> bins has the same meaning, i.e. the shape of the output image -> pos0 and pos1 are the same a x and y; they should contain the coordinated in (rho,theta) of the center of input pixels. -> delta_pos0 and delta_pos1 are the distance in (rho,theta) from pixel center to the pixel corner. This is of course new compared to numpy -> weights is as in numpy the intensity of input pixels. All array need to have the same size. The output is: -> your transformed image (maybe it's transposed), i.e. weighted_histogram/unweighted_histogram -> axis in rho (or theta) 1D -> axis in theta (or rho) 1D -> weighted histogram (2D) -> unweighted histogram (2D) splitBBox should be stand alone: just use this setup.py to compile the module (needs Cython) ######################################## try: from setuptools import setup except ImportError: from distutils.core import setup from distutils.core import Extension from Cython.Distutils import build_ext # for numpy from numpy.distutils.misc_util import get_numpy_include_dirs split_ext = Extension("splitBBox", include_dirs=get_numpy_include_dirs(), sources=['splitBBox.pyx'], setup(name='histogram', version="0.3.0", author="Jerome Kieffer", author_email="[hidden email]", description='test for azim int', ext_modules=[split_ext], cmdclass={'build_ext': build_ext}, ) ##################################################### Nota: I am working on X-ray diffraction so Your "rho" could be named "2theta", tth or q. Your "theta" is usually named "chi" in my comments. Different sciences; different conventions :) I hope this helps, cheers, -- Jérôme Kieffer On-Line Data analysis / Software Group ISDD / ESRF tel +33 476 882 445 _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
Free forum by Nabble | Edit this page |