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

12

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.

mga
źródło
próbowałeś przekonwertować siatkę na punkty, a następnie dodać pola X, Y i sortować?
Jakub Sisak GeoGraphics
Czy wartości rastrowe są zmiennoprzecinkowe czy całkowite?
whuber
@Jakub - Nie, nie mam. Prawdopodobnie zainteresuje mnie tylko około 1% punktów, więc nie wiem, czy warto dodać pola x, y dla wszystkich punktów, a następnie sortować. Może jeśli nie ma bardziej wydajnej opcji?
MGA
@ whuber - Wartości rastrowe są zmiennoprzecinkowe.
MGA
@mga, warto spróbować. Konwersja jest dość szybka, a dodanie XY jest również domyślnym narzędziem. Usuwanie niechcianych rekordów jest prostą operacją i wszystkie mogą być powiązane w jeden model. Po prostu pomysł.
Jakub Sisak GeoGraphics

Odpowiedzi:

5

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.

djq
źródło
1
+1 Dla tego odniesienia. Po przeczytaniu rastra Rmożesz użyć standardowych Rfunkcji lub getValuesmetody dostępu do wartości komórek. Stamtąd łatwo jest zidentyfikować najwyższe wartości i ich lokalizacje.
whuber
1
Dzięki twojej rekomendacji właśnie to zrobiłem. Korzystanie z pakietu rastrowego w języku R było bardzo proste w porównaniu z wypróbowaniem tego w ArcGIS. Nadal korzystałem z innych analiz przestrzennych w R i byłem bardzo zadowolony z wyników. Dobra rada!
MGA
8

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.

Whuber
źródło
7

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
Mike T.
źródło
+1 Dziękujemy za udostępnienie tego. Nie widzę niemożności obsługi NoData jako ograniczenia: po prostu przekonwertuj wszystkie NoData na skrajnie ujemne wartości przed kontynuowaniem. Należy również zauważyć, że jeśli nastąpi jakiekolwiek ponowne rzutowanie siatki , odpowiedzi prawdopodobnie zmienią się z powodu ponownego próbkowania siatki, więc zwykle nie chce się, aby powtórzenie miało miejsce automatycznie podczas takich obliczeń. Zamiast tego zgłoszone współrzędne można później ponownie rzutować. Zatem twoje rozwiązanie jest całkowicie ogólne.
whuber
Obsługa NODATA może zostać zaimplementowana poprzez pobranie wartości z rastra 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ą jest argsortjakoś usunięcie wartości NODATA z analizy lub po prostu pominięcie ich w pętli for, jeśli są one NaN / maskowane.
Mike T