I think this script does what you want, without iterating over each pixel in myArray:
import numpy as np
import scipy.ndimage as nd
import matplotlib.pyplot as plt
# Create a function that generates a circle (inside a square-shaped array) with
# different index values for each unique pixel and the corresponding distance to
# the center.
def make_circles(index):
array_shape = (radius*2+1, radius*2+1)
circle_index = np.zeros(array_shape) * np.nan
circle_distance = np.zeros(array_shape) * np.nan
a,b = (radius, radius)
nx,ny = circle_index.shape
y,x = np.ogrid[-a:nx-a,-b:ny-b]
mask = x*x + y*y <= radius*radius
mask[a,b] = False
pixels_in_mask = np.sum(mask)
circle_index[mask] = np.arange(0, pixels_in_mask, 1)
circle_distance[mask] = np.sqrt(x*x + y*y)[mask]
return circle_index, circle_distance
# Generate some dummy data
x = 100
y = 100
my_array = np.random.randint(2, size=(x, y))
# specify the radius in which to search for pixels in my_array with value = 1.
radius = 10
# Create circles.
circle_index, circle_distance = make_circles(radius)
plt.figure(1)
plt.title("Each pixel has a unique index number")
plt.imshow(circle_index)
plt.colorbar()
# Determine the maximum amount of pixels for which the distance
# needs to be calculated per pixel.
distances_per_pixel = np.int(np.nanmax(circle_index) + 1)
# Create an array in which the results will be stored.
distances = np.zeros((distances_per_pixel, x, y))
# Determine the distances from every pixel to the current relative pixel as
# defined by circle_index.
for position in circle_index[np.isfinite(circle_index)]:
distance = circle_distance[circle_index == position][0]
mask = np.fliplr(np.flipud(np.where(circle_index != position, 0, 1)))
distance_map = nd.filters.convolve(my_array, mask, mode = 'constant', cval = 0.0) * distance
distance_map[distance_map == 0.] = np.nan
distances[np.int(position), :, :] = distance_map
# Each layer in the distances array corresponds to the "unique index number" as shown
# in figure 1. e.g. distances[3,:,:] shows a map that gives the distance to the circle_pixel
# with number 3, IF that pixel has value 1 in my_array.
# Calculate some other things, for example the shortest distance within a pixels
# radius to a pixel with my_array = 1.
shortest_distances = np.nanmin(distances, axis = 0)
longest_distances = np.nanmax(distances, axis = 0)
mean_distances = np.nanmean(distances, axis = 0)
# Plot a map showing the shortest distances.
plt.figure(2)
plt.imshow(shortest_distances, vmin = 0, vmax = radius)
plt.colorbar()
Be carefull though, as you increase the radius of you search area around each pixel, the size of the "distances" array will grow fast. Requiring more and more RAM-memory.
Best Answer
Adding 'MERGE_ALG=ADD' to the list of options to gdal.RasterizeLayer() did the trick, thanks to kyle. An example:
This is documented in http://www.gdal.org/gdal__alg_8h.html.