# Find maximum of the difference between each cell of a raster and its surrounding cells

I would like to do the following for each cell (let's call it cell_i) of a raster:

1. Define a neighboring region around cell_i (for instance, a region defined by the 8 surrounding cells of cell_i).
2. For each cell in the region defined in (1), find the difference between its value and the value of cell_i.
3. For all the differences computed in (2), find the minimum.
4. Store the minimum value from (3) in a new raster in the same spatial location as cell_i.

Are there existing tools that can do this?

If not, can you provide pointers on how to script this (ideally in ArcPy)?

I think that I have the algorithm down (for instance, I would know how to implement this if the raster was a matrix in matlab or 2D array in C#).

What are the commands needed, such as iterating through the original raster and saving values in a new raster?

• Check the focal statistics tool in ArcGIS, it could solve your problem. – ami Dec 15 '16 at 9:45
• You can use GDAL and numpy modules in a python script (see my answer). – xunilk Dec 16 '16 at 0:09
• I think both suggestions are ok. I will build on the one that involves scripting because at some point I may want to change the type of calculation (i.e. instead of finding differences, do something else). – user57029 Dec 16 '16 at 5:38

This workaround can be done with GDAL python module. I used a small raster to make easier the corroboration of minimum values of the difference between each raster cell and its neighboring region around it. As the values in the periphery of raster has not 8 surrounding cells, next code add a "ring" of -999 values to facilitate the evaluation of differences (incorporated in 'new_list' list) in this area. This 'new_list' was used for determining minimum values before they be incorporated in 'raster' list. This 'raster' list was reshaped as one numpy array, with the same number of rows and columns that the original raster, and used for obtaining resulting raster.

``````from osgeo import gdal, osr
import numpy as np

driver = gdal.GetDriverByName('GTiff')
filename = "/home/zeito/pyqgis_data/test.tif"
dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

cols = dataset.RasterXSize
rows = dataset.RasterYSize

#inserting ring of -999 around data values
data = np.insert(data, 0, -999, axis = 0)
data = np.insert(data, 0, -999, axis = 1)
data = np.insert(data, len(data), -999, axis = 0)
data = np.insert(data, len(data), -999, axis = 1)

print data

raster = [ ]

for i in range(1, rows + 1):
for j in range(1, cols + 1):
val = data[i][j]
new_list = [ data[i-1][j-1], data[i-1][j], data[i-1][j+1],
data[i][j-1], data[i][j+1],
data[i+1][j-1], data[i+1][j], data[i+1][j+1]]

print val
print new_list

#list with differences
new_list = [ val - item for item in new_list if item != -999]

#determining minimum
raster.append(np.min(new_list))

#creating numpy array with minimum values by cell
raster = np.asarray(np.reshape(raster, (rows, cols)))

transform = dataset.GetGeoTransform()

# Create gtif file
driver = gdal.GetDriverByName("GTiff")

output_file = "/home/zeito/pyqgis_data/minimum_differences.tif"

dst_ds = driver.Create(output_file,
cols,
rows,
1,
gdal.GDT_Float32)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( raster )

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(transform)

wkt = dataset.GetProjection()

# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )

#Close output raster dataset
dataset = None
dst_ds = None
``````

Above code was run at the Python Console of QGIS and, there were printed the corroboration values (see next image). Resulting raster was obtained as expected (by using Value Tool plugin to corroborate cell raster values): • Thank you, this gives me a great overview of how to code this myself. – user57029 Dec 16 '16 at 5:37

I agree, focal statistics is probably the way to go :

http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/focal-statistics.htm

Use the "Minimum" setting to get what you are looking for, except that will include the value of the central cell itself if this happens to be the lowest value. Is that acceptable?