¿Contando valores de píxeles consecutivos para un conjunto de rásteres con ArcGIS Spatial Analyst?

23

Estoy usando ArcGIS 10 con Spatial Analyst.

Tengo un conjunto de rásteres (8 en total) que solo contienen un 1 o un 0 para cada celda. Cada ráster representa un valor de datos de años diferentes. Por razones de sake del año 1 al año 8.

Puedo agregar todos los rásteres juntos, lo que me dará una grilla final con valores que van de 0 a 8. Un 8 que indica que la celda era constantemente 1 para el conjunto de rásteres (todos los años).

Me gustaría encontrar para cada celda el número consecutivo más largo de 1.

Entonces, por ejemplo, la cuadrícula total puede registrar para una sola celda un valor de digamos 5 pero en las 8 cuadrículas esa celda tiene el mayor número consecutivo de 1 igual a 3. O otra forma de expresar esto es durante 3 años que la celda fue un 1 entonces comenzó a oscilar entre ceros y unos.

Mis habilidades de procesamiento de ráster no son tan buenas como mis habilidades de procesamiento de vectores y he echado un buen vistazo al archivo de ayuda de ESRI, pero no puedo entender cómo se podría lograr esto usando las herramientas de procesamiento geográfico disponibles.

¿Algunas ideas?

Hornbydd
fuente
1
Eso es realmente un análisis bastante bueno. Como de costumbre, hay más de una forma de hacer lo que estás tratando de hacer. Sin embargo, creo que necesitarás hacer algo de programación para recorrer todas las combinaciones.
MLowry
1
Un comentario general (inspirado en este comentario de @MLowry): por favor, vote las preguntas cada vez que parezcan interesantes o estén claramente articuladas. Las buenas preguntas impulsan todo en nuestro sitio; ¡Por favor hagamos lo que podamos para recompensar a los que les piden!
whuber

Respuestas:

14

Como se trata de una operación local, veamos cómo hacerlo para una sola celda: Map Algebra se encargará del resto.

Primero tenga en cuenta que el orden de los rásteres obviamente es importante. Por lo tanto, las estadísticas de celdas de un solo disparo, como una suma de celdas, no lo harán.

Si tuviéramos que encontrar una secuencia como 01110101 en una celda determinada, la procesaríamos de principio a fin y

  1. Comience con un recuento en cero.

  2. Incremente el conteo cada vez que encontremos un 1.

  3. Restablezca el recuento cada vez que encontremos un 0, después de guardar el último recuento .

  4. Al final, tome el recuento máximo guardado (incluido el recuento final).

El paso 1 se implementa con una cuadrícula cero constante. Los pasos 2 y 3 dependen de lo que encontremos: por lo tanto, es una operación condicional . El paso 4 claramente es un máximo local. Codificaremos esto, entonces, un poco más formalmente como:

count = 0
result = 0
For each value:
    If (value==1):
        count=count+1
    else
        result = max(result, count)
        count=0
result = max(result, count)

Es mejor hacerlo con un script de Python cuando tiene muchas cuadrículas, pero con ocho no es oneroso desenrollar el bucle y escribir los pasos a mano. Esto revela un pequeño problema: result=max(longest,count)es una especie de "efecto secundario" que es difícil de codificar con operaciones ráster. (Pero se puede hacer, como se muestra en la segunda solución a continuación). También es ineficiente, porque agrega un cálculo adicional en cada paso. Por lo tanto, modificamos un poco el enfoque, con el objetivo de posponer la maxoperación hasta el final. Esto requerirá guardar un recuento separado en cada etapa.

Al pasar por este proceso, también encontré un acceso directo para el primer paso. Esto lleva a la siguiente solución, que aunque es un poco larga e intensiva en RAM, es simple e implica pasos ejecutados rápidamente:

result1 = "grid1"
result2 = con("grid2"==1, "result1"+1, 0)
result3 = con("grid3"==1, "result2"+1, 0)
result4 = con("grid4"==1, "result3"+1, 0)
result5 = con("grid5"==1, "result4"+1, 0)
result6 = con("grid6"==1, "result5"+1, 0)
result7 = con("grid7"==1, "result6"+1, 0)
result8 = con("grid8"==1, "result7"+1, 0)
CellStatistics(["result1", "result2", "result3", "result4", "result5", "result6", "result7" "result8"], "max")

La sintaxis real varía con su versión de ArcMap. (Por ejemplo, CellStatisticscreo que es nuevo en la versión 10, pero siempre ha estado disponible una operación máxima local).

En el ejemplo con la entrada 01110101, la secuencia de cuadrículas "resultado *" contendría los valores 0, 1, 2, 3, 0, 1, 0, 1, por lo que al final CellStatisticsdevolvería 3, la longitud de la cadena más larga de 1's.

Si la RAM es escasa, la solución se puede modificar para reutilizar los resultados intermedios, a un costo de aproximadamente duplicar el tiempo de ejecución:

result = "grid1"
temp = con("grid2"==1, "result"+1, 0)
result = CellStatistics[["temp", "result"], "max"]
temp = con("grid3"==1, "temp"+1, 0)
result = CellStatistics[["temp", "result"], "max"]
...
temp = con("grid8"==1, "temp"+1, 0)
CellStatistics[["temp", "result"], "max"]

En el ejemplo con la entrada 01110101, los valores ("temp", "result") serían (NoData, 0) después de la primera línea y después de cada par de operaciones ("Con", "CellStatistics") los valores serían (1 , 1), (2, 2), (3, 3), (0, 3), (1, 3), (0, 3), (1, 3). Una vez más, el valor final es 3.

El patrón regular de las expresiones de Álgebra de mapas en cualquiera de las soluciones indica cómo codificar el algoritmo como un bucle en un script, cambiando los índices según corresponda con cada iteración.

whuber
fuente
Puede haber un error tipográfico en el bloque de código de tipo: count = count = 1 probablemente debería ser count = count + 1
MLowry
1
@ML Gracias (¡buenos ojos!): Ya está arreglado. Es difícil hacer que el pseudocódigo sea absolutamente correcto; La revisión humana es un activo real en la búsqueda de errores. Además, aunque no probé las soluciones en ArcGIS, implementé la primera solución en R, por lo que tengo la seguridad de que el enfoque es correcto.
whuber
1
"whuber" una vez más, eres el hombre que sabe! ¡Que Dios nos ayude al resto de nosotros si alguna vez te atropella un autobús! Su enfoque inicial de Python fue la dirección en la que estaba pensando, pero sé que con los rásteres a menudo se puede hacer todo mucho más ingenioso que usted ha demostrado. ¡Si se encuentra en Gran Bretaña, sería un honor comprarle una pinta de nuestra mejor cerveza plana a temperatura ambiente! :)
Hornbydd
Gracias, Duncan: ¡pero mira la excelente solución de Andy Harfoot!
whuber
14

Solo charlando sobre esto y preguntándome si podría abordar el problema tratando las cuadrículas de entrada como una secuencia binaria. Esto le permitiría combinarlos para obtener un número entero de resumen único para la secuencia, es decir, 01110101 = 117. Este valor podría reclasificarse para obtener el número máximo de 1 consecutivos.

Aquí hay un ejemplo que muestra una forma de combinar ocho cuadrículas:

2*(2*(2*(2*(2*(2*(2*"g8" + "g7") + "g6") + "g5") + "g4") + "g3") + "g2") + "g1"

Las operaciones bit a bit también podrían ser puestas en servicio para este paso. Alternativamente, puede usar combinar seguido de un cálculo de campo. (El cálculo del campo tendrá una expresión similar a la anterior).

La tabla de reclasificación debe proporcionar longitudes de ejecución máximas para todos los valores entre 00000000B = 0 y 11111111B = 255. En orden, aquí están:

0, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 6, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 4, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 7, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 6, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 7, 8

Este enfoque está limitado a aproximadamente 20 cuadrículas en ArcGIS: usar más de esto puede crear una tabla de atributos difícil de manejar. ( Combinese limita específicamente a 20 cuadrículas).

Andy Harfoot
fuente
Un enorme +1: esta es una muy buena idea. (La única limitación es que cuando hay más de 31 cuadrículas involucradas, se quedará sin bits para usar). Me he tomado la libertad de desarrollar un poco su idea para que otros puedan ver lo fácil que es implementarla.
whuber
3

¿Has pensado en cambiar los valores de 0 y 1 a valores con la potencia de 2 (1,2,4,8,16,32). Cuando combine las 8 cuadrículas, obtendrá valores únicos para cada celda que le proporcionarán información consecutiva (es decir: un valor de 3 significa año 1 y 2, donde un valor de 54 sería de 6 a 8 años).

Solo un pensamiento

Ryan Garnett
fuente
Esto es precisamente lo que sugirió @Andy Harfoot unas horas antes, Ryan. :-)
whuber
Gracias y perdon Estoy leyendo esto en mi teléfono de vacaciones.
Ryan Garnett