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.
Odpowiedzi:
Jeśli z przyjemnością używasz R , istnieje pakiet o nazwie raster . Możesz czytać w rastrze za pomocą następującego polecenia:
Następnie, gdy przejrzysz go (wpisując
test
), zobaczysz następujące informacje: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.
źródło
R
możesz użyć standardowychR
funkcji lubgetValues
metody dostępu do wartości komórek. Stamtąd łatwo jest zidentyfikować najwyższe wartości i ich lokalizacje.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.
źródło
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
argsort
i 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)Pokazuje mój testowy plik rastrowy:
źródło
NODATA = rast_band.GetNoDataValue()
, a następnie za pomocą wartości NaN (rast[rast == NODATA] = np.nan
) lub tablicy zamaskowanej (rast = np.ma.array(rast, mask=(rast == NODATA))
). Bardziej skomplikowaną sztuczką jestargsort
jakoś usunięcie wartości NODATA z analizy lub po prostu pominięcie ich w pętli for, jeśli są one NaN / maskowane.