Looping dentro de r.mapcalc en GRASS GIS

Question

Tengo un mapa de clasificación de la cubierta terrestre que he creado a partir de imágenes multiespectrales de 2 m de resolución y estoy intentando reducir los resultados a un mapa de clases de 20 m de resolución.

El clima en cuestión es árido, por lo que hay muchas zonas con escasos arbustos o árboles. Por ello, una clase de resolución de 2 m puede tener un par de píxeles para un solo arbusto. No me interesa saber dónde está cada uno de los árboles o arbustos, sino dónde hay zonas de alta concentración de una clase determinada (es decir, quiero crear un mapa de clases con una resolución de 20 m que agregue las clases de alta resolución en clases como "arbustos dispersos" o "arbustos densos" en función de la densidad de clases de arbustos de 2 m x 2 m que se encuentren en la cuadrícula de 20 m x 20 m).

Usando r.mapcalc en GRASS GIS Sé que puedes referirte a celdas vecinas usando el formato map[1,-2] y hay muchas funciones útiles disponibles en r.mapcalc. Sin embargo, mi problema es que cuando estoy reduciendo la resolución a ~20m desde una resolución de 2m hay ~100 vecinos para analizar y tendría que dirigirme a cada uno específicamente ya que no hay manera que he encontrado para anidar un bucle for dentro de una llamada a r.mapcalc.

¿Alguien tiene alguna sugerencia sobre una forma de recopilar estadísticas de las celdas que rodean una celda del mapa dada, y cambiar la celda en cuestión en función de sus vecinos de la forma que describo?

Answer 1

Como sugirió @xunilk, encontré en GDAL la herramienta que buscaba. También utilicé Counter para crear una forma de histograma que me ayudó a completar mi análisis.

for i in xrange(0,iRange,resampleInterval):
    for j in xrange(0,jRange,resampleInterval):
        scanline = band.ReadRaster(i, j, windowSize, windowSize, windowSize, windowSize, band.DataType)
        neighbors = struct.unpack(fmtTypes[BandType] * windowSizeSquared, scanline)
        neighborData = Counter(neighbors)
        statsTuple = neighborData.most_common() # Returns all unique items and their counts [(4,101),(1,50)]
        numNeighbors = sum(neighborData.values()) # sum(neighborData.values()) yields number of neighbors
        # findClass assigns classID based on percent of total area each class makes up
        classID = findClass(statsTuple, numNeighbors, classID, 'Roads')
        classID = findClass(statsTuple, numNeighbors, classID, 'Water')

        classID = findClass(statsTuple, numNeighbors, classID, 'Forest')
        classID = findClass(statsTuple, numNeighbors, classID, 'Shrub')
        # If mixture does not fit into one of the above buckets, set class equal to the mode of the window
        if not classID:
            classID = statsTuple[0][0] # Set classID equal to mode of neighbors
        classArray.append((i+1,j+1,classID)) #save x,y coordinate with classID