Cerco alternativa alle statistiche focali ArcGIS in Python open source

Question

Sto lottando per eseguire l'aggregazione dei pixel del raster in python open source come fa la funzione di statistica Focal di ArcGIS, vorrei creare una finestra rettangolare 5 x 5 su cui la funzione del programma calcolerà la media del pixel centrale utilizzando pixel vicini che cadono all'interno della finestra definita. I miei valori raster di input sono in formato float 0 - 1. Qualcuno può suggerire un modo possibile per farlo in Python?

Ho provato il codice seguente, non funziona

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

Codice modificato in base al suggerimento di @Neprin, tuttavia, vorrei modificarlo in base alla struttura del file, per favore aiutatemi

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

un'altra alternativa è usare 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

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

Lo r.neighborsstrumento in GRASS è simile alle statistiche focali in ArcGIS. Ciascuno consente le statistiche calcolate all'interno di una finestra mobile.

r.neighbors - Rende il valore di ciascuna categoria di celle una funzione dei valori di categoria assegnati alle celle circostanti e memorizza i nuovi valori di cella in un livello di mappa raster di output.

Answer 3

Potresti usare ndimage.convolve di scipy :

from scipy.ndimage import convolve

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

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

Answer 4

Un'altra soluzione è utilizzare Orfeo Toolbox. Ha una funzione chiamata BandMathX che esegue la statistica focale in base a qualsiasi vicinanza e una grande varietà di funzioni, o l'applicazione di levigatura con meno scelte ma con la funzione media. Il livellamento è più facile da usare ma BandMathX è più flessibile. Ecco un esempio sull'uso dell'applicazione Smoothing, con maggiori dettagli qui su come installare l'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()