¿Cómo amortiguar los píxeles ráster por sus valores?

28

Los píxeles a la izquierda representan ubicaciones de árboles y sus radios de copa asociados (es decir, valores de píxeles que van de 2 a 5). Me gustaría amortiguar estos píxeles de trama por su valor de radio de corona. La imagen a la derecha es lo que espero lograr usando solo métodos de procesamiento de ráster .

Inicialmente, pensaría usar una suma focal circular en ArcGIS, aunque la configuración de vecindario es un valor fijo, que no tendría en cuenta el radio de corona de tamaño variable.

¿Cuál es un buen método para "almacenar" los píxeles por sus valores?

ingrese la descripción de la imagen aquí

Aaron
fuente
2
¿Intentó convertir el ráster en puntos, luego almacenar en búfer por campo y luego convertir nuevamente en ráster?
2
Ayuda a darse cuenta de que esta es una operación no local , porque esa no localidad muestra que existen límites inherentes a cómo se puede llevar a cabo el cálculo. Por ejemplo, su salida cambiaría radicalmente en casi todas partes si solo un píxel aislado en la entrada cambiara a un valor grande. Por lo tanto, si conoce restricciones en los valores de entrada, compártalos, porque eso puede conducir a soluciones mejoradas. Por ejemplo, ¿estarán todos sus valores de entrada siempre en el conjunto {2,3,4}?
whuber
@Dan Patterson Así es como se me ocurrió la imagen a la derecha. Sin embargo, estoy tratando de evitar las operaciones de vectores por completo y evitar esos pasos.
Aaron
2
@whuber Este conjunto de datos representa árboles con diámetros de copa variables. Dado eso, las mediciones del radio de la copa del árbol pueden variar de 1 a 10. También debo agregar que la salida almacenada solo necesita ser 0 para ausencia de corona y 1 para presencia de corona.
Aaron
1
Vale gracias. Parece que creó la salida de ejemplo al unir los 3 búferes de los puntos con el valor 3, los 4 búferes de los puntos con el valor 4 y los 5 búferes de los puntos con el valor 5. (Parece haber olvidado para procesar los puntos con valor 2.) Ese proceso no solo responde a su pregunta, sino que (creo) es la solución más simple usando las herramientas disponibles en Spatial Analyst.
whuber

Respuestas:

14

Aquí hay una solución raster pura para Python 2.7usar numpyy scipy:

import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt

#create tree location matrix with values indicating crown radius
A = np.zeros((120,320))
A[60,40] = 1
A[60,80] = 2
A[60,120] = 3
A[60,160] = 4
A[60,200] = 5
A[60,240] = 6
A[60,280] = 7

#plot tree locations
fig = plt.figure()
plt.imshow(A, interpolation='none')
plt.colorbar()

#find unique values
unique_vals = np.unique(A)
unique_vals = unique_vals[unique_vals > 0]

# create circular kernel
def createKernel(radius):
    kernel = np.zeros((2*radius+1, 2*radius+1))
    y,x = np.ogrid[-radius:radius+1, -radius:radius+1]
    mask = x**2 + y**2 <= radius**2
    kernel[mask] = 1
    return kernel

#apply binary dilation sequentially to each unique crown radius value 
C = np.zeros(A.shape).astype(bool)   
for k, radius in enumerate(unique_vals):  
    B = ndimage.morphology.binary_dilation(A == unique_vals[k], structure=createKernel(radius))
    C = C | B #combine masks

#plot resulting mask   
fig = plt.figure()
plt.imshow(C, interpolation='none')
plt.show()

Entrada: ingrese la descripción de la imagen aquí

Salida: ingrese la descripción de la imagen aquí

jatobat
fuente
1
¡+1 para el enfoque de dilatación! Funciona con puntos cercanos también.
Antonio Falciano
Este es un gran ejemplo de por qué ese viejo esquema de color de jet era terrible. Esto se ve mucho más claro con viridis.
nada101
8

Enfoque basado en vectores

Esta tarea se puede realizar en tres pasos:

Nota: el uso del campo de búfer evita el cálculo de un búfer para cada valor de radio de corona.


Enfoque basado en ráster

Al evitar la solución basada en vectores, este problema sugiere utilizar un tipo de autómatas celulares basados ​​en los vecinos más cercanos. Suponiendo que todos los píxeles negros son ceros, los píxeles son cuadrados y su tamaño es igual a 1 (o, alternativamente, están escalados oportunamente), las reglas a adoptar son muy simples:

  1. Si el valor de píxel ( VALUE) es mayor que 1, su valor se convierte VALUE-1y luego considera sus píxeles circundantes. Si sus valores son menores que VALUE-1, estos píxeles nacen o crecen y su valor se vuelve VALUE-1. De lo contrario, estos píxeles sobreviven y no se modifican.
  2. Si VALUE<=1no hace nada (¡el píxel está muerto!).

Estas reglas deben aplicarse hasta que todos los píxeles estén muertos, es decir, sus valores son iguales a 0 o 1. Entonces N-1, ¿dónde Nestá el valor máximo que tiene en el ráster de entrada? Este enfoque se puede implementar fácilmente con un poco de Python y numpy.

Antonio Falciano
fuente
1
Gracias por la respuesta afalciano. Este método es cómo creé la imagen a la derecha y utilizo un enfoque vectorial, uno que estoy tratando de evitar.
Aaron
1
Ok Aaron, aquí hay un enfoque basado en ráster ahora. Espero que esto ayude.
Antonio Falciano
7

Otra opción sería crear rásteres separados para cada valor de píxel, en este caso 4 rásteres, con una condición. Luego expanda los rásteres por un recuento de píxeles correspondiente al valor del ráster (posiblemente iterando sobre una lista de valores). Por último, une los rásteres (algebraicos o espaciales) para crear un ráster binario para las coronas de los árboles.

HDunn
fuente
1
Esta idea es la correcta. Los detalles se pueden mejorar: (1) la selección crea un indicador binario (0,1) de los árboles de un radio de copa determinado. (2) Una suma focal de esa selección, usando una vecindad circular del radio dado, es rápida de calcular usando la FFT. (3) Agregar las sumas focales (puntiagudas) y compararlas con 0 da el búfer deseado.
whuber
7

Es una pregunta difícil hacer esto en la trama porque no tiene la oportunidad de usar el valor del píxel para definir el tamaño del búfer. Por lo tanto, necesitaría hacer el filtro focal para cada valor, como ya dijo.

Aquí hay una posible respuesta para hacerlo con solo 3 filtros (no pude encontrar menos), pero no perfectamente como lo mencionó Whuber: sus amortiguadores se truncarán cuando los árboles estén cerca unos de otros.

1) EDITAR: asignación euclidiana (esto no resuelve completamente el problema, ya que corta los amortiguadores en la vecindad de árboles más pequeños, pero es mejor que los artefactos de mi primera solución).

2) distancia euclidiana alrededor de cada píxel

3) calculadora ráster (álgebra de mapas) con una declaración condicional

Con("allocation_result" > ("distance_result" / pixel_size) , 1 , 0)

Tenga en cuenta que puede ajustar la condición según sus necesidades en términos de radio (con o sin el píxel central)

radouxju
fuente
+1 Este es un enfoque creativo. Probaré para ver si es posible ampliar usando este enfoque.
Aaron
2
El enfoque de distancia euclidiana no funcionará, porque calcula solo la distancia al árbol más cercano , que no es necesariamente la distancia a un árbol cuya corona cubre el punto.
whuber
2

¿Me pregunto por qué no usas la herramienta de expansión de ArcGIS ?

import arcpy
from arcpy.sa import *

raster_in  = r'c:\test.tif'
raster_out = r'c:\test_out.tif'

outExpand1 = Expand(raster_in, 2, 2)
outExpand2 = Expand(outExpand1, 3, 3)
outExpand3 = Expand(outExpand3, 4, 4)
outExpand4 = Expand(outExpand4, 5, 5)

outExpand4.save(raster_out)

En caso de superposición: el último expandcomando cubrirá los anteriores.

Señor che
fuente
2

Si tiene la posición del píxel, el radio y el algoritmo del círculo del punto medio (una variante del Bresenham Alg.) Le dan una pista. En mi opinión, es fácil crear un polígono a partir de este enfoque y creo que es fácil implementarlo en Python. Una unión de este conjunto de polígonos le proporciona el área de cobertura.

huckfinn
fuente
Sé que no es cuestión de pregunta, pero, ¿quieres saber más sobre primitivas gráficas y relleno de polígonos de líneas de escaneo? Para cirles es muy fácil. La combinación convexa es una palabra de moda y así sucesivamente ...
huckfinn
¿Cómo se aplicaría esto usando operaciones raster básicas?
whuber
Si intenta manejar esto en el espacio ráster, determine los puntos del círculo, clasificándolos por la y o xy llene el espacio con una línea recta (línea de exploración) para completar el círculo. En el enfoque triangular, si construye el círculo mediante una aproximación de sectores tringulares e intenta llenar el triángulo, necesita una prueba si el punto está dentro o fuera (combinación convexa) y es al revés. Y en el enfoque "SIG", construir polígonos (polígonos orientados en sentido horario) y hacer una unión es el tercero (IMO, el más costoso).
huckfinn
Para ser claros: Y en el enfoque de "SIG" ... hacer una operación algebrarica como unión, intersección, contacto ... es la tercera OMI la más costosa.
huckfinn