Mencari alternatif statistik fokus ArcGIS dalam Python open source

Question

Saya berjuang untuk melakukan agregasi piksel raster di python open-source seperti yang mirip dengan fungsi statistik ArcGIS Focal, saya ingin membuat jendela persegi panjang 5 x 5 di mana fungsi program akan menghitung rata-rata piksel pusat menggunakan piksel tetangga jatuh di dalam jendela yang ditentukan. Nilai raster masukan saya dalam format float 0 - 1. Adakah yang bisa menyarankan, cara yang mungkin untuk melakukannya dengan python?

Saya mencoba kode di bawah ini, tidak berfungsi

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

Kode yang dimodifikasi berdasarkan saran @Neprin, namun, saya ingin mengubahnya berdasarkan struktur file saya, Tolong bantu dalam hal ini.

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

Alternatif lain adalah menggunakan opencv2 Image Filtering ( tautan ):

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

The r.neighborsalat dalam GRASS mirip dengan statistik fokus di ArcGIS. Masing-masing memungkinkan statistik dihitung dalam jendela bergerak.

r.neighbours - Menjadikan setiap nilai kategori sel sebagai fungsi dari nilai kategori yang ditetapkan ke sel di sekitarnya, dan menyimpan nilai sel baru dalam lapisan peta raster keluaran.

Answer 3

Anda bisa menggunakan scipy's ndimage.convolve :

from scipy.ndimage import convolve

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

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

Answer 4

Solusi lain adalah dengan menggunakan Orfeo Toolbox. Ini memiliki fungsi yang disebut BandMathX yang melakukan statistik fokus berdasarkan lingkungan mana pun dan berbagai macam fungsi, atau aplikasi penghalusan dengan lebih sedikit pilihan tetapi dengan fungsi rata-rata. Smoothing lebih mudah digunakan tetapi BandMathX lebih fleksibel. Berikut adalah contoh penggunaan aplikasi Smoothing, dengan detail lebih lanjut di sini tentang cara menginstal 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")
app.SetParameterString("type.mean.radius",2)

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