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

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?

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ą.

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):

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

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()).

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:

  

• Przefiltrowany:

  

Ł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.