Какие инструменты растрового сглаживания / обобщения доступны?

46

У меня есть ЦМР, который я хотел бы сгладить или обобщить, чтобы убрать топографические крайности (отрубить пики и заполнить долины). В идеале мне также хотелось бы иметь контроль над радиусом или уровнем «размытости». В конце мне понадобится набор растров от слегка размытых до действительно размытых. (Теоретически, размытым будет постоянный растр среднего арифметического всех значений).

Могу ли я использовать какие-либо инструменты или методы (основанные на Esri, GDAL, GRASS)? Нужно ли мне дома испечь мою собственную процедуру размытия по Гауссу ? Могу ли я использовать фильтр нижних частот (например, фильтр ArcGIS ), и если да, нужно ли запускать его несколько раз, чтобы получить эффект большого радиуса?

Майк Т
источник
Как насчет простого экспорта растра в ячейку большего размера? Не приведет ли это также к приглушению крайностей?
1
Да, это также уменьшило бы экстремумы (при условии, что неявная повторная выборка включает некоторую форму усреднения), но это ужасный способ сглаживания матрицы высот: вы бы создали небольшое количество больших блоков. Кстати, для этого обычно не нужно экспортировать растр; агрегации , а также передискретизации к другому cellsize основные операции , обычно находятся в растровом на основе программного обеспечения.
whuber

Ответы:

29

Размытие по Гауссу - это только средневзвешенное фокусное значение. Вы можете воссоздать его с высокой точностью с помощью последовательности коротких расстояний по окружности (невзвешенных): это приложение центральной предельной теоремы .

У вас есть много вариантов. «Фильтр» слишком ограничен - он только для 3х3 окрестностей - так что не беспокойтесь об этом. Наилучший вариант для больших ЦМР - перенести вычисления за пределы ArcGIS в среду, в которой используются быстрые преобразования Фурье: они выполняют те же фокусные вычисления, но (для сравнения) они делают это невероятно быстро. (У GRASS есть модуль FFT . Он предназначен для обработки изображений, но вы можете использовать его для своей ЦМР, если сможете с разумной точностью изменить его масштаб до диапазона 0..255.) За исключением этого, как минимум два решения : стоит учесть:

  1. Создайте набор весов окрестностей для аппроксимации размытия по Гауссу для значительного соседства. Используйте последовательные проходы этого размытия, чтобы создать свою последовательность более плавных ЦМР.

    (Веса вычисляются как exp (-d ^ 2 / (2r)), где d - расстояние (в ячейках, если хотите), а r - эффективный радиус (также в ячейках). Они должны быть вычислены в пределах круга, продолжающегося по крайней мере до 3r . После этого разделите каждый вес на сумму их всех так, чтобы в конце они составили 1)

  2. Кроме того, забудьте о взвешивании; просто запустите круговое фокусное среднее значение несколько раз. Я сделал именно это для изучения того, как производные сетки (такие как уклон и аспект) изменяются с разрешением матрицы высот.

Оба метода будут работать хорошо, и после первых нескольких проходов будет мало выбора между этими двумя, но есть и убывающая отдача: эффективный радиус из n последовательных средних значений фокуса (все используют одинаковый размер окрестности) только (приблизительно) квадратный корень из n раз радиуса среднего очага. Таким образом, для огромного количества размытия вы захотите начать сначала с окрестности большого радиуса. Если вы используете невзвешенное среднее значение, выполните 5-6 проходов по ЦМР. Если вы используете веса, приблизительно равные гауссову, вам потребуется только один проход, но вы должны создать матрицу весов.

Этот подход действительно имеет среднее арифметическое DEM в качестве предельного значения.

Whuber
источник
1
Если у ваших данных есть всплески, вы можете попробовать медианный фильтр ( en.wikipedia.org/wiki/Median_filter ), прежде чем применять более общее размытие, как предлагает whuber.
MerseyViking
@Mersey Это отличное предложение. Я никогда не видел ЦМР с локальными выбросами, но опять же мне никогда не приходилось обрабатывать необработанные ЦМР (такие как необработанные результаты LIDAR). Вы не можете делать медианные фильтры с БПФ, но вам (обычно) нужно только 3 х 3 соседства, так что в любом случае это быстрая операция.
whuber
Спасибо, что. Я должен признать, что когда-либо использовал только предварительно обработанные данные LiDAR, но есть некоторые существенные всплески в данных SRTM, которые выиграли бы от медианного фильтра. Они имеют тенденцию иметь ширину 2 или 3 образца, поэтому потребуется больший медианный фильтр.
MerseyViking
@Mersey Вы все еще в порядке с большим медианным фильтром 5 x 5 или 7 x 7. Однако, если вы рассматриваете (скажем) фильтр 101 x 101, будьте готовы ждать! Вы также предлагаете уточнить важный момент: очень полезно выполнить предварительный анализ матрицы высот, прежде чем что-либо делать. Это будет включать в себя определение шипов (локальные выбросы) и характеристику их размеров и размеров. Вы должны быть уверены, что это действительно артефакты (а не какое-то реальное явление), прежде чем вы начнете стирать их с помощью фильтра!
whuber
1
+1 для БПФ на данных высоты. Я фактически сделал эту работу в траве для 32-битных данных NED, чтобы удалить двунаправленное чередование. В конце концов, это также было проблематично, поскольку в нем вновь появился эффект террасирования, который поражает многие другие ЦМР, полученные из контуров.
Джей Гварнери
43

Я изучаю подход SciPy signal.convolve (основанный на этой кулинарной книге ), и у меня действительно хороший успех со следующим фрагментом:

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

Я использую это в другой функции, которая читает / записывает файлы GeoTIFF float32 через GDAL (нет необходимости изменять масштаб до 0-255 байт для обработки изображений), и я использую попытку размера пикселя (например, 2, 5, 20), и это имеет действительно хороший вывод (визуализируется в ArcGIS с 1: 1 пикселем и постоянным минимальным / максимальным диапазоном):

Gaussian DTM

Примечание: этот ответ был обновлен для использования гораздо более быстрой функции обработки signal.fftconvolve на основе FFT .

Майк Т
источник
1
+1 Отличное решение! Я точно не знаю, но стоит поспорить, что signal.convolve использует БПФ.
whuber
Я искал некоторый размытый код для инструмента автоматического сшивания, который я пишу, и наткнулся на это. Хорошая работа @MikeToews!
Раги Язер Бурхум
@RagiYaserBurhum Хотелось бы узнать больше о вашем инструменте. MikeToews Отличный ответ и высоко ценимый фрагмент кода.
Джей Лора
@JayLaura Ничего особенного, просто написание инструмента для автоматического сшивания изображений, которые я сделал с друзьями с помощью воздушного шара. Использование классов Orfeo Toolbox orfeo-toolbox.org/SoftwareGuide/…
Раги Язер Бурхум
2
@whuber после пересмотра этой процедуры, он не использовал FFT, но сейчас и намного быстрее.
Майк Т,
4

Это может быть комментарием к отличному ответу MikeT , если он не был слишком длинным и слишком сложным. Я много играл с ним и создал плагин QGIS под названием FFT Convolution Filters (пока что «на экспериментальной» стадии), основанный на его функции. Помимо сглаживания, плагин также может делать резкие края, вычитая сглаженный растр из исходного.

Я немного улучшил функцию Майка в процессе:

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)

Проверки достоверности довольно очевидны, но важно то, что они приводятся на плавание и обратно. До этого функция делала целочисленные массивы черными (только нули) из-за деления на сумму значений ( g / g.sum()).

Павел В.
источник
3

В QGIS я легко добился хороших результатов, используя фильтрацию изображений Orfeo Toolbox . Это разумно быстро и пакетный режим работает нормально. Доступны гауссовы, средние или анизотропные диффузии.

Обратите внимание, что Radiusотносится к числу ячеек, а не расстояние.

Вот пример использования Smoothing (гауссовский) :

  • сырье:

    Нет фильтра

  • Отфильтрованный:

    фильтр

Tactopoda
источник
1

Хорошее решение для размытия по Гауссу и классной анимации. Что касается упомянутого выше инструмента «Фильтр Esri», то это просто инструмент «Фокальная статистика» Esri, жестко запрограммированный до размера 3x3. Инструмент Focal Statistics дает вам гораздо больше возможностей по форме вашего движущегося фильтра, размеру и статистике, которую вы хотите запустить. http://desktop.arcgis.com/en/arcmap/latest/tools/spatial-analyst-toolbox/focal-statistics.htm

Вы также можете создать «нерегулярный» фильтр, в котором вы передадите свой собственный текстовый файл с весами для использования в каждой ячейке. Текстовый файл имеет столько строк, сколько вы хотите в области фильтра, с разделенными пробелами значениями для столбцов. Я думаю, вы всегда должны использовать нечетное количество строк и столбцов, поэтому ваша целевая ячейка находится посередине.

Я создал таблицу Excel для игры с разными весами, которую просто копирую / вставляю в этот файл. Это должно привести к тем же результатам, что и выше, если вы корректируете формулы.

Дэвид А
источник