Jakie narzędzia do wygładzania / generowania rastra są dostępne?

46

Mam DEM, który chciałbym wygładzić lub uogólnić w celu usunięcia ekstremów topograficznych (odciąć szczyty i wypełnić doliny). Idealnie chciałbym mieć kontrolę nad promieniem lub poziomem „rozmycia”. Na koniec będę potrzebować zestawu rastrów od lekko rozmytych do naprawdę rozmytych. (Teoretycznie najbardziej rozmyty byłby stały raster średniej arytmetycznej wszystkich wartości).

Czy są jakieś narzędzia lub metody, których mogę użyć (w oparciu o Esri, GDAL, GRASS)? Czy muszę upiec w domu własną rutynową rozmycie Gaussa ? Czy mogę zastosować filtr dolnoprzepustowy (np. Filtr ArcGIS ), a jeśli tak, czy musiałbym go uruchamiać kilka razy, aby uzyskać efekt dużego promienia?

Mike T.
źródło
A co z eksportowaniem rastra do większego rozmiaru komórki? Czy nie doprowadziłoby to również do wyciszenia skrajności?
1
Tak, zmniejszyłoby to również skrajności (zakładając, że niejawne ponowne próbkowanie wymaga pewnej formy uśredniania), ale jest to okropny sposób na wygładzenie DEM: utworzyłbyś niewielką liczbę dużych bloków. BTW, zwykle nie trzeba eksportować rastra, aby to zrobić; agregacja, a także ponowne próbkowanie do innej wielkości komórki to podstawowe operacje zwykle spotykane w oprogramowaniu rastrowym.
whuber

Odpowiedzi:

29

Rozmycie gaussowskie to tylko ważona średnia ogniskowa. Możesz go odtworzyć z dużą dokładnością dzięki sekwencji krótkiego sąsiedztwa kołowego (nieważonego): oznacza to zastosowanie Twierdzenia o granicy centralnej .

Masz duży wybór. „Filtr” jest zbyt ograniczony - dotyczy tylko dzielnic 3 x 3 - więc nie przejmuj się nim. Najlepszą opcją dla dużych DEM jest przeniesienie obliczeń poza ArcGIS do środowiska, które korzysta z szybkich transformacji Fouriera: wykonują te same obliczenia ogniskowe, ale (w porównaniu) robią to niesamowicie szybko. (GRASS ma moduł FFT . Jest przeznaczony do przetwarzania obrazu, ale możesz go wypchnąć do pracy dla swojego DEM, jeśli możesz przeskalować go z rozsądną precyzją do zakresu 0..255.) Pomijając to, przynajmniej dwa rozwiązania są warte rozważenia:

  1. Utwórz zestaw wag sąsiedztwa, aby przybliżyć rozmycie gaussowskie dla sporego sąsiedztwa. Użyj kolejnych przejść tego rozmycia, aby utworzyć sekwencję coraz płynniejszych DEM.

    (Wagi są obliczane jako exp (-d ^ 2 / (2r)), gdzie d to odległość (w komórkach, jeśli chcesz), a r to promień efektywny (również w komórkach). do co najmniej 3r . Po wykonaniu tej czynności podziel każdą wagę przez ich sumę, aby na koniec sumować do 1)

  2. Ewentualnie zapomnij o ważeniu; po prostu kilkakrotnie użyj okrągłej wartości ogniskowej. Zrobiłem to dokładnie w celu zbadania, jak siatki pochodne (takie jak nachylenie i aspekt) zmieniają się wraz z rozdzielczością DEM.

Obie metody będą działać dobrze, a po pierwszych kilku przejściach nie będzie już więcej do wyboru między tymi dwoma, ale maleją zwroty: efektywny promień n kolejnych środków ogniskowych (wszystkie wykorzystujące ten sam rozmiar sąsiedztwa) jest tylko (w przybliżeniu) pierwiastek kwadratowy n razy promień średniej ogniskowej. Tak więc, w przypadku ogromnego rozmycia, będziesz chciał zacząć od nowa z sąsiedztwem o dużym promieniu. Jeśli użyjesz nieważonej średniej ogniskowej, przeprowadź 5-6 przejść nad DEM. Jeśli używasz odważników zbliżonych do gaussowskiego, potrzebujesz tylko jednego przejścia: ale musisz utworzyć macierz wag.

To podejście rzeczywiście ma średnią arytmetyczną DEM jako wartość graniczną.

Whuber
źródło
1
Jeśli twoje dane mają skoki, możesz najpierw wypróbować filtr mediany ( en.wikipedia.org/wiki/Median_filter ), zanim zastosujesz bardziej ogólne rozmycie, jak sugeruje whuber.
MerseyViking
@Mersey To doskonała sugestia. Nigdy nie widziałem DEM z lokalnymi wartościami odstającymi, ale z drugiej strony nigdy nie musiałem przetwarzać surowego DEM (takiego jak surowe wyniki LIDAR). Nie możesz zrobić filtrów mediany za pomocą FFT, ale potrzebujesz tylko (zazwyczaj) sąsiedztwa 3 x 3, więc i tak jest to szybka operacja.
whuber
Dzięki, kurwa. Muszę przyznać, że używałem tylko wstępnie przetworzonych danych LiDAR, ale istnieją pewne znaczące wzrosty w danych SRTM, które skorzystałyby na filtrze mediany. Zwykle mają szerokość 2 lub 3 próbek, więc potrzebny byłby większy filtr mediany.
MerseyViking
@Mersey Nadal jesteś w porządku z większym medianą filtra 5 x 5 lub 7 x 7. Jeśli zastanawiasz się (powiedzmy) filtr 101 x 101, przygotuj się jednak na czekanie! Sugerujesz także ważną kwestię, którą warto rozwinąć: bardzo dobrze jest przeprowadzić analizę eksploracyjną DEM przed zrobieniem czegokolwiek. Obejmowałoby to identyfikację szczytów (lokalne wartości odstające) oraz charakteryzowanie ich wielkości i zakresu. Chcesz się upewnić, że to naprawdę artefakty (a nie jakieś prawdziwe zjawisko), zanim zaczniesz wymazywać je za pomocą filtra!
whuber
1
+1 dla FFT na danych wysokości. Naprawdę sprawiłem, że działał on na trawie dla 32-bitowych danych NED, aby usunąć dwukierunkowe paski. W końcu było to również problematyczne, ponieważ ponownie wprowadzono efekt tarasowania, który nęka wiele innych DEM pochodzących z konturów.
Jay Guarneri
43

Eksplorowałem podejście Signal.convolve SciPy (oparte na tej książce kucharskiej ) i odnoszę naprawdę niezły sukces dzięki poniższemu fragmentowi:

import numpy as np
from scipy.signal import fftconvolve

def gaussian_blur(in_array, size):
    # expand in_array to fit edge of kernel
    padded_array = np.pad(in_array, size, 'symmetric')
    # build kernel
    x, y = np.mgrid[-size:size + 1, -size:size + 1]
    g = np.exp(-(x**2 / float(size) + y**2 / float(size)))
    g = (g / g.sum()).astype(in_array.dtype)
    # do the Gaussian blur
    return fftconvolve(padded_array, g, mode='valid')

Używam tego w innej funkcji, która odczytuje / zapisuje pliki GeoTIFF typu float32 przez GDAL (nie trzeba przeskalowywać do 0-255 bajtów w celu przetworzenia obrazu), i próbowałem użyć wielkości pikseli (np. 2, 5, 20) i ma naprawdę ładne wyjście (wizualizowane w ArcGIS z pikselem 1: 1 i stałym zakresem min / max):

Gaussian DTM

Uwaga: ta odpowiedź została zaktualizowana, aby używać znacznie szybszej funkcji przetwarzania sygnału opartego na FFT .fftconvolve .

Mike T.
źródło
1
+1 Dobre rozwiązanie! Nie wiem na pewno, ale dobrze jest założyć, że signal.convolve używa FFT.
whuber
Szukałem rozmazanego kodu dla narzędzia do automatycznego łączenia, które piszę i natknąłem się na to. Dobra robota @MikeToews!
Ragi Yaser Burhum
@RagiYaserBurhum Chciałbym dowiedzieć się więcej o swoim narzędziu. MikeToews Świetna odpowiedź i doceniony fragment kodu.
Jay Laura,
@JayLaura Nic specjalnego, tylko napisanie narzędzia do automatycznego łączenia niektórych zdjęć, które zrobiłem z przyjaciółmi z balonem. Korzystanie z klas Orfeo Toolbox orfeo-toolbox.org/SoftwareGuide/…
Ragi Yaser Burhum
2
@ Whuber po przejrzeniu tej procedury nie korzystał z FFT, ale jest teraz i jest o wiele szybszy.
Mike T
4

Może to być komentarz do doskonałej odpowiedzi MikeT , jeśli nie byłby zbyt długi i zbyt skomplikowany. Dużo się z tym bawiłem i stworzyłem wtyczkę QGIS o nazwie FFT Convolution Filters (na etapie „eksperymentalnym”) w oparciu o jego funkcję. Oprócz wygładzania, wtyczka może również wyostrzać krawędzie, odejmując wygładzony raster od oryginalnego.

Trochę zaktualizowałem funkcję Mike'a:

def __gaussian_blur1d(self, in_array, size):
        #check validity
        try:
            if 0 in in_array.shape:
                raise Exception("Null array can't be processed!")
        except TypeError:
            raise Exception("Null array can't be processed!")
        # expand in_array to fit edge of kernel
        padded_array = np.pad(in_array, size, 'symmetric').astype(float)
        # build kernel
        x, y = np.mgrid[-size:size + 1, -size:size + 1]
        g = np.exp(-(x**2 / float(size) + y**2 / float(size)))
        g = (g / g.sum()).astype(float)
        # do the Gaussian blur
        out_array = fftconvolve(padded_array, g, mode='valid')
        return out_array.astype(in_array.dtype)

Kontrole ważności są dość oczywiste, ale ważne jest rzutowanie na pływak i powrót. Wcześniej funkcja zamieniała tablice liczb całkowitych na czarne (tylko zera), ze względu na podzielenie przez sumę wartości ( g / g.sum()).

Pavel V.
źródło
3

W QGIS łatwo uzyskałem dobre wyniki dzięki filtrowaniu obrazów Orfeo Toolbox . Jest dość szybki, a tryb wsadowy działa dobrze. Dostępne są dyfuzje Gaussa, średnie lub anizotropowe.

Uwaga: Radiusodnosi się do liczby komórek, a nie odległości.

Oto przykład użycia wygładzania (gaussa) :

  • Surowy:

    Bez filtra

  • Przefiltrowany:

    filtr

Tactopoda
źródło
1

Ładne rozwiązanie dla rozmycia Gaussa i fajnej animacji. Jeśli chodzi o wspomniane powyżej narzędzie Esri Filter, jest to po prostu narzędzie Esri „Focal Statistics” zakodowane na stałe do rozmiaru 3x3. Narzędzie Ogniskowe statystyki daje o wiele więcej opcji dotyczących kształtu ruchomego filtra, rozmiaru i statystyki, którą chcesz uruchomić. http://desktop.arcgis.com/en/arcmap/latest/tools/spatial-analyst-toolbox/focal-statistics.htm

Możesz także utworzyć „nieregularny” filtr, w którym przekazujesz swój własny plik tekstowy z wagami do użycia dla każdej komórki. Plik tekstowy ma tyle wierszy, ile chcesz w obszarze filtru, z wartościami oddzielonymi spacjami dla kolumn. Myślę, że zawsze powinieneś używać nieparzystej liczby wierszy i kolumn, więc twoja komórka docelowa jest w środku.

Utworzyłem arkusz kalkulacyjny programu Excel, aby grać z różnymi wagami, które po prostu kopiuję / wklejam do tego pliku. Powinien osiągnąć takie same wyniki jak powyżej, jeśli dostosujesz formuły.

David A.
źródło