# Iterating through every raster cell and calculating distance to other cells using Python?

I have a binary raster (.tif) in ArcGIS that shows builtup areas (value=1) and non builtup areas (value=0) like this: In order to make a proximity analysis, I need to iterate through the raster with arcpy/python and for each cell (i) get the distance to EVERY other cell (j) with value=1 within in a given radius around the current cell (i) and write all those distance values to a list/array.

Is there a performant way to do this using Python?

I know there is the EucDistance tool in ArcGIS: but it may not be the smartest, to create several thousand new rasters.

My plan was to use a RasterToNumpyArray and then calculate on the numpy array

So far I tried the following Python code:

``````import numpy
from scipy.spatial import distance

def dist_euclid(cell_a, cell_b):
dist = distance.euclidean(cell_a, cell_b)
return dist

'''Masking is important to speed up processes'''
a, b = coords, coords
nx, ny = array.shape
y, x = numpy.ogrid[-a:nx - a, -b:ny - b]
# create mask where each value in radius around coords is true and everything else is false

# reshape mask to rectangular shape
indices = numpy.meshgrid(numpy.arange(min(i), max(i) + 1),
numpy.arange(min(j), max(j) + 1),
indexing='ij')

def calcCellDistancesForEachCellInRaster(in_array, buffer_size):
'''For each cell make slice sub-array around the cell and calculate distances to every other cell (value=1)'''
dist_arr = []
for (source_x, source_y), source_value in numpy.ndenumerate(in_array):
if source_value == 1:
for (target_x, target_y), target_value in numpy.ndenumerate(masked_array):
if target_value == 1 and ([source_x, source_y] != [target_x, target_y]):
dist = dist_euclid([source_x, source_y], [target_x, target_y])
if dist <= buffer_size:
dist_arr.append(dist)
# print("dist from %s to %s = %s" % ([source_x, source_y], [target_x, target_y], dist))
else:
pass
else:
pass
return dist_arr

# Create a simple array from scratch using random values
myArray = numpy.random.randint(2, size=(100, 100))
myArray.shape = (100, 100)

dist_arr = calcCellDistancesForEachCellInRaster(myArray, 2)
``````
• For a question about performance of Python code I think you will also need to describe the hardware and OS configuration on which you are attaining that performance. – PolyGeo May 4 '19 at 22:39
• What was the difference in performance between your EucDistance and other Python code? Did they produce identical results from identical inputs? – PolyGeo May 6 '19 at 11:19
• The arcpy EucDistance does not iterate over every cell of the array and so far I did not find any other code that does so. Hence, I could not test other codes – feinheitsbrei May 6 '19 at 11:30
• An idea, convert your raster to a point dataset then call the generate near table tool. See how that performs? – Hornbydd May 10 '19 at 12:56
• What kind of output you need? Is it a raster of minimum distance extracted from dist_arr? Does buffer vary over area of interest? – FelixIP Jun 5 '19 at 22:10

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):

circle_index = np.zeros(array_shape) * np.nan
circle_distance = np.zeros(array_shape) * np.nan

nx,ny = circle_index.shape
y,x = np.ogrid[-a:nx-a,-b:ny-b]

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.

# Create circles.

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]

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.

• Your code does not "make" the circles around each pixel of the source array (my_array), right? In my case I would have to do the for-loop for every pixel, or am I taking this the wrong way? – feinheitsbrei May 20 '19 at 20:38
• well it does (kindoff), the "nd.filters.convolve" puts the circle on top of each pixel. But this is very fast and more or less independent of your size of my_array. What does have an influence on the speed (and on the for-loop) is the size of the circle. I.e. the more pixels are inside the circle, the longer the for-loop will be be and the larger the distances array will become. – Bert Coerver May 21 '19 at 7:08
• Ok, it works very fast but (as you said) it runs me into a Memory Error, when testing it with a realistic raster – feinheitsbrei May 21 '19 at 18:41
• well just to be clear, the memory error is probably cause by your radius (so not your "realistic raster"). What you could try to do is completely removing the " distances" array, and instead saving the "distance_map" variable to a Gtiff or netcdf or something. Then later you can do your analysis on them, depending on what you want to do with them. – Bert Coerver May 22 '19 at 7:24
• In my case it is already the np.zeros (distances = np.zeros((distances_per_pixel, x, y))) that is causing the error. Because x and y have to be the same as the source raster dimensions right? – feinheitsbrei May 22 '19 at 12:59