Procurando alternativas de estatísticas focais ArcGIS em Python de código aberto

Question

Estou lutando para fazer a agregação de pixels de raster em python de código aberto semelhante à função de estatísticas do ArcGIS Focal, gostaria de fazer uma janela retangular 5 x 5 na qual a função do programa calculará a média do pixel central usando pixels vizinhos caindo dentro da janela definida. Meus valores raster de entrada estão no formato float 0-1. Alguém pode sugerir uma maneira possível de fazer isso em python?

Tentei o código abaixo, não está funcionando

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

start_time_script = time.clock()

path_ras=r'D:\Firm_SM\F1A/'

for rasterfile in glob.glob(os.path.join(path_ras,'*.tif')):
    rasterfile_name=str(rasterfile[rasterfile.find('IMG'):rasterfile.find('.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()
#srs.ImportFromEPSG(4726)

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)
resized_array.SetGeoTransform((ds_xform[0],ds_xform[1]*2,ds_xform[2],ds_xform[3],ds_xform[4],ds_xform[5]*2))
resized_array.SetProjection(srs.ExportToWkt())
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')

band.WriteArray(zero_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')  

Código MOdificado com base na sugestão de @Neprin, no entanto, gostaria de modificá-lo com base na minha estrutura de arquivo, por favor ajude nisto

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
  

Answer 1

outra alternativa é usar a filtragem de imagem opencv2 ( 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

plt.figure(figsize=(15,10))
plt.subplot(121)
sns.heatmap(img2[5:-5, 5:-5], annot=True, cbar=False)
plt.title('original')
plt.subplot(122)
sns.heatmap(dst[5:-5, 5:-5], annot=True, cbar=False)
plt.title('focal')

Answer 2

A r.neighborsferramenta no GRASS é semelhante às estatísticas focais no ArcGIS. Cada um permite estatísticas calculadas dentro de uma janela móvel.

r.neighbors - torna cada valor de categoria de célula uma função dos valores de categoria atribuídos às células ao seu redor e armazena novos valores de célula em uma camada de mapa raster de saída.

Answer 3

Você pode usar o ndimage.convolve de scipy :

from scipy.ndimage import convolve

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

focal_mean = convolve(ds_array, weights) / np.sum(weights)

Answer 4

Outra solução é usar o Orfeo Toolbox. Ele tem uma função chamada BandMathX que executa a estatística focal com base em qualquer vizinhança e uma grande variedade de funções, ou o aplicativo de suavização com menos opções, mas com a função média. A suavização é mais fácil de usar, mas BandMathX é mais flexível. Aqui está um exemplo de uso do aplicativo Smoothing, com mais detalhes aqui sobre como instalar a API Python.

# 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")
app.SetParameterString("type.mean.radius",2)

# This will execute the application and save the output file
app.ExecuteAndWriteOutput()