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.
Respuestas:
Si está contento de usar R , hay un paquete llamado raster . Puede leer un ráster con el siguiente comando:
Luego, cuando lo veas (escribiendo
test
), puedes ver la siguiente información: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.
fuente
R
, puede usarR
funciones estándar o elgetValues
método para acceder a los valores de las celdas. A partir de ahí, es sencillo identificar los valores más altos y sus ubicaciones.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.
fuente
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
argsort
y 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)Muestra lo siguiente para mi archivo ráster de prueba:
fuente
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 esargsort
eliminar 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.