Welche Tools zur Rasterglättung / Generalisierung sind verfügbar?

46

Ich habe ein DEM, das ich glätten oder verallgemeinern möchte, um topografische Extreme zu entfernen (Spitzen abschneiden und Täler auffüllen). Idealerweise hätte ich auch gerne die Kontrolle über den Radius oder den Grad der "Unschärfe". Am Ende brauche ich eine Reihe von Rastern, die von leicht verschwommen bis sehr verschwommen reichen. (Theoretisch wäre die Unschärfe ein konstantes Raster des arithmetischen Mittels aller Werte.)

Gibt es Tools oder Methoden, die ich verwenden kann (basierend auf Esri, GDAL, GRASS)? Muss ich zu Hause meine eigene Gaußsche Unschärferoutine backen ? Könnte ich einen Tiefpassfilter (z. B. den ArcGIS- Filter ) verwenden und müsste er in diesem Fall einige Male ausgeführt werden, um den Effekt eines großen Radius zu erzielen?

Mike T
quelle
Wie wäre es, wenn Sie nur das Raster in eine größere Zelle exportieren? Würde dies nicht auch zu einer Stummschaltung der Extreme führen?
1
Ja, das würde auch Extreme reduzieren (vorausgesetzt, dass das implizite Resampling eine Art Mittelwertbildung beinhaltet), aber es ist ein schrecklicher Weg, ein DEM zu glätten : Sie würden eine kleine Anzahl großer Blöcke erstellen. Übrigens muss man normalerweise kein Raster exportieren, um dies zu tun; Aggregation sowie Resampling auf eine andere Zellengröße sind grundlegende Vorgänge, die normalerweise in rasterbasierter Software zu finden sind.
Whuber

Antworten:

29

Die Gaußsche Unschärfe ist nur ein gewichteter Mittelwert. Sie können es mit einer Folge von kreisförmigen Nahbereichen (ungewichtet) mit hoher Genauigkeit nachbilden: Dies ist eine Anwendung des zentralen Grenzwertsatzes .

Sie haben viele Möglichkeiten. "Filter" ist zu begrenzt - es ist nur für 3 x 3 Nachbarschaften - also kümmere dich nicht darum. Die beste Option für große DEMs besteht darin, die Berechnung außerhalb von ArcGIS in eine Umgebung zu übertragen, in der schnelle Fouriertransformationen verwendet werden: Sie führen dieselben Fokusberechnungen durch, aber (im Vergleich) sind sie blitzschnell. (GRASS verfügt über ein FFT-Modul . Es ist für die Bildverarbeitung vorgesehen, Sie können es jedoch möglicherweise für Ihr DEM in Betrieb nehmen, wenn Sie es mit angemessener Genauigkeit in den Bereich 0..255 skalieren können.) Wenn dies nicht möglich ist, gibt es mindestens zwei Lösungen bedenkenswert:

  1. Erstellen Sie eine Reihe von Nachbarschaftsgewichten, um eine Gaußsche Unschärfe für eine beträchtliche Nachbarschaft zu approximieren. Verwenden Sie aufeinanderfolgende Durchläufe dieser Unschärfe, um eine Sequenz immer glatterer DEMs zu erstellen.

    (Die Gewichte werden als exp (-d ^ 2 / (2r)) berechnet, wobei d der Abstand (in Zellen, wenn Sie möchten) und r der effektive Radius (auch in Zellen) ist. Sie müssen innerhalb eines sich erstreckenden Kreises berechnet werden bis mindestens 3r dividieren. Danach dividieren Sie jedes Gewicht durch die Summe aller, so dass am Ende die Summe 1 ergibt .)

  2. Alternativ vergessen Sie die Gewichtung; Führen Sie einfach wiederholt einen kreisförmigen Mittelwert aus. Ich habe genau dies getan, um zu untersuchen, wie sich abgeleitete Gitter (wie Neigung und Aspekt) mit der Auflösung eines DEM ändern.

Beide Methoden werden gut funktionieren, und nach den ersten Durchgängen wird es wenig geben, um zwischen den beiden zu wählen, aber es gibt abnehmende Renditen: Der effektive Radius von n aufeinanderfolgenden Fokusmitteln (alle mit derselben Nachbarschaftsgröße) ist nur (ungefähr) der Quadratwurzel des n- fachen Radius des Mittelwerts. Bei großen Unschärfen sollten Sie daher mit einer Nachbarschaft mit großem Radius von vorne beginnen. Wenn Sie einen ungewichteten Mittelwert verwenden, durchlaufen Sie den DEM in 5-6 Durchläufen. Wenn Sie Gewichte verwenden, die ungefähr Gauß'sch sind, benötigen Sie nur einen Durchgang. Sie müssen jedoch die Gewichtsmatrix erstellen.

Dieser Ansatz hat in der Tat das arithmetische Mittel des DEM als Grenzwert.

whuber
quelle
1
Wenn Ihre Daten Spitzen aufweisen, können Sie zuerst einen Medianfilter ( en.wikipedia.org/wiki/Median_filter ) ausprobieren, bevor Sie eine allgemeinere Unschärfe anwenden, wie von whuber vorgeschlagen.
MerseyViking
@ Jersey Das ist ein ausgezeichneter Vorschlag. Ich habe noch nie ein DEM mit lokalen Ausreißern gesehen, musste aber auch noch nie ein unformatiertes DEM (wie z. B. unformatierte LIDAR-Ergebnisse) verarbeiten. Sie können mit FFT keine Medianfilter erstellen, benötigen jedoch (normalerweise) nur eine 3 x 3-Umgebung, sodass die Operation schnell ausgeführt werden kann.
whuber
Vielen Dank. Ich muss zugeben, dass ich bisher nur vorverarbeitete LiDAR-Daten verwendet habe, aber es gibt einige signifikante Spitzen in SRTM-Daten, die von einem Medianfilter profitieren würden. Sie sind jedoch in der Regel 2 oder 3 Samples breit, sodass ein größerer Medianfilter erforderlich wäre.
MerseyViking
@Mersey Ein größerer Medianfilter von 5 x 5 oder 7 x 7 ist immer noch in Ordnung. Wenn Sie jedoch über einen 101 x 101-Filter nachdenken (sagen wir), sollten Sie warten! Sie schlagen auch einen wichtigen Punkt vor, der näher erläutert werden sollte: Es ist eine sehr gute Idee, eine explorative Analyse des DEM durchzuführen, bevor Sie etwas unternehmen. Dies würde das Identifizieren von Spitzen (lokalen Ausreißern) und das Charakterisieren ihrer Größen und Ausmaße umfassen. Sie möchten sichergehen, dass es sich um echte Artefakte handelt (und nicht um ein echtes Phänomen), bevor Sie sie mit einem Filter auslöschen!
whuber
1
+1 für FFT auf Höhendaten. Ich habe diese Arbeit im Gras für 32-Bit-NED-Daten gemacht, um bidirektionales Striping zu entfernen. Dies war letztendlich auch problematisch, weil der Terrassierungseffekt, der viele andere von Konturen abgeleitete DEMs plagt, wieder eingeführt wurde.
Jay Guarneri
43

Ich habe den signal.convolve- Ansatz von SciPy (basierend auf diesem Kochbuch ) untersucht und habe mit dem folgenden Snippet einige wirklich schöne Erfolge erzielt:

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

Ich verwende dies in einer anderen Funktion, die float32-GeoTIFFs über GDAL liest / schreibt (für die Bildverarbeitung muss keine Neuskalierung auf 0-255 Byte vorgenommen werden), und ich habe versucht, Pixelgrößen (z. B. 2, 5, 20) zu verwenden, und dies ist der Fall Sehr schöne Ausgabe (in ArcGIS mit 1: 1 Pixel und konstantem Min / Max-Bereich dargestellt):

Gaußsche DTM

Hinweis: Diese Antwort wurde aktualisiert, um eine viel schnellere FFT-basierte Signal - Fft - Convolve- Verarbeitungsfunktion zu verwenden.

Mike T
quelle
1
+1 Schöne Lösung! Ich weiß es nicht genau, aber es ist eine gute Wette, dass signal.convolve FFTs verwendet.
whuber
Ich war auf der Suche nach einem Code für ein automatisches Stitching-Tool, das ich schreibe, und bin darauf gestoßen. Gute Arbeit @MikeToews!
Ragi Yaser Burhum
@RagiYaserBurhum Würde gerne mehr über Ihr Werkzeug erfahren. MikeToews Tolle Antwort und viel geschätzter Code-Ausschnitt.
Jay Laura
@JayLaura Nichts Besonderes, ich schreibe nur ein Tool, um einige Bilder, die ich mit einigen Freunden mit einem Ballon aufgenommen habe, automatisch zusammenzufügen. Verwenden der Orfeo Toolbox-Klassen orfeo-toolbox.org/SoftwareGuide/…
Ragi Yaser Burhum
2
@whuber nach Überarbeitung dieser Routine, es war nicht mit FFT, aber es ist jetzt und ist so viel schneller.
Mike T
4

Dies könnte ein Kommentar zu MikeTs hervorragender Antwort sein , wenn es nicht zu lang und zu komplex wäre. Ich habe viel damit gespielt und ein QGIS-Plugin namens FFT Convolution Filters (noch im "experimentellen" Stadium) basierend auf seiner Funktion erstellt. Neben dem Glätten kann das Plug-in auch Kanten schärfen, indem das geglättete Raster vom ursprünglichen subtrahiert wird.

Ich habe Mikes Funktion ein wenig verbessert:

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)

Die Gültigkeitsprüfungen sind ganz selbstverständlich, aber wichtig ist, dass Sie schweben und zurück. Zuvor hat die Funktion Integer-Arrays aufgrund der Division durch die Summe der Werte ( g / g.sum()) schwarz (nur Nullen ) gemacht.

Pavel V.
quelle
3

In QGIS habe ich mit der Orfeo Toolbox- Bildfilterung problemlos gute Ergebnisse erzielt . Es ist relativ schnell und der Batch-Modus funktioniert einwandfrei. Gaußsche, mittlere oder anisotrope Diffusionen sind verfügbar.

Beachten Sie, dass sich die RadiusAngabe auf die Anzahl der Zellen bezieht, nicht auf die Entfernung.

Hier ist ein Beispiel für die Verwendung von Glättung (Gauß) :

  • Roh:

    Kein Filter

  • Gefiltert:

    Filter

Tactopoda
quelle
1

Schöne Lösung für die Gaußsche Unschärfe und coole Animation. In Bezug auf das oben erwähnte Esri-Filter-Tool handelt es sich im Grunde genommen nur um das Esri-Tool "Focal Statistics", das auf eine Größe von 3 x 3 fest codiert ist. Mit dem Werkzeug "Fokusstatistik" haben Sie viel mehr Optionen für die Form Ihres beweglichen Filters, die Größe und die Statistik, die Sie ausführen möchten. http://desktop.arcgis.com/de/arcmap/latest/tools/spatial-analyst-toolbox/focal-statistics.htm

Sie können auch einen "unregelmäßigen" Filter erstellen, in dem Sie Ihre eigene Textdatei mit Gewichten übergeben, die für jede Zelle verwendet werden. Die Textdatei enthält beliebig viele Zeilen in Ihrem Filterbereich mit durch Leerzeichen getrennten Werten für die Spalten. Ich denke, Sie sollten immer eine ungerade Anzahl von Zeilen und Spalten verwenden, damit sich Ihre Zielzelle in der Mitte befindet.

Ich habe eine Excel-Tabelle erstellt, um mit verschiedenen Gewichten zu spielen, die ich einfach kopiere / in diese Datei einfüge. Es sollte die gleichen Ergebnisse wie oben erzielen, wenn Sie die Formeln anpassen.

David A
quelle