Quels sont les outils de lissage / généralisation disponibles?

46

J'ai un MNT que j'aimerais lisser ou généraliser pour éliminer les extrêmes topographiques (pics d'abattage et creux de remplissage). Idéalement, j'aimerais aussi pouvoir contrôler le rayon ou le niveau de "flou". À la fin, j'aurai besoin d'un ensemble de rasters allant du plus flou au plus flou. (Théoriquement, le plus flou serait un raster constant de la moyenne arithmétique de toutes les valeurs).

Existe-t-il des outils ou des méthodes que je peux utiliser (basés sur Esri, GDAL, GRASS)? Dois-je cuire à la maison ma propre routine de flou gaussien ? Pourrais-je utiliser un filtre passe-bas (par exemple, le filtre d'ArcGIS ), et si oui, devrais-je l'exécuter plusieurs fois pour obtenir l'effet d'un grand rayon?

Mike T
la source
Qu'en est-il de simplement exporter le raster vers une plus grande taille de cellule? Cela n'aboutirait-il pas également à une mise en sourdine des extrêmes?
1
Oui, cela réduirait également les extrêmes (en supposant que le rééchantillonnage implicite implique une certaine forme de calcul de la moyenne), mais c’est une très mauvaise façon de lisser un DEM: vous créez un petit nombre de gros blocs. En passant, il n’est généralement pas nécessaire d’exporter un raster pour le faire; L'agrégation ainsi que le ré-échantillonnage sur une taille de cellule différente sont des opérations de base que l'on trouve généralement dans les logiciels à base de raster.
whuber

Réponses:

29

Le flou gaussien est juste une moyenne focale pondérée. Vous pouvez le recréer avec une grande précision avec une séquence de voisinage circulaire à courte distance (non pondéré): il s’agit d’une application du théorème de la limite centrale .

Vous avez beaucoup de choix. Le "filtre" est trop limité - il ne concerne que les quartiers 3 x 3 - ne vous embêtez pas. La meilleure option pour les grands MNA consiste à intégrer le calcul en dehors d'ArcGIS à un environnement utilisant les transformations rapides de Fourier: ils effectuent les mêmes calculs focaux, mais (en comparaison), ils le font à une vitesse fulgurante. (GRASS possède un module FFT . Il est destiné au traitement des images, mais vous pourrez peut-être le mettre en service pour votre DEM si vous pouvez le redimensionner avec une précision raisonnable dans la plage 0..255.) À part cela, il existe au moins deux solutions. Vaut la peine d'être considéré:

  1. Créez un ensemble de pondérations de voisinage qui se rapprochent d'un flou gaussien pour un voisinage non négligeable. Utilisez les passes successives de ce flou pour créer votre séquence de MNT toujours plus lisses.

    (Les poids sont calculés sous la forme exp (-d ^ 2 / (2r)) où d est la distance (dans les cellules si vous préférez) et r le rayon effectif (également dans les cellules). Ils doivent être calculés dans un cercle sur au moins l' 3r . Après cela, diviser chaque poids par la somme de tous donc à la fin de leur somme 1.)

  2. Sinon, oubliez la pondération; juste exécuter une moyenne focale circulaire à plusieurs reprises. C'est exactement ce que j'ai fait pour étudier comment les grilles dérivées (comme la pente et l'aspect) changent avec la résolution d'un MNT.

Les deux méthodes fonctionneront bien et après les premiers passages, il y aura peu de choix, mais les rendements diminuent: le rayon effectif de n moyennes focales successives (toutes utilisant la même taille de voisinage) n'est que (approximativement) le plus proche. racine carrée de n fois le rayon de la moyenne focale. Ainsi, pour de grandes quantités de flou, vous voudrez recommencer avec un voisinage de grand rayon. Si vous utilisez une moyenne focale non pondérée, exécutez 5-6 passages sur le DEM. Si vous utilisez des poids approximativement gaussiens, vous n'avez besoin que d'un seul passage: vous devez créer la matrice de poids.

Cette approche a en effet la moyenne arithmétique du DEM comme valeur limite.

whuber
la source
1
Si vos données présentent des pics, vous pouvez commencer par essayer un filtre médian ( en.wikipedia.org/wiki/Median_filter ) avant d'appliquer un flou plus général, comme suggéré par whuber.
MerseyViking
@ Jersey C'est une excellente suggestion. Je n'ai jamais vu de DEM avec des valeurs aberrantes locales, mais je n'ai encore jamais eu à traiter de DEM brut (comme les résultats LIDAR bruts). Vous ne pouvez pas faire de filtres médians avec FFT, mais vous avez seulement (généralement) besoin d'un voisinage 3 x 3, donc c'est une opération rapide quand même.
whuber
Merci beaucoup. Je dois avouer que je n’ai jamais utilisé que des données LiDAR prétraitées, mais il existe quelques pics significatifs dans les données SRTM qui bénéficieraient d’un filtre médian. Cependant, ils ont tendance à avoir une largeur de 2 ou 3 échantillons. Un filtre médian plus grand serait donc nécessaire.
MerseyViking
@Mersey Vous êtes toujours d'accord avec un filtre médian plus grand de 5 x 5 ou 7 x 7. Si vous envisagez un filtre de 101 x 101, attendez-vous à attendre! Vous suggérez également un point important qui mérite d’être précisé: c’est une très bonne idée de procéder à une analyse exploratoire du DEM avant de faire quoi que ce soit. Cela inclurait l'identification des pics (points aberrants locaux) et la caractérisation de leurs tailles et de leurs étendues. Vous voulez vous assurer qu'ils sont vraiment des artefacts (et non un phénomène réel) avant de les éliminer avec un filtre!
whuber
1
+1 pour la FFT sur les données d'altitude. En fait, j'ai fait ce travail dans l'herbe pour les données NED 32 bits afin d'éliminer les bandes bidirectionnelles. Au final, cela a également posé problème, car il a réintroduit l’effet de terrassement qui sévit dans de nombreux autres MNA dérivés des courbes de niveau.
Jay Guarneri
43

J'explore l' approche signal.convolution de SciPy (basée sur ce livre de recettes ) et rencontre un vif succès avec l'extrait suivant:

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

J'utilise ceci dans une autre fonction qui lit / écrit des GeoTIFFs float32 via GDAL (inutile de redimensionner à 0-255 octets pour le traitement des images), et j'ai déjà essayé de tenter des tailles de pixels (par exemple, 2, 5, 20) et très belle sortie (visualisée dans ArcGIS avec 1: 1 pixel et une plage min / max constante):

DTM gaussien

Remarque: cette réponse a été mise à jour pour utiliser une fonction de traitement signal.fftconvolve beaucoup plus rapide basée sur la FFT .

Mike T
la source
1
+1 belle solution! Je ne sais pas avec certitude, mais il y a fort à parier que signal.convolve utilise des FFT.
whuber
Je cherchais un code flou pour un outil de couture automatique que j'écrivais et suis tombé par hasard sur cela. Beau travail @MikeToews!
Ragi Yaser Burhum
@RagiYaserBurhum J'aimerais en savoir plus sur votre outil. MikeToews Excellente réponse et extrait de code très apprécié.
Jay Laura
@JayLaura Rien de spécial, juste d'écrire un outil pour coudre automatiquement des images que j'ai prises avec des amis avec un ballon. Utilisation des classes d'Orfeo Toolbox orfeo-toolbox.org/SoftwareGuide/…
Ragi Yaser Burhum
2
@whuber lors de la révision de cette routine, il n'utilisait pas la FFT, mais c'est maintenant et c'est beaucoup plus rapide.
Mike T
4

Cela pourrait être un commentaire sur l'excellente réponse de MikeT , si ce n'était pas trop long et trop complexe. J'ai beaucoup joué avec elle et créé un plugin QGIS nommé Filtres de convolution FFT (au stade "expérimental") basé sur sa fonction. Outre le lissage, le plug-in peut également accentuer les contours en soustrayant le raster lissé de celui d'origine.

J'ai mis à niveau la fonction de Mike un peu plus tard:

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)

Les vérifications de validité vont de soi, mais ce qui est important, c’est de lancer en arrière. Auparavant, la fonction rendait les tableaux de nombres entiers noirs (zéros uniquement), en raison de la division par la somme des valeurs ( g / g.sum()).

Pavel V.
la source
3

Sous QGIS, j'ai facilement obtenu de bons résultats en utilisant le filtrage d'images Orfeo Toolbox . C'est assez rapide et le mode batch fonctionne bien. Des diffusions gaussiennes, moyennes ou anisotropes sont disponibles.

Notez que se Radiusréfère au nombre de cellules, pas à la distance.

Voici un exemple utilisant Smoothing (gaussian) :

  • Brut:

    Pas de filtre

  • Filtré:

    filtre

Tactopoda
la source
1

Belle solution pour le flou gaussien et une animation cool. En ce qui concerne l'outil de filtrage Esri mentionné ci-dessus, il s'agit essentiellement de l'outil Esri "Statistiques focales" codé en dur à une taille de 3x3. L'outil Statistiques focales vous donne beaucoup plus d'options sur la forme de votre filtre en mouvement, la taille et les statistiques que vous souhaitez exécuter. http://desktop.arcgis.com/fr/arcmap/latest/tools/spatial-analyst-toolbox/focal-statistics.htm

Vous pouvez également créer un filtre "irrégulier" dans lequel vous transmettez votre propre fichier texte avec des poids à utiliser pour chaque cellule. Le fichier texte contient autant de lignes que vous le souhaitez dans votre zone de filtre, avec des valeurs délimitées par des espaces pour les colonnes. J'imagine que vous devriez toujours utiliser un nombre impair de lignes et de colonnes afin que votre cellule cible se trouve au milieu.

J'ai créé un tableur Excel pour jouer avec différents poids que je viens de copier / coller dans ce fichier. Vous devriez obtenir les mêmes résultats que ci-dessus si vous ajustez les formules.

David A
la source