Suche nach einer alternativen ArcGIS-Statistik für Statistiken in Open Source Python

Question

Ich habe Probleme, die Pixelaggregation von Raster in Open-Source-Python ähnlich wie bei der ArcGIS Focal-Statistikfunktion durchzuführen. Ich möchte ein rechteckiges 5 x 5-Fenster erstellen, in dem die Programmfunktion den Mittelwert des mittleren Pixels mithilfe berechnet Nachbarpixel, die in das definierte Fenster fallen. Meine Eingabe-Rasterwerte haben das Float-Format 0 - 1. Kann mir jemand einen möglichen Weg vorschlagen, dies in Python zu tun?

Ich habe den folgenden Code ausprobiert, er funktioniert nicht

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

Modifizierter Code basierend auf dem @ Neprin-Vorschlag. Ich möchte ihn jedoch basierend auf meiner Dateistruktur ändern. Bitte helfen Sie dabei

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

Eine andere Alternative ist die Verwendung von 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

Das r.neighborsTool in GRASS ähnelt der Fokusstatistik in ArcGIS. Jedes erlaubt Statistiken, die in einem sich bewegenden Fenster berechnet werden.

r.neighbors - Macht jeden Zellenkategoriewert zu einer Funktion der Kategoriewerte, die den Zellen um ihn herum zugewiesen sind, und speichert neue Zellenwerte in einem Ausgabe-Raster-Map-Layer.

Answer 3

Sie könnten scipys ndimage.convolve verwenden :

from scipy.ndimage import convolve

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

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

Answer 4

Eine andere Lösung ist die Verwendung der Orfeo Toolbox. Es hat eine Funktion namens BandMathX, die den Fokusstatus basierend auf einer beliebigen Nachbarschaft und einer Vielzahl von Funktionen ausführt, oder die Glättungsanwendung mit weniger Auswahlmöglichkeiten, aber mit der mittleren Funktion. Das Glätten ist einfacher zu verwenden, aber BandMathX ist flexibler. Hier ist ein Beispiel für die Verwendung der Glättungs Anwendung, mit mehr Details hier , wie der Python - API zu installieren.

# 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()