Estadística de píxeles vecinos, en bloques 3x3, mediante ndimage del módulo scipy

Cuando se requiere determinar las estadísticas de píxeles vecinos en un ráster, el efecto a considerar cuando se recorre cada uno de los píxeles es el de borde. En los bordes, por ejemplo, un bloque de 3x3 estará truncado y hay que tomarlo en consideración al determinar las estadísticas.

Con Python y numpy una forma de soslayar esto es añadiendo un "anillo" de valores que se considerán como NoData; tal como se presenta en el código siguiente:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

import numpy as np matrix = [[13, 21, 13, 8], [ 5, 10, 22, 14], [21, 33, 9, 0], [ 0, 0, 0, 0]] matrix = np.insert(matrix, 0, -999, axis=0) matrix = np.insert(matrix, 0, -999, axis=1) matrix = np.insert(matrix, len(matrix[0]), -999, axis=1) matrix = np.insert(matrix, len(matrix), -999, axis=0) print matrix means = [ [] for i in range(len(matrix)-2)] for i in range(len(matrix) - 2): for j in range(len(matrix[0]) - 2): sub_matrix = [ matrix[i, j], matrix[i, j+1], matrix[i, j+2], matrix[i+1, j], matrix[i+1, j+2], matrix[i+2, j], matrix[i+2, j+1], matrix[i+2, j+2] ] sub_matrix = [ element for element in sub_matrix if element != -999 ] means[i].append(np.mean(sub_matrix)) print np.array(means)

La ejecución del código anterior produce como resultado final un arreglo de medias que puede ser convertido, si se desea, en ráster.

1 2 3 4 5 6 7 8 9 10 11

[[-999 -999 -999 -999 -999 -999] [-999 13 21 13 8 -999] [-999 5 10 22 14 -999] [-999 21 33 9 0 -999] [-999 0 0 0 0 -999] [-999 -999 -999 -999 -999 -999]] [[ 12. 12.6 15. 16.33333333] [ 19.6 17.125 13.5 10.4 ] [ 9.6 8.375 9.875 9. ] [ 18. 12.6 8.4 3. ]]

No obstante, varias líneas y todo el loop del código anterior pueden ser condensadas en una sola línea de código si se emplea el método 'generic_filter' del módulo ndimage de scipy. Esto permite determinar un filtro multidimensional para las funciones que acepta como segundo parámetro. Como se quiere excluir los efectos de bordes entonces los valores truncados se consideran como NaN y eso se señala explícitamente en 'cval' como np.NaN. Por tanto, todas las posibles funciones a considerar corresponden a las nanfunciones de numpy.

Usando como filtro la palabra 'nan' y el módulo de expresiones regulares (re) se puede tener una visión panorámica de cuales se pueden usar; tal como se presenta a continuación:

1 2 3 4

['isnan', 'nan', 'nan_to_num', 'nanargmax', 'nanargmin', 'nancumprod', 'nancumsum', 'nanmax', 'nanmean', 'nanmedian', 'nanmin', 'nanpercentile', 'nanprod', 'nanstd', 'nansum', 'nanvar']

En el código siguiente se va a emplear la función 'nanmean' para determinar las medias por bloques de 3x3 píxeles en un ráster de 20 x20 y que excluirá automáticamente los NaN en los bordes sin producir errores en tiempo de ejecución. El ráster de medias se produce como resultado final.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

from osgeo import gdal, osr import numpy as np from scipy import ndimage path = '/home/zeito/pyqgis_data/aleatorio.tif' dataset = gdal.Open(path) transform = dataset.GetGeoTransform() band = dataset.GetRasterBand(1) cols = dataset.RasterXSize rows = dataset.RasterYSize data = band.ReadAsArray() data = data.astype(np.float64) result = ndimage.generic_filter(data, np.nanmean, size = 3, mode='constant', cval=np.NaN) # Create gtif file driver = gdal.GetDriverByName("GTiff") output_file = "/home/zeito/pyqgis_data/mean_raster.tif" dst_ds = driver.Create(output_file, cols, rows, 1, gdal.GDT_Float32) #writting output raster dst_ds.GetRasterBand(1).WriteArray( result ) #setting extension of output raster # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution dst_ds.SetGeoTransform(transform) wkt = dataset.GetProjection() # setting spatial reference of output raster srs = osr.SpatialReference() srs.ImportFromWkt(wkt) dst_ds.SetProjection( srs.ExportToWkt() ) #Close output raster dataset dataset = None dst_ds = None

Después de ejecutado el código anterior en la Python Console de QGIS, el ráster resultante se observa a continuación. Los valores de los píxeles son tal como se esperaban.