I am struggling to do the pixel aggregating of raster in open-source python as similar to the ArcGIS Focal statistics function does, I would like to make a 5 x 5 rectangular window on which the program's function will calculate the mean of the center pixel using neighbor pixels falling inside the defined window. My input raster values are in float format 0 - 1. Please can anyone suggest, a possible way to do it in python?

I tried the below code, it's not working

import time 
import glob
import os
import gdal
import osr
import numpy as np 

start_time_script = time.clock()


for rasterfile in glob.glob(os.path.join(path_ras,'*.tif')):

print ('Processing:'+ ' ' + str(rasterfile_name))

ds = gdal.Open(rasterfile,gdal.GA_ReadOnly)
ds_xform = ds.GetGeoTransform()

print (ds_xform)

ds_driver = gdal.GetDriverByName('Gtiff')
srs = osr.SpatialReference()

ds_array = ds.ReadAsArray()

sz = ds_array.itemsize

print ('This is the size of the neighbourhood:' + ' ' + str(sz))

h,w = ds_array.shape

print ('This is the size of the Array:' + ' ' + str(h) + ' ' + str(w))

bh, bw = 5,5

shape = (h/bh, w/bw, bh, bw)

print ('This is the new shape of the Array:' + ' ' + str(shape))

strides = sz*np.array([w*bh,bw,w,1])

blocks = np.lib.stride_tricks.as_strided(ds_array,shape=shape,strides=strides)

resized_array = ds_driver.Create(rasterfile_name + '_resized_to_52m.tif',shape[1],shape[0],1,gdal.GDT_Float32)
band = resized_array.GetRasterBand(1)

zero_array = np.zeros([shape[0],shape[1]],dtype=np.float32)

print ('I start calculations using neighbourhood')
start_time_blocks = time.clock()

for i in xrange(len(blocks)):
    for j in xrange(len(blocks[i])):

        zero_array[i][j] = np.mean(blocks[i][j])

print ('I finished calculations and I am going to write the new array')


end_time_blocks = time.clock() - start_time_blocks

print ('Image Processed for:' + ' ' + str(end_time_blocks) + 'seconds' + '\n')

end_time = time.clock() - start_time_script
print ('Program ran for: ' + str(end_time) + 'seconds')  

MOdified code based on @Neprin suggestion, however, I would like to modify it based on my file structure, Please help on this

import numpy as np
import gdal
import cv2
import matplotlib.pyplot as plt
import seaborn as sns

img = gdal.Open('20180305.tif').ReadAsArray() # i have multiple raster i.e.20180305, 20180306, 20180305 so on  
 # i want put give the path of folder where i kept my input raster 
img2 = np.zeros(np.array(img.shape) + 10)
img2[5:-5,5:-5] = img  # fix edge interpolation
kernel = np.ones((5,5),np.float32)
dst = cv2.filter2D(img2,-1,kernel)/25

# Save the output raster in same name as input with projection
  • 1
    Have you looked at r.neighbors in GRASS?: grass.osgeo.org/grass78/manuals/r.neighbors.html
    – Aaron
    Commented Sep 9, 2020 at 15:23
  • Yes, just I looked at the function. It is the same as ArcGIS focal statistics do. can I get a script for that where I can apply it for multiple rasters? Thanks
    – SWAT
    Commented Sep 9, 2020 at 15:26

4 Answers 4


You could use scipy's ndimage.convolve:

from scipy.ndimage import convolve

weights = np.ones((5, 5))

focal_mean = convolve(ds_array, weights) / np.sum(weights)
  • This is a very interesting suggestion because it means that the focal stats can happen in pure python with one additional plugin. It's a shame the detail on convolve is severely lacking. Id like to use this, but I need to have a good idea about what's happening to place some faith in the results.
    – anakaine
    Commented Oct 14, 2021 at 10:36
  • Additionally, could this be modified to work with a radius instead of a square window?
    – anakaine
    Commented Oct 14, 2021 at 11:02
  • 1
    Yes, just modify weights to be a "circle" of ones with zeros elsewhere. Wikipedia may help if the scipy docs aren't sufficient
    – mikewatt
    Commented Oct 14, 2021 at 16:35

another alternative is using opencv2 Image Filtering (link):

import numpy as np
import cv2
import matplotlib.pyplot as plt
import seaborn as sns

img = np.diag([1, 1, 1, 1, 1, 1, 1]).astype('float')
img2 = np.zeros(np.array(img.shape) + 10)
img2[5:-5,5:-5] = img  # fix edge interpolation
kernel = np.ones((5,5),np.float32)
dst = cv2.filter2D(img2,-1,kernel)/25

sns.heatmap(img2[5:-5, 5:-5], annot=True, cbar=False)
sns.heatmap(dst[5:-5, 5:-5], annot=True, cbar=False)


  • I have tried it, i hope it will work, however, I would like to modify it based on my requirement, please help on on this, i have posted my modified code in above
    – SWAT
    Commented Sep 10, 2020 at 7:49
  • Can this be modified to suit a circular pattern (radius)? Also, what happens at the edges without the edge fix? Will the edge fix affect the dimensions of the output?
    – anakaine
    Commented Oct 14, 2021 at 11:01

The r.neighbors tool in GRASS is similar to focal statistics in ArcGIS. Each allows statistics calculated within a moving window.

r.neighbors - Makes each cell category value a function of the category values assigned to the cells around it, and stores new cell values in an output raster map layer.


Another solution is to use Orfeo Toolbox. It has a function called BandMathX which performs the focal stat based on any neighborhood and a large variety of function, or the smoothing application with less choices but with the mean function. Smoothing is easier to use but BandMathX is more flexible. Here is an example on the use of the Smoothing application, with more details here on how to install the Python API.

# The python module providing access to OTB applications is otbApplication
import otbApplication as otb

# Let's create the application with codename "Smoothing"
app = otb.Registry.CreateApplication("Smoothing")

# We set its parameters
app.SetParameterString("in", "my_input_image.tif")
app.SetParameterString("type", "mean")
app.SetParameterString("out", "my_output_image.tif")

# This will execute the application and save the output file

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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