¿Encuentra ubicaciones de valores más altos en ráster utilizando ArcGIS Desktop?

12

Con ArcGIS 10, tengo un ráster donde me gustaría encontrar el píxel con el valor máximo en el ráster y devolver su ubicación (centro del píxel) en grados decimales. Me gustaría repetir este proceso devolviendo la ubicación del segundo valor más alto del ráster, luego el tercero, y así sucesivamente para que al final tenga una lista de N ubicaciones que tienen los valores más altos en el ráster en orden.

Me imagino que esto podría hacerse más fácilmente usando un script Python, pero estoy abierto a otras ideas si hay una mejor manera.

mga
fuente
¿Has intentado convertir la cuadrícula en puntos y luego agregar campos X, Y y ordenar?
Jakub Sisak GeoGraphics
¿Los valores ráster son flotantes o enteros?
whuber
@Jakub - No, no lo he hecho. Probablemente solo me interesará el 1% más o menos de los puntos, así que no sé si vale la pena agregar los campos x, y para todos los puntos y luego ordenarlos. Tal vez si no hay una opción más eficiente?
Jueves
@whuber: los valores ráster son flotantes.
mga
@mga, vale la pena intentarlo. La conversión es bastante rápida y agregar XY también es una herramienta predeterminada. Eliminar registros no deseados es una operación sencilla y todo podría agruparse en un solo modelo. Solo una idea.
Jakub Sisak GeoGraphics

Respuestas:

5

Si está contento de usar R , hay un paquete llamado raster . Puede leer un ráster con el siguiente comando:

install.packages('raster')
library(raster)
test <- raster('F:/myraster')

Luego, cuando lo veas (escribiendo test), puedes ver la siguiente información:

class       : RasterLayer 
dimensions  : 494, 427, 210938  (nrow, ncol, ncell)
resolution  : 200, 200  (x, y)
extent      : 1022155, 1107555, 1220237, 1319037  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
values      : F:/myraster 
min value   : 0 
max value   : 1 

Puede haber mejores formas de manipular el ráster, pero una forma de encontrar la información que desea podría ser encontrar el valor más alto y obtener su ubicación de matriz, y luego agregarlo a las extensiones más bajas.

djq
fuente
1
+1 para esta referencia. Una vez que haya leído el ráster R, puede usar Rfunciones estándar o el getValuesmétodo para acceder a los valores de las celdas. A partir de ahí, es sencillo identificar los valores más altos y sus ubicaciones.
whuber
1
Gracias a tu recomendación, esto es lo que terminé haciendo. Usar el paquete ráster en R ha sido muy fácil en comparación con probar esto en ArcGIS. También he seguido usando otros análisis espaciales en R y he estado muy contento con los resultados. ¡Buen consejo!
mga
8

La respuesta puede ser obtenida por la combinación de una rejilla de indicador de la parte superior 1% de los valores con rejillas de latitud y longitud. El truco radica en crear esta cuadrícula de indicadores, porque ArcGIS (¡todavía! ¡Después de 40 años!) No tiene un procedimiento para clasificar los datos ráster.

Una solución para los rásteres de punto flotante es iterativa, pero afortunadamente rápida . Sea n el número de celdas de datos. La distribución acumulativa empírica de valores consiste en todos los pares (z, n (z)) donde z es un valor en la cuadrícula yn (z) es el número de celdas en la cuadrícula con valores menores o iguales a z . Obtenemos una curva que conecta (-infinito, 0) a (+ infinito, n) fuera de la secuencia de estos vértices ordenada por z . De este modo, define una función f , donde (z, f (z)) siempre se encuentra en la curva. Desea encontrar un punto (z0, 0.99 * n) en esta curva.

En otras palabras, la tarea es encontrar un cero de f (z) - (1-0.01) * n . Haga esto con cualquier rutina de búsqueda cero (que puede manejar funciones arbitrarias: esta no es diferenciable). El más simple, que a menudo es eficiente, es adivinar y verificar: inicialmente sabe que z0 se encuentra entre el valor mínimo zMin y el máximo zMax. Adivina cualquier valor razonable estrictamente entre estos dos. Si la suposición es demasiado baja, establezca zMin = z0; de lo contrario, establezca zMax = z0. Ahora repite. Usted convergerá rápidamente a la solución; estás lo suficientemente cerca cuando zMax y zMin están lo suficientemente cerca. Para ser conservador, elija el valor final de zMin como la solución: podría recoger algunos puntos adicionales que puede descartar más tarde. Para enfoques más sofisticados, vea el Capítulo 9 de Recetas Numéricas (el enlace va a una versión gratuita anterior).

Al mirar hacia atrás en este algoritmo, se revela que necesita realizar solo dos tipos de operaciones ráster : (1) seleccionar todas las celdas menores o iguales a algún valor objetivo y (2) contar las celdas seleccionadas. Esas son algunas de las operaciones más simples y rápidas. (2) se puede obtener como un recuento zonal o leyendo un registro de la tabla de atributos de la cuadrícula de selección.

whuber
fuente
7

Lo hice hace algún tiempo, aunque mi solución está usando GDAL (por lo tanto, esto no es solo para ArcGIS). Creo que puede obtener una matriz NumPy de un ráster en ArcGIS 10, pero no estoy seguro. NumPy proporciona indexación de matriz simple y potente, como argsorty otros. Este ejemplo no maneja NODATA ni transforma coordenadas de proyectado a lat / long (pero esto no es difícil de hacer con osgeo.osr, provisto con GDAL)

import numpy as np
from osgeo import gdal

# Open raster file, and get GeoTransform
rast_src = gdal.Open(rast_fname)
rast_gt = rast_src.GetGeoTransform()

def get_xy(r, c):
    '''Get (x, y) raster centre coordinate at row, column'''
    x0, dx, rx, y0, ry, dy = rast_gt
    return(x0 + r*dx + dx/2.0, y0 + c*dy + dy/2.0)

# Get first raster band
rast_band = rast_src.GetRasterBand(1)

# Retrieve as NumPy array to do the serious work
rast = rast_band.ReadAsArray()

# Sort raster pixels from highest to lowest
sorted_ind = rast.argsort(axis=None)[::-1]

# Show highest top 10 values
for ind in sorted_ind[:10]:
    # Get row, column for index
    r, c = np.unravel_index(ind, rast.shape)
    # Get [projected] X and Y coordinates
    x, y = get_xy(r, c)
    print('[%3i, %3i] (%.3f, %.3f) = %.3f'%
          (r, c, x, y, rast[r, c]))

Muestra lo siguiente para mi archivo ráster de prueba:

[467, 169] (2813700.000, 6353100.000) = 844.538
[467, 168] (2813700.000, 6353200.000) = 841.067
[469, 168] (2813900.000, 6353200.000) = 840.705
[468, 168] (2813800.000, 6353200.000) = 840.192
[470, 167] (2814000.000, 6353300.000) = 837.063
[468, 169] (2813800.000, 6353100.000) = 837.063
[482, 166] (2815200.000, 6353400.000) = 833.038
[469, 167] (2813900.000, 6353300.000) = 832.825
[451, 181] (2812100.000, 6351900.000) = 828.064
[469, 169] (2813900.000, 6353100.000) = 827.514
Mike T
fuente
+1 Gracias por compartir esto. No veo una incapacidad para manejar NoData como una limitación: simplemente convierta todos los NoData a valores extremadamente negativos antes de continuar. Tenga en cuenta, también, que si se produce una reproyección de la cuadrícula , las respuestas probablemente cambiarán debido al remuestreo de la cuadrícula, por lo que normalmente no se desea que la reproyección ocurra automáticamente durante un cálculo como este. En cambio, las coordenadas informadas pueden reproyectarse después. Por lo tanto, su solución es perfectamente general.
whuber
El manejo de NODATA se puede implementar obteniendo primero el valor del ráster NODATA = rast_band.GetNoDataValue(), luego usando un valor NaN ( rast[rast == NODATA] = np.nan) o usando una matriz enmascarada ( rast = np.ma.array(rast, mask=(rast == NODATA))). El truco más complicado es argsorteliminar de alguna manera los valores NODATA del análisis, o simplemente pasarlos por alto en el bucle for si están enmascarados o con NaN.
Mike T