[Numpy-discussion] Rasterizing points onto an array

[Nathan Bell]

On Sun, May 31, 2009 at 8:39 PM, Thomas Robitaille
<thomas.robitaille at gmail.com> wrote:
> Hi,
>
> I have a set of n points with real coordinates between 0 and 1, given
> by two numpy arrays x and y, with a value at each point represented by
> a third array z. I am trying to then rasterize the points onto a grid
> of size npix*npix. So I can start by converting x and y to integer
> pixel coordinates ix and iy. But my question is, is there an efficient
> way to add z[i] to the pixel given by (xi[i],yi[i])? Below is what I
> am doing at the moment, but the for loop becomes very inefficient for
> large n. I would imagine that there is a way to do this without using
> a loop?
>
> ---
>
> import numpy as np
>
> n = 10000000
> x = np.random.random(n)
> y = np.random.random(n)
> z = np.random.random(n)
>
> npix = 100
> ix = np.array(x*float(npix),int)
> iy = np.array(y*float(npix),int)
>
> image = np.zeros((npix,npix))
> for i in range(len(ix)):
>     image[ix[i],iy[i]] = image[ix[i],iy[i]] + z[i]
>

There might be a better way, but histogram2d() seems like a good fit:

import numpy as np

n = 1000000
x = np.random.random(n)
y = np.random.random(n)
z = np.random.random(n)

npix = 100

bins = np.linspace(0, 1.0, npix + 1)
image = np.histogram2d(x, y, bins=bins, weights=z)[0]

-- 
Nathan Bell wnbell at gmail.com
http://www.wnbell.com/