¿Cálculo de estadísticas focales para vecindario especial?

18

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

Fondo: tengo tres rásteres binarios, cada uno representa un único tipo de interés de vegetación. 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 datos de entrada para un modelo de hábitat, y las áreas fuera de mi área de estudio no se pueden considerar hábitat posible, están urbanizadas. Incluirlos me daría estadísticas erróneas. Entonces, lo que estoy buscando hacer es escribir un código en Python que elija un vecindario que represente las n celdas más cercanas ( n determinado por el número de celdas necesarias para cubrir un área igual al tamaño de mi vecindario deseado) que cumpla con mis criterios. El criterio es que no caen dentro de un área urbanizada. Estoy pensando que debería usarse algún tipo de autómata celular. Aunque 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 todas las áreas fuera de mi área de estudio a cero (o ignoro NoData), obtendré una estadística que representa aproximadamente la mitad de la cobertura de área en la que estoy interesado. Entonces, porcentaje de cobertura en un área de ~ 10 km ^ 2 , en lugar de 20 km ^ 2 área. Ya que estoy estudiando 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 será de ayuda.

"PROGRESO" A PARTIR DEL 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_cov

Desafortunadamente, es INCREÍBLEMENTE lento.

    
pregunta CSB 16.10.2014 - 21:06

1 respuesta

0

El código anterior es la respuesta eventual e imperfecta que surgió para este problema. Al final, pensé que el mejor enfoque era utilizar un vecindario circular y calcular el área que se intersecaba con mi área "disponible". (Un vecindario circular daría las celdas más cercanas de todos modos, así que no hay necesidad de volverse demasiado sofisticado con los autómatas celulares). Si el área era demasiado pequeña, solo crecí el círculo hasta que no fue así.

Funcionó bien pero, como señalé, fue muy lento. Vea este hilo para sugerencias sobre cómo acelerarlo. Maximizar el rendimiento del código para Shapely . Seguí las sugerencias, que llevan a este hilo Comprensión del uso de índices espaciales . No terminé aplicando un r-tree 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, y así lo hice. Tal vez lo termine algún día ...

    
respondido por el CSB 23.09.2015 - 20:14

Lea otras preguntas en las etiquetas