Jak buforować piksele rastrowe według ich wartości?

28

Piksele po lewej stronie reprezentują lokalizacje drzewa i związane z nimi promienie korony (tj. Wartości pikseli w zakresie od 2 do 5). Chciałbym buforować te piksele rastrowe według ich wartości promienia korony. Obraz po prawej stronie jest tym, co mam nadzieję osiągnąć przy użyciu tylko metod przetwarzania rastrowego .

Początkowo pomyślałbym o zastosowaniu okrągłej sumy ogniskowej w ArcGIS, chociaż ustawienie sąsiedztwa jest stałą wartością, która nie uwzględniałaby promienia korony o zmiennej wielkości.

Jaka jest dobra metoda „buforowania” pikseli według ich wartości?

wprowadź opis zdjęcia tutaj

Aaron
źródło
2
Czy próbowałeś przekonwertować raster na punkty, następnie buforować według pola, a następnie przekonwertować z powrotem na raster?
2
Pomaga uświadomić sobie, że jest to operacja nielokalna , ponieważ ta nielokalność pokazuje, że istnieją nieodłączne ograniczenia dotyczące sposobu wykonywania obliczeń. Na przykład twój wynik zmieniłby się radykalnie prawie wszędzie, gdyby tylko jeden izolowany piksel na wejściu miałby zmienić się na dużą wartość. Dlatego jeśli znasz ograniczenia wartości wejściowych, udostępnij je, ponieważ może to prowadzić do ulepszonych rozwiązań. Na przykład, czy wszystkie wartości wejściowe zawsze będą znajdować się w zbiorze {2,3,4}?
whuber
@ Dan Patterson W ten sposób wymyśliłem obraz po prawej stronie. Staram się jednak całkowicie unikać operacji wektorowych i tych kroków.
Aaron
2
@whuber Ten zestaw danych reprezentuje drzewa o zmiennej średnicy korony. Biorąc to pod uwagę, wymiary promienia korony drzewa mogą realistycznie różnić się od 1-10. Powinienem również dodać, że zbuforowane wyjście musi wynosić tylko 0 dla braku korony i 1 dla obecności korony.
Aaron
1
Ok dzięki. Wygląda na to, że utworzyłeś przykładowy wynik, łącząc 3 bufory punktów o wartości 3, 4 bufory punktów o wartości 4 i 5 buforów punktów o wartości 5. (Wygląda na to, że zapomniałeś przetwarzać punkty o wartości 2.) Ten proces nie tylko odpowiada na twoje pytanie, ale (wierzę) jest to najprostsze rozwiązanie z wykorzystaniem narzędzi dostępnych w Spatial Analyst.
whuber

Odpowiedzi:

14

Oto czyste rozwiązanie rastrowe w Python 2.7użyciu numpyi scipy:

import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt

#create tree location matrix with values indicating crown radius
A = np.zeros((120,320))
A[60,40] = 1
A[60,80] = 2
A[60,120] = 3
A[60,160] = 4
A[60,200] = 5
A[60,240] = 6
A[60,280] = 7

#plot tree locations
fig = plt.figure()
plt.imshow(A, interpolation='none')
plt.colorbar()

#find unique values
unique_vals = np.unique(A)
unique_vals = unique_vals[unique_vals > 0]

# create circular kernel
def createKernel(radius):
    kernel = np.zeros((2*radius+1, 2*radius+1))
    y,x = np.ogrid[-radius:radius+1, -radius:radius+1]
    mask = x**2 + y**2 <= radius**2
    kernel[mask] = 1
    return kernel

#apply binary dilation sequentially to each unique crown radius value 
C = np.zeros(A.shape).astype(bool)   
for k, radius in enumerate(unique_vals):  
    B = ndimage.morphology.binary_dilation(A == unique_vals[k], structure=createKernel(radius))
    C = C | B #combine masks

#plot resulting mask   
fig = plt.figure()
plt.imshow(C, interpolation='none')
plt.show()

Wkład: wprowadź opis zdjęcia tutaj

Wydajność: wprowadź opis zdjęcia tutaj

jatobat
źródło
1
+1 za podejście dylatacyjne! Działa również w pobliżu punktów.
Antonio Falciano
To świetny przykład tego, dlaczego ten stary odrzutowy barwnik był okropny. Z viridis wygląda to znacznie wyraźniej.
naught101
8

Podejście wektorowe

To zadanie można wykonać w trzech krokach:

Uwaga: użycie pola bufora pozwala uniknąć obliczenia bufora dla każdej wartości promienia korony.


Podejście oparte na rastrze

Unikając rozwiązania wektorowego, problem ten sugeruje zastosowanie rodzaju automatów komórkowych opartych na najbliższych sąsiadach. Zakładając, że wszystkie czarne piksele są zerami, piksele są podniesione do kwadratu, a ich rozmiar jest równy 1 (lub, alternatywnie, są odpowiednio skalowane), zasady do przyjęcia są bardzo proste:

  1. Jeśli wartość piksela ( VALUE) jest większa niż 1, jej wartość staje się, VALUE-1a następnie rozważ otaczające ją piksele. Jeśli ich wartości są mniejsze niż VALUE-1, piksele te rodzą się lub rosną, a ich wartość staje się VALUE-1. W przeciwnym razie piksele te przetrwają i pozostaną niezmienione.
  2. Jeśli VALUE<=1nic nie rób (piksel nie żyje!).

Reguły te muszą być stosowane, dopóki wszystkie piksele nie będą martwe, tzn. Ich wartości będą równe 0 lub 1. Więc N-1razy, gdzie Njest maksymalna wartość, którą masz w rastrze wejściowym. To podejście można dość łatwo zaimplementować za pomocą odrobiny Pythona i numpy.

Antonio Falciano
źródło
1
Dzięki za odpowiedź afalciano. W ten sposób stworzyłem obraz po prawej stronie i wykorzystuję podejście wektorowe - takiego, którego staram się unikać.
Aaron
1
Ok Aaron, oto teraz podejście oparte na rastrze. Mam nadzieję że to pomoże.
Antonio Falciano
7

Inną opcją byłoby utworzenie osobnych rastrów dla każdej wartości piksela, w tym przypadku 4 rastrów, z warunkiem. Następnie rozwiń rastry o liczbę pikseli odpowiadającą wartości rastra (ewentualnie poprzez iterację po liście wartości). Na koniec dołącz do rastrów (algebraicznych lub przestrzennych), aby utworzyć jeden binarny raster koron drzew.

HDunn
źródło
1
Ten pomysł jest właściwy. Szczegóły można poprawić: (1) wybór tworzy binarny (0,1) wskaźnik drzew o danym promieniu korony. (2) Suma ogniskowa tego wyboru - przy użyciu okrągłego sąsiedztwa o danym promieniu - jest szybka do obliczenia przy użyciu FFT. (3) Dodanie sum ogniskowych (punktowo) i porównanie z 0 daje pożądany bufor.
whuber
7

Jest to trudne pytanie w rastrze, ponieważ nie masz możliwości wykorzystania wartości piksela do zdefiniowania rozmiaru bufora. Dlatego musisz zrobić filtr ogniskowy dla każdej wartości, jak już powiedziałeś.

Oto możliwa odpowiedź, aby to zrobić za pomocą tylko 3 filtrów (nie mogłem znaleźć mniej), ale nie idealnie, jak wspomniał Whuber: twoje bufory zostaną obcięte, gdy drzewa będą blisko siebie.

1) EDYCJA: Alokacja euklidesowa (nie rozwiązuje to całkowicie problemu, ponieważ tnie bufory w pobliżu mniejszych drzew, ale jest lepsza niż artefakty z mojego pierwszego rozwiązania).

2) odległość euklidesowa wokół każdego piksela

3) kalkulator rastrowy (algebra mapy) z instrukcją warunkową

Con("allocation_result" > ("distance_result" / pixel_size) , 1 , 0)

Pamiętaj, że możesz dostosować warunki w zależności od potrzeb w zakresie promienia (z lub bez centralnego piksela)

radouxju
źródło
+1 To jest kreatywne podejście. Spróbuję sprawdzić, czy możliwe jest zwiększenie skali za pomocą tego podejścia.
Aaron
2
Euklidesowa metoda odległości nie działa, ponieważ oblicza tylko odległość do najbliższego drzewa, która niekoniecznie jest odległością do drzewa, którego korona pokrywa punkt.
whuber
2

Zastanawiasz się, dlaczego nie korzystasz z narzędzia rozwijania ArcGIS ?

import arcpy
from arcpy.sa import *

raster_in  = r'c:\test.tif'
raster_out = r'c:\test_out.tif'

outExpand1 = Expand(raster_in, 2, 2)
outExpand2 = Expand(outExpand1, 3, 3)
outExpand3 = Expand(outExpand3, 4, 4)
outExpand4 = Expand(outExpand4, 5, 5)

outExpand4.save(raster_out)

W przypadku nakładania się: najnowsze expandpolecenie obejmie poprzednie.

Mr. Che
źródło
2

Jeśli masz pozycję w pikselach, promień i algorytm okręgu punktu środkowego (wariant Bresenham Alg.) Daje ci wskazówkę. IMO łatwo jest stworzyć wielokąt z tego podejścia i myślę, że łatwo jest to zaimplementować w Pythonie. Połączenie tego zestawu wielokątów zapewnia obszar pokrycia.

huckfinn
źródło
Wiem, że to nie jest pytanie, ale czy chcesz dowiedzieć się więcej o prymitywach graficznych i wypełnieniu wielokąta linii skanowania? Dla cirles jest to bardzo łatwe. Wypukła kombinacja jest
modnym
Jak można to zastosować przy użyciu podstawowych operacji rastrowych?
whuber
Jeśli spróbujesz poradzić sobie z tym w przestrzeni rastrowej, określ punkty koła, sortuj je według y lub x i wypełnij przestrzeń prostą linią (linia skanowania), która jest w stanie wypełnić okrąg. W podejściu trójkątnym, jeśli budujesz okrąg przez przybliżenie sektorów trójgłowych i próbujesz wypełnić trójkąt, potrzebujesz testu, czy punkt jest wewnątrz lub na zewnątrz (kombinacja wypukła) i jest odwrotnie. A w podejściu „GIS” budowanie wielokątów (wielokąty zorientowane zgodnie z ruchem wskazówek zegara) i tworzenie związku jest trzecim (IMO najdroższym obliczeniowo).
huckfinn
Żeby było jasne: w podejściu „GIS” ... wykonaj operację algebraryczną, taką jak zjednoczenie, przecięcie, dotknięcie .... jest trzecim IMO najdroższym obliczeniowo.
huckfinn