2

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
2

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

data = band.ReadAsArray()

#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[0]), -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).

enter image description here

Resulting raster was obtained as expected (by using Value Tool plugin to corroborate cell raster values):

enter image description here

  • Thank you, this gives me a great overview of how to code this myself. – user57029 Dec 16 '16 at 5:37
1

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?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.