Quiero calcular las estadísticas focales de cada celda de un raster, dentro de una vecindad de un criterio especificado.
Antecedentes - Tengo tres rásters binarios, cada uno de los cuales representa un único tipo de vegetación de interés. Me gustaría calcular el porcentaje de cobertura de cada tipo de vegetación dentro de (por ejemplo) 20 km^2 de cualquier celda de mi área de estudio (suma/total de celdas en el vecindario). El problema es que no puedo utilizar un simple círculo o cuadrado alrededor de cada celda porque, si lo hiciera, el área de búsqueda utilizada para calcular la suma incorporaría áreas fuera de mi área de estudio. Esta excepción es importante porque las estadísticas se utilizarán como datos de entrada para un modelo de hábitat, y las zonas situadas fuera de mi zona de estudio no pueden considerarse posibles hábitats: están urbanizadas. Incluirlas me daría estadísticas erróneas. Por lo tanto, lo que quiero hacer es escribir un código en python que elija un barrio que represente la zona de estudio. n celdas más cercanas ( n determinado por el número de celdas necesarias para cubrir un área igual al tamaño de mi barrio deseado) que cumplen mis criterios. El criterio es que no estén dentro de una zona urbanizada. Creo que habría que utilizar algún tipo de autómata celular. Sin embargo, nunca he trabajado con AC.
Supongo que lo que me gustaría es algo así como un código de inicio, o un punto en la dirección correcta.
RESPONDER AL COMENTARIO A CONTINUACIÓN:
Supongamos que calculo esta estadística para una celda situada en el límite de mi zona de estudio. Si asigno el valor cero a todas las zonas situadas fuera de mi zona de estudio (o ignoro NoData), obtendré una estadística que representa aproximadamente la mitad de la cobertura de área que me interesa. Es decir, el porcentaje de cobertura en un área de ~10 km^2, en lugar de un área de 20 km^2. Como estoy estudiando el tamaño del área de distribución, esto es importante. El vecindario tiene que cambiar de forma, ya que así es como el animal ve/utiliza el paisaje. Si necesitan 20 km^2, van a cambiar la forma o su territorio de origen. Si no marco ignorar NoData, la salida de la celda será NoData - y NoData no es de ayuda.
"PROGRESO" A 24/10/2014
Aquí está el código que he encontrado hasta ahora usando Shapely y Fiona:
import numpy as np
import pprint
import shapely
from shapely.geometry import*
import fiona
from fiona import collection
import math
traps = fiona.open('C:/Users/Curtis/Documents/ArcGIS/GIS_Data/occurrence/ss_occ.shp', 'r')
study_area = fiona.open('C:/Users/Curtis/Documents/ArcGIS/GIS_Data/Study_Area.shp', 'r')
for i in study_area: #for every record in 'study_area'
sa = shape(i['geometry']) #make a variable called 'sa' that is a polygon
grassland = fiona.open('C:/Users/Curtis/Documents/ArcGIS/GIS_Data/land_cover/polys_for_aa/class3_aa.shp', 'r')
pol = grassland.next()
gl = MultiPolygon([shape(pol['geometry']) for pol in grassland])
areaKM2 = 20
with traps as input:
r = (math.sqrt(areaKM2/math.pi))*1000
for point in input:
pt = shape(point['geometry'])
pt_buff = pt.buffer(r)
avail_area = pt_buff.intersection(sa).area
# works to here
while avail_area < areaKM2:
r += 10
pt_buff = pt.buffer(r)
avail_area = pt_buff.intersection(sa).area
perc_cov = pt_buff.intersection(gl).area//areaKM2
print perc_covPor desgracia, es INCREÍBLEMENTE lento.
