Znaleźć lokalizacje o najwyższych wartościach w rastrze za pomocą ArcGIS Desktop?

Korzystając z ArcGIS 10, mam raster, w którym chciałbym znaleźć piksel o maksymalnej wartości w rastrze i zwrócić jego lokalizację (środek piksela) w stopniach dziesiętnych. Chciałbym powtórzyć ten proces, zwracając lokalizację drugiej najwyższej wartości rastra, a następnie trzeciej i tak dalej, tak że w końcu mam listę N lokalizacji, które mają najwyższe wartości w rastrze w kolejności.

Wyobrażam sobie, że najłatwiej można to zrobić za pomocą skryptu Python, ale jestem otwarty na inne pomysły, jeśli istnieje lepszy sposób.

Jeśli z przyjemnością używasz R , istnieje pakiet o nazwie raster . Możesz czytać w rastrze za pomocą następującego polecenia:

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

Następnie, gdy przejrzysz go (wpisując test), zobaczysz następujące informacje:

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 

Mogą istnieć lepsze sposoby manipulowania rastrem, ale jednym ze sposobów znalezienia potrzebnych informacji może być znalezienie najwyższej wartości i uzyskanie jej położenia macierzy, a następnie dodanie jej do niższych zakresów.

Odpowiedź można uzyskać , łącząc siatkę wskaźnikową z górnego 1% wartości z siatkami szerokości i długości geograficznej. Sztuczka polega na stworzeniu tej siatki wskaźników, ponieważ ArcGIS (wciąż! Po 40 latach!) Nie ma procedury uszeregowania danych rastrowych.

Jedno rozwiązanie dla rastrów zmiennoprzecinkowych jest iteracyjne, ale na szczęście szybkie . Niech n będzie liczbą komórek danych. Empiryczny skumulowany rozkład wartości obejmuje wszystkie pary (Z, n (z)), w którym z ma wartość w siatce i n (z) jest liczbą komórek w sieci o wartości mniejszej niż lub równej Z . Otrzymujemy krzywą łączącą (-nieskończoność, 0) z (+ nieskończoność, n) z sekwencji tych wierzchołków uporządkowanych według z . Definiuje w ten sposób funkcję f , gdzie (z, f (z)) zawsze leży na krzywej. Chcesz znaleźć punkt (z0, 0,99 * n) na tej krzywej.

Innymi słowy, zadaniem jest znalezienie zera f (z) - (1-0.01) * n . Zrób to za pomocą dowolnej procedury znajdowania zera (która może obsłużyć dowolne funkcje: tej nie można rozróżnić). Najprostszym, często wydajnym, jest odgadywanie i sprawdzanie: początkowo wiadomo, że z0 leży między wartością minimalną zMin a maksymalną zMax. Zgadnij jakąkolwiek rozsądną wartość ściśle między tymi dwoma. Jeśli prawdopodobieństwo jest zbyt niskie, ustaw zMin = z0; w przeciwnym razie ustaw zMax = z0. Teraz powtórz. Szybko przejdziesz do rozwiązania; jesteś wystarczająco blisko, gdy zMax i zMin są wystarczająco blisko. Aby zachować ostrożność, wybierz ostateczną wartość zMin jako rozwiązanie: może zebrać kilka dodatkowych punktów, które możesz odrzucić później. Bardziej wyrafinowane podejście znajduje się w rozdziale 9 przepisów numerycznych (link prowadzi do starszej bezpłatnej wersji).

Patrząc wstecz na ten algorytm, ujawniasz, że musisz wykonać tylko dwa rodzaje operacji rastrowych : (1) zaznacz wszystkie komórki mniejsze lub równe pewnej wartości docelowej i (2) policz wybrane komórki. Są to jedne z najprostszych i najszybszych operacji na rynku. (2) można uzyskać jako liczbę strefową lub przez odczytanie jednego rekordu z tabeli atrybutów siatki wyboru.

Zrobiłem to jakiś czas temu, chociaż moje rozwiązanie korzysta z GDAL (więc nie dotyczy to tylko ArcGIS). Myślę, że możesz uzyskać tablicę NumPy z rastra w ArcGIS 10, ale nie jestem pewien. NumPy zapewnia proste i wydajne indeksowanie tablic, takie jak argsorti inne. Ten przykład nie obsługuje NODATA ani nie przekształca współrzędnych z rzutowanych na długość / długość (ale nie jest to trudne w przypadku osgeo.osr, dostarczonego z 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]))

Pokazuje mój testowy plik rastrowy:

[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