Calcul des statistiques focales pour un quartier spécial?

Je cherche à calculer des statistiques focales pour chaque cellule d'un raster, dans un voisinage d'un critère spécifié.

Contexte - J'ai trois rasters binaires, chacun représentant un seul type de végétation d'intérêt. J'aimerais calculer le pourcentage de couverture de chaque type de végétation dans (par exemple) 20 km ^ 2 de n'importe quelle cellule de ma zone d'étude (somme / total des cellules dans le voisinage). Le problème est que je ne peux pas utiliser un simple cercle ou voisinage carré autour de chaque cellule parce que, si je le faisais, la zone de recherche utilisée pour calculer la somme incorporerait des zones en dehors de ma zone d'étude. Cette exception est importante car les statistiques seront utilisées comme intrants pour un modèle d'habitat, et les zones à l'extérieur de ma zone d'étude ne peuvent pas être considérées comme un habitat possible - elles sont urbanisées. Les inclure me donnerait des statistiques erronées. Alors, ce que je 'n déterminé par le nombre de cellules nécessaires pour couvrir une zone égale à la taille de quartier souhaitée) qui répondent à mes critères. Les critères étant qu'ils ne relèvent pas d'une zone urbanisée. Je pense qu'une certaine forme d'automates cellulaires devrait être utilisée. Je n'ai cependant jamais travaillé avec CA.

Je suppose que ce que j'aimerais, c'est quelque chose comme du code de démarrage, ou un point dans la bonne direction.

RÉPONSE AU COMMENTAIRE CI-DESSOUS:

Disons que je calcule cette statistique pour une cellule à la limite de mon site d'étude. Si j'attribue à zéro toutes les zones en dehors de ma zone d'étude (ou si j'ignore NoData), j'obtiendrai une statistique qui représente environ la moitié de la couverture surfacique qui m'intéresse. Donc, pourcentage de couverture dans une zone de ~ 10 km ^ 2 , au lieu de 20 km ^ 2. Puisque j'étudie les tailles de domaine domestique, c'est important. Le quartier doit changer de forme, car c'est ainsi que l'animal voit / utilise le paysage. S'ils ont besoin de 20 km ^ 2, ils changeront la forme ou leur territoire d'origine. Si je ne coche pas ignorer NoData, la sortie de la cellule sera NoData - et NoData n'est d'aucune aide.

"PROGRÈS" AU 24/10/2014

Voici le code que j'ai trouvé jusqu'à présent en utilisant Shapely et 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_cov

Malheureusement, c'est incroyablement lent.

Le code ci-dessus est la réponse éventuelle et imparfaite que j'ai trouvée pour ce problème. En fin de compte, j'ai pensé que la meilleure approche était d'utiliser un voisinage circulaire et de calculer la zone intersectant ma zone "disponible". (Un voisinage circulaire donnerait de toute façon les n ~ cellules les plus proches - donc, pas besoin de devenir trop sophistiqué avec les automates cellulaires.) Si la zone était trop petite, j'ai simplement agrandi le cercle jusqu'à ce qu'il ne le soit pas.

Cela a bien fonctionné mais, comme je l'ai noté, c'était très lent. Voir ce fil pour des suggestions sur la façon d'accélérer. Optimisation des performances du code pour Shapely . J'ai suivi les suggestions qui ont mené à ce fil Comprendre l'utilisation des index spatiaux . Je n'ai pas fini par appliquer un r-tree à la fin, car je n'ai finalement jamais utilisé le code. J'ai découvert que je pouvais aborder le problème sous un angle entièrement différent et me faire gagner beaucoup de temps / d'énergie, alors je l'ai fait. Je finirai peut-être un jour ...