Модуль r.series.filter ГИС GRASS
| Создано в | 
|
Веб ГИС для вашей организации по доступной цене |
Назначение модуля
Модуль r.series.filter -- это модуль для геоинформационной системы GRASS, предназначенный для обработки временных рядов растровых данных, в первую очередь вегетационных индексов. Модуль производит очистку от шумов временных рядов с использованием различных фильтров (медианый фильтр, фильтр Савицкого-Голея). Модуль может использоваться для очистки от шумов временных рядов произвольной природы, однако одной из важных особенностей модуля является поддерка фильтрации данных вегетационных индексов. Для обработки таких рядов модуль реализует процедуру фильтрации, предложенную в статье "Chen J. et al. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter //Remote sensing of Environment. – 2004. – Т. 91. – №. 3. – С. 332-344.". Основная идея этого фильтра опирается на то, что шумы в значениях вегетационных индексов обычно индуцируются условиями съемки и атмосферными явлениями, которые занижают зачения индекса. В очень упрощенном виде можно сказать, что предлагаемая в указанной статье методика производит отбраковку низких значений индекса и производит фильтрацию по верхним значениям анализируемого временного ряда.
Установка
Модуль был написан на языке Python и предназначен для работы в GRASS GIS седьмой версии. В модуле используются библиотека scipy.signal, входящая в состав пакета SciPy. Поэтому перед началом работы в интерпретатор Python, используемый в GRASS, следует установить SciPy.
После установки SciPy появится возможность установить сам модуль r.series.filter. Данный модуль доступен в репозитории модулей GRASS ADDONS, поэтому его установка производится автоматически, если использовать модуль g.extension (вызываемый из графического интерфейса пользователя или командной строки). Например, можно использовать следущюую команду:
g.extension r.series.filter
В результате модуль будет автоматически скачан из репозитория и установлен.
Особенности работы и примеры
Параметры модуля и их назначение
Модуль находится в состоянии разработки, это означает, что в него будут добавляться новые возможности. На момент написания статьи модуль использует следующие параметры и ключи (показаны наиболее важные):
r.series.filter [-cu] input=string[,string,...] result_prefix=string
[method=string] [winsize=value] [order=value] [opt_points=value]
[diff_penalty=value] [deriv_penalty=value] [iterations=value]
Flags:
-c Try to find optimal parameters for filtering
-u Fit the result curve by upper boundary
Parameters:
input Raster names of equally spaced time series.
result_prefix Prefix for raster names of filtered X(t)
method Used method
default: savgol
winsize Length of running window for the filter
default: 9
order Order of the Savitzky-Golay filter
default: 2
opt_points Count of random points used for parameter optimization
default: 50
diff_penalty Penalty for difference between original and filtered signals
default: 1.0
deriv_penalty Penalty for big derivates of the filtered signal
default: 1.0
iterations Number of iterations
default: 1
Параметры общего назначения
Два параметра input и result_prefix отвечают за перечень входных и выходных растров. Модуль ожидает, что на вход ему будет передан список растров, перечисленных через запятую в порядке возрастания (или убывания) времени. Например:
input=ndvi1,ndvi2,ndvi3,ndvi4,....
При этом временной ряд должен быть равноотстоящим, т.е. растры представляют собой измерения какой-либо величины (или состояния территории), сделанные через равные промежутки времени. Но на практике не всегда возможно иметь дело с равноотстоящими рядами, часть каких-либо измерений может отсутсвовать по объетивным причинам. В таком случае следует вставить пустой растр, состоящий из значений NULL (нет данных), который легко сгенерировать, воспользовавшись растровым калькулятором:
r.mapcalc expression="dummy = null()"
В самих растрах так же могут быть области, с отсутсвующими в них данными, такие участки не требуют специальной предварительной обработки.
Перед началом расчетов такие пропущенные значения будут предварительно заполнены данными, полученными линейной интерполяцией (во времени, пространственная интерполяция не используется).
Результаты будут записаны в новые растры, названия которых начинаются с указанного result_prefix. Например, в случае, если result_prefix="res.", будут получены растры: res.ndvi1,res.ndvi2,res.ndvi3,res.ndvi4,....
Таким образом команда со значениями, заданными по умолчанию может выглядеть следующим образом:
r.series.filter input=ndvi1,ndvi2,ndvi3,ndvi4,.... result_prefix="res."
Методы фильтрации
На момент написания статьи в модуле реализовано два метода фильтрации:
- медианный фильтр (параметр method=median);
- фильтр Савицкого-Голея (method=savgol).
Каждый из этих фильтров работает в скользящем окне и рассчитывает прогнозируемое значение исследуемой величины для центральной точки окна. Поэтому модуль использует еще один параметр winsize, задающий ширину скользящего окна. Этот параметр должен быть нечетным числом.
Кроме того у метода Савцикого-Голея есть еще один параметр -- порядок полинома (order), который используется при аппроксимации временного ряда внутри скользящего окна.
Естественно, что ширина окна и порядок полинома должны быть таковы, что позволят рассчитать коэффициенты полинома. В случае, если порядок полиома превышает допустимый для данной величины окна, будет выведено сообщение об ошибке.
Например, для использования фильтра Савицкого-Голея с порядком полинома, равным трем, и шириной окна 17 отсчетов, нужно использовать следующую команду:
r.series.filter input=.... result_prefix=... method=savgol order=3 wunsiz=17
Иногда есть смысл применить последовательно сразу несколько итераций фильтрации (особенно это актуально для медианного фильтра): сначала произвести фильтрацию исходных данных, затем применить этот же фильтр к результам фильтрации и т.д. Поэтому пользователь может уаказать параметр iterations, отвечающий за число итераций, например:
r.series.filter input=.... result_prefix=... method=median iterations=5 winsize=21
Подбор оптимальных значений параметров
Часто бывает сложно понять, какие значения параметров являются "хорошими" для конкретного набора растров.
Понятно, что эти величины зависят как и от форм анализируемых кривых, так и от целей исследователя, но в довольно обшем случае можно сказать, что качество подгонки может выражаться в виде некоторой функции от точности и "гладкости" результирующей кривой. Точность рассчитать легко как средний модуль разности между исходной кривой и результирующей, а гладкость можно оценить средним модулем разности между соседними значениями результирующей кривой. К сожалению, обычно бывает так, что гладкая и "красивая" кривая является слишком грубой и неточной, и наоборот -- точная кривая обычно повторяет все флуктуации в данных, в том числе и случайные шумы. Поэтому пользователь должен решить, что ему более важно, точность или гладкость.
На этом строится процедура подбора параметров фильтрации: пользователь задает два числа -- штрафы за неточность и негладкость результирующей кривой (параметры модуля diff_penalty и deriv_penalty). Эти штрафы показывают, насколько пользователю важна точность в ущерб гладкости. Затем запускается процедура подбора параметров, работающая следующим образом:
- Выбираются начальные значения параметра ширины окна для медианного фильра (ширины окна и порядка полинома для фильтра Савицкого-Голея).
- Выбирается opt_points случайных пикселей.
- Для выбранных пикселей производится фильтрация с заданными параметрами.
- Для каждого выбранного пикселя вычилсяется E - средний модуль разности между исходной кривой и результирующей, а также D - средний модуль разности между соседними значениями результирующей кривой.
- Для каждого пикселя находится взвешенная сумма штрафов: diff_penalty * E + deriv_penalty * D, а затем расчитывается суммарный штраф по всем пикселям. Этот суммарный штраф является оценкой качества подгонки кривых для заданных параметров фильтра.
Затем запускается процедура поиска таких параметров, при которых суммарный штраф будет минимальным. После нахождения таких параметров производится фильтрация всего временного растрового ряда.
Процедура подбора параметров включается флагом "-c", например:
r.series.filter -с input=.... result_prefix=... method=median
Фильтрация вегетационных индексов
Как отмечалось выше, в статье "Chen J. et al. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter //Remote sensing of Environment. – 2004. – Т. 91. – №. 3. – С. 332-344." предлагается метод фильтрации временных рядов вегетационных индексов, исходя из предположения, что атмосферные явления обычно занижают значения индексов. В исходной статье медод предлагался для фильра Савицкого-Голея, но с таким же успехом его можно применять и к другим фильтрам. Поэтому в модуле r.series.filter пользователь может включить флагом "-u" данную модификацию фильтра:
r.series.filter -u input=.... result_prefix=... method=median winsize=15
Предлагаемый в статье метод является итерационным, он действует на постепенно уменьшающихся величинах ширины скользящего окна, поэтому параметр winsize в данном случае будет означать начальную ширину окна.
Более подробно об функционировании фильтра см. в указанной выше статье.