Cálculo de estadísticas focales para vecindad especial?

18

Estoy buscando calcular estadísticas focales para cada celda de un ráster, dentro de un vecindario de un criterio específico.

Antecedentes: tengo tres rásteres 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 en mi área de estudio (suma / total de celdas en el vecindario). El problema es que no puedo usar un círculo simple o un vecindario 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 insumos para un modelo de hábitat, y las áreas fuera de mi área de estudio no pueden considerarse hábitat posible, están urbanizadas. Incluirlos me daría estadísticas erróneas. Entonces, lo que yo 'n determinado por el número de celdas necesarias para cubrir un área igual al tamaño de mi vecindario deseado) que cumplan con mis criterios. El criterio es que no caen dentro de un área urbanizada. Estoy pensando que debería usarse alguna forma de autómata celular. Sin embargo, nunca he trabajado con CA.

Supongo que lo que me gustaría es algo como el código de inicio, o un punto en la dirección correcta.


RESPUESTA AL COMENTARIO A CONTINUACIÓN:

Digamos que estoy calculando esta estadística para una celda en el límite de mi sitio de estudio. Si asigno a cero todas las áreas fuera de mi área de estudio (o ignoro NoData), obtendré una estadística que representa aproximadamente la mitad de la cobertura regional que me interesa. Entonces, porcentaje de cobertura en un área de ~ 10 km ^ 2 , en lugar de 20 km ^ 2 de área. Como estoy estudiando los tamaños de rango de casa, esto es importante. El vecindario tiene que cambiar de forma, ya que así es como el animal ve / usa el paisaje. Si necesitan 20 km ^ 2, cambiarán la forma o su territorio de origen. Si no selecciono ignorar NoData, la salida de la celda será NoData, y NoData no es de ayuda.


"PROGRESO" AL 24/10/2014

Aquí está el código que he creado 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_cov

Desafortunadamente, es INCREÍBLEMENTE lento.

CSB
fuente
1
Ese es un problema interesante. Puede establecer todas las celdas fuera de su área de estudio en NoData, pero no sé cómo va a lograr que un vecindario se adapte y mantenga el mismo tamaño de 20 km2 (tendría que cambiar de forma).
jbchurchill
@CSB jbchurchill tiene razón, lo mejor que puede hacer aquí es asignar valores NoData fuera de su área de estudio. La herramienta de Estadísticas Focales puede tratar esos valores de noda adecuadamente. Consulte 'Procesar celdas de NoData' aquí resources.arcgis.com/en/help/main/10.1/index.html#//…
WhiteboxDev
@WhiteboxDev: su sugerencia no resolverá mi problema. Editaré lo anterior y explicaré por qué eso no funcionará.
CSB
¿Has visto esta publicación, que trata sobre el uso de estadísticas focales con un radio variable ( gis.stackexchange.com/questions/34306/… )? Este parece ser su problema: las celdas en el borde deben tener un radio grande y considerar solo un vecindario semicircular. Por supuesto, dependiendo del tamaño de su celda, es posible que deba crear muchos, muchos rásteres para elegir.
floema
1
@CSB Te encontrarás con efectos de borde independientemente de si usas NoData y un vecindario reducido o si cambias la forma / ubicación de tu vecindario para asegurar el tamaño. Al menos con el primero, no estará sobremuestreando / representando datos de borde cercano de una manera no transparente. Esto es parte del infame problema de la unidad areal modificable.
WhiteboxDev

Respuestas:

0

El código anterior es la respuesta eventual e imperfecta que se me ocurrió para este problema. Al final, pensé que el mejor enfoque era usar un vecindario circular y calcular el área que cruza mi área "disponible". (De todos modos, un vecindario circular daría las n ~ celdas más cercanas, por lo tanto, no es necesario ser demasiado elegante con Autómatas Celulares). Si el área era demasiado pequeña, simplemente crecí el círculo hasta que no lo fue.

Funcionó bien pero, como señalé, fue muy lento. Consulte este hilo para obtener sugerencias sobre cómo acelerarlo. Maximizando el rendimiento del código para Shapely . Seguí las sugerencias, que conducen a este hilo Comprensión del uso de índices espaciales . No terminé aplicando un árbol r al final, porque en realidad nunca terminé usando el código. Descubrí que podía abordar el problema desde un ángulo completamente diferente y ahorrarme mucho tiempo / energía, así que lo hice. Quizás lo termine algún día ...

CSB
fuente
Al leer su código, parece que hay una buena posibilidad de que el uso de un índice espacial pueda acelerar el código, a menudo dramáticamente. Hacer intersecciones contra un MultiPolygon como ese es muy lento.
Snorfalorpagus
@Snorfalorpagus Sí, si miras la respuesta, señalo otros dos hilos relacionados con esta pregunta. Ambos discuten el uso de un índice espacial.
CSB