Классический анализ данных ДЗЗ происходит по на основе обработки яркостей пикселей в нескольких спектральных каналах. Другой подход связан с анализом динамики яркостей в течении длительгого периода. Суть этих методов основана на том, что различные классы Земной поверхности (здесь, в первую очередь, речь идет о растительности) обладают различными фенологическими характеристиками. Эти характеристики порождают различия во временных рядах NDVI, поэтому вместо анализа яркостей пикселей на отдельных снимках можно анализировать динамику NDVI.
Пример такого подхода приводится в Lhermitte S. et al. Hierarchical image segmentation based on similarity of NDVI time series //Remote Sensing of Environment. – 2008. – Т. 112. – №. 2. – С. 506-521. В статье предлагается найти разложение в ряд Фурье для временного ряда NDVI. Коэффициенты ряда Фурье однозначно описывают исходный ряд и поэтому классификацию Земной поверхности можно можно произвести на базе этих коэффициентов.
Ниже производится эта процедура в очень упрощенном виде, поскольку в цитируемой работе анализируется не только временные, но и пространственные вариации NDVI. Здесь в качестве данных для анализа используется продукт MOD13 (16-ти дневные композиты).
Эта статья строится следующим образом:
Дана функция $f(x)$, определенная на интервале $x\in (-\pi, pi]$. Известно, что ее можно разложить в ряд, состоящий из тригонометрических функций (ограничения на $f(x)$ очень слабые, на практике любая функция будет им удовлетворять):
$$ f(x) = \frac{a_0}2 + \sum_{n=0}^\infty (a_n\cos nx + b_n \sin nx) $$Легко показать, что коэффициенты ряда могут быть вычисленны по формулам:
$$ a_n = \frac1{\pi} \int_{-\pi}^\pi f(x) \cos nx dx, $$для $n \ge 0$ и $$ b_n = \frac1{\pi} \int_{-\pi}^\pi f(x) \sin nx dx $$
для $n >0$.
Исходный интервал $x\in (-\pi, \pi]$ может быть легко промасштабирован на произвольный симметричный отрезов $[-L, L]$. Кроме того, формула разложения часто записывается в более простом виде в комплексной форме: $$ f(x) = \sum_{-\infty}^{\infty} c_n e^{in\pi x / L}, x\in [-L, L]. $$ (Последнее выражение базируется на формуле Эйлера: $e^{i\phi} = \cos \phi + i \sin \phi$. Нужно обратить внимание на то, что хотя теперь ряд Фурье стал комплексным, сама функция $f(x)$ по прежнему остается действительной).
Преобразование Фурье -- это преобразование функции, заданное на всей числовой оси $x\in (-\pi, \pi)$ следующим образом:
$$ F(k) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-i kx} f(x) dx. $$Обратное преобразование (переход от $F(k)$ к $f(x)$) задается формулой:
$$ f(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-i kx} F(k) dk. $$Таким образом преобразование Фурье возвращает результат, представленный в виде комплексных чисел. При анализе сигналов при помощи преобразования Фурье важную роль играют понятия амплитуды (модуля) и частоты (аргумента) комплексного числа.
Для вычисления коэффициентов преобразования Фурье существует классический алгоритм, который носит название "Быстрое преборазование Фурье". Вот его основные особенности:
In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft, fftshift, fftfreq
from tools import get_modis_ndvi_frame
In [2]:
water_files = ['data/MODIS/point_1_NDVI.csv', 'data/MODIS/point_2_NDVI.csv', 'data/MODIS/point_3_NDVI.csv']
water_frames = map(get_modis_ndvi_frame, water_files)
plot(water_frames[0].DayNum, water_frames[0].NDVI,
water_frames[1].DayNum, water_frames[1].NDVI,
water_frames[2].DayNum, water_frames[2].NDVI)
Out[2]:
Откровенно говоря я не ожидал, что вода будет так сильно меняться в годичном цикле. (Возможно, я перепутал названия классов, т.к. данные по NDVI были получены из продуктов MODIS, а названия классов при интерпретации Landsat. Возможно, дело в том, что в ячейка растра MODIS хранит "интегральный" сигнал нескольких "реальных" классов. Может быть, вода "цветет". Короче, много чего может быть, но для нижеследующих примеров это не так важно -- важно показать принцип работы.)
Найдем разложение в ряд Фурье для первой точки и выведем первые 5 коэффициентов разложения:
In [3]:
water0 = water_frames[0].NDVI
water0_fft = fft(water0)
water0_fft[:5]
Out[3]:
Как видим, коэффициенты являются комплексными числами. Произведем обратное преобразование Фурье, убедимся, что на выходе мы получаем первоначальный ряд:
In [4]:
print np.array(water0[:5]) # Исходный сигнал
print ifft(fft(water0))[:5] # Преобразование Фурье (туда-сюда)
print ifft(water0_fft).real[:5] # Действительная часть сигнала
Здесь видно, что на самом деле после обратного преобразования Фурье получаются комплексные числа, но мнимая часть их очень близка к нулю (появление ненулевой комплексной части объясняется погрешностями вычислений). Поэтому эти мнимые части можно отбросить:
In [5]:
plot(water0)
Out[5]:
In [6]:
water0_fft = fft(water0)
water0_restored = ifft(water0_fft).real
plot(water0_restored)
Out[6]:
Функция fft возвращает "сдвинутые" вдоль числовой оси коэффициенты (подробнее об этом см. http://docs.scipy.org/doc/numpy/reference/routines.fft.html#background-information): в элементе water0_fft[0] хранится константа (среднее значение сигнала), water0_fft[1:n/2] коэффициенты, соответствующие положительным частотам, а water0_fft[n/2+1:] -- коэффициенты, соответствующие отрицательным частотам. Для того, чтобы получить коэффициенты, следующие в "естественном" порядке, можно воспользоваться функцией fftshift, которая переупорядочивает эти коэффициенты. Альтернативный путь -- воспользоваться функцией fftfreq, которая, наборот, подстраивает массив частот под порядок коэффициентов.
Пусть a -- исходный сигнал, зависящий от времени, A = fft(A), тогда np.abs(A) возвращает амплитудный спектр сигнала, np.abs(A)**2 -- энергетический спект сигнала, а np.angle(A) -- фазовый спектр.
(Хорошее практическое введение в тему есть тут: http://habrahabr.ru/post/112068/, правда оно рассматривает MATLAB в качестве языка реализации).
Построим амплитудный и энергетический спектры:
In [7]:
# Амплитудный спектр
freq = fftfreq(water0.shape[0])
water0_fft = fft(water0)
idx = argsort(freq) # индексы, чтобы частоты на графике шли по возрастанию
plot(freq[idx], abs(water0_fft)[idx], '-o')
Out[7]:
In [8]:
# Энергетический спектр
plot(freq[idx], abs(water0_fft[idx])**2, '-o')
Out[8]:
После построения энергетического спектра становится очевидным, что почити вся энергия сигнала содержится в паре-тройке компонент (гармоник). В данном случае гармоники с высокой энергией -- это $A_0$ (константа), а также $A_{3}$ и $A_{-3}$ (в силу симметрии эти два коэффициента -- одна и та же гармоника). Эти гармоники соответствуют выраженным частотам, так, например, гармоники с индексом "3" ($A_{3}$ и $A_{-3}$) сигнализируют о том, что в наших данных присутствуют три хорошо выраженных "волны" (цикла).
Остальные гармоники несут очень незначительную энергию и могут быть опущены. Поэтому "занулим" все коэффициенты разложения, кроме $A_0$, $A_{3}$ и $A_{-3}$, и восстановим сигнал из полученной последовательности. (Для простоты реализации обнулим все коэффициенты, модуль которых меньше определенного порога).
In [9]:
# Зануляем
water0_fft = fft(water0)
water0_fft[abs(water0_fft)< 50000] = 0
# Востанавливаем
water0_filt = ifft(water0_fft).real
plot(water0_filt)
Out[9]:
На самом деле, учитывая, что мы оставили только три коэффициента ($A_0$, $A_3$, $A_{-3}$, причем, в силу симметрии, по сути имеем только два коэффициента), мы получили на выходе синусоиду, которая приближает наш сигнал: $A_0 + A_3\sin(...)$. Покажем для наглядности на одном графике исходный и восстановленный сигнал:
In [10]:
t = range(len(water0))
plot(t, water0, t, water0_filt)
Out[10]:
Если оставить большее количество коэффициентов, получим иное приближение:
In [11]:
# Зануляем
water0_fft = fft(water0)
water0_fft[abs(water0_fft)< 20000] = 0
# Востанавливаем
water0_filt = ifft(water0_fft).real
plot(t, water0, t, water0_filt)
Out[11]:
Как было сказано выше, идея, лежащая в основе, -- проста: основную форму кривой описывают коэффициенты разложения, которые несут наибольшую энергию, остальные коэффициенты описывают "второстепенные" компоненты сигнала. Поэтому можно выделить несколько основных коэффициентов и использовать их в качестве признаков при классификации.
Проверим это на практике:
In [12]:
water_files = ['data/MODIS/point_1_NDVI.csv', 'data/MODIS/point_2_NDVI.csv', 'data/MODIS/point_3_NDVI.csv']
water_frames = map(get_modis_ndvi_frame, water_files)
plot(water_frames[0].DayNum, water_frames[0].NDVI,
water_frames[1].DayNum, water_frames[1].NDVI,
water_frames[2].DayNum, water_frames[2].NDVI)
Out[12]:
In [13]:
pine_swamp_files = ['data/MODIS/point_4_NDVI.csv', 'data/MODIS/point_5_NDVI.csv', 'data/MODIS/point_6_NDVI.csv']
pine_swamp_frames = map(get_modis_ndvi_frame, pine_swamp_files)
plot(pine_swamp_frames[0].DayNum, pine_swamp_frames[0].NDVI,
pine_swamp_frames[1].DayNum, pine_swamp_frames[1].NDVI,
pine_swamp_frames[2].DayNum, pine_swamp_frames[2].NDVI)
Out[13]:
In [14]:
rostrata_files = ['data/MODIS/point_7_NDVI.csv', 'data/MODIS/point_8_NDVI.csv', 'data/MODIS/point_9_NDVI.csv']
rostrata_frames = map(get_modis_ndvi_frame, rostrata_files)
plot(rostrata_frames[0].DayNum, rostrata_frames[0].NDVI,
rostrata_frames[1].DayNum, rostrata_frames[1].NDVI,
rostrata_frames[2].DayNum, rostrata_frames[2].NDVI)
Out[14]:
In [15]:
pine_files = ['data/MODIS/point_10_NDVI.csv', 'data/MODIS/point_11_NDVI.csv', 'data/MODIS/point_12_NDVI.csv']
pine_frames = map(get_modis_ndvi_frame, pine_files)
plot(pine_frames[0].DayNum, pine_frames[0].NDVI,
pine_frames[1].DayNum, pine_frames[1].NDVI,
pine_frames[2].DayNum, pine_frames[2].NDVI)
Out[15]:
In [16]:
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(water_frames[0].DayNum, water_frames[0].NDVI, label='Water')
plt.plot(pine_swamp_frames[0].DayNum, pine_swamp_frames[0].NDVI, label='Pine+swamp')
plt.plot(rostrata_frames[0].DayNum, rostrata_frames[0].NDVI, label='Rostrata')
plt.plot(pine_frames[0].DayNum, pine_frames[0].NDVI, label='Pine')
plt.legend(loc=4)
plt.show()
Как видим, единственный объект, который хорошо отличается от остальных -- вода. Водные объекты имеют более низкий NDVI и более короткий "вегетационный период", если так вообще можно говорить о водных объектах :). У всех остальных объектов вегетация начинается и заканчивается приблизительно в одно время.
У сосняков относительно высокий NDVI в зимний период.
Отличается ли один тип болота от другого, сказать сложно: график осокового болота лежит чуть ниже, чем график сосново-кустарничкого болота, но насколько это различие значимо статистически по одной точке не скажешь -- нужно привлечь больше статистики.
Поместим все имеющиеся точки болот на один график:
In [17]:
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(pine_swamp_frames[0].DayNum, pine_swamp_frames[0].NDVI, 'r', label='Pine+swamp')
plt.plot(pine_swamp_frames[1].DayNum, pine_swamp_frames[1].NDVI, 'r')
plt.plot(pine_swamp_frames[2].DayNum, pine_swamp_frames[2].NDVI, 'r')
plt.plot(rostrata_frames[0].DayNum, rostrata_frames[0].NDVI, 'b', label='Rostrata')
plt.plot(rostrata_frames[1].DayNum, rostrata_frames[1].NDVI, 'b')
plt.plot(rostrata_frames[2].DayNum, rostrata_frames[2].NDVI, 'b')
plt.legend(loc=1)
plt.show()
Видимо, основное отличие их состоит в том, что в зимний период сосновое болото имеет более высокий NDVI.
Для работы над разными рядами нам обязательно нужно, чтобы данные были за один и тот же период времени. Убедимся в этом:
In [18]:
for num in [0, 1, 2]:
for frames in [water_frames, pine_swamp_frames,
rostrata_frames, pine_frames]:
print np.all(water_frames[0].DayNum == frames[1].DayNum),
Выполним преобразование Фурье
In [19]:
water0 = water_frames[0].NDVI
water1 = water_frames[1].NDVI
water2 = water_frames[2].NDVI
water0_fft = fft(water0)
water1_fft = fft(water1)
water2_fft = fft(water2)
Построим графики амплитудных спектров
In [20]:
freq = fftfreq(water0.shape[0])
idx = argsort(freq) # индексы, чтобы частоты на графике шли по возрастанию
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(freq[idx], abs(water0_fft[idx]), '-o', label='0')
plt.plot(freq[idx], abs(water1_fft[idx]), '-o', label='1')
plt.plot(freq[idx], abs(water2_fft[idx]), '-o', label='2')
plt.legend(loc=1)
plt.show()
Подберем порог обнуления коэффициентов разложения, при котором сохряняются основные особенности сигнала и запомним оставшиеся коэффициенты:
In [21]:
water0_fft = fft(water0)
water1_fft = fft(water1)
water2_fft = fft(water2)
tol = 20000
water0_fft[abs(water0_fft)< tol] = 0
water1_fft[abs(water1_fft)< tol] = 0
water2_fft[abs(water2_fft)< tol] = 0
water0_filt = ifft(water0_fft).real
water1_filt = ifft(water1_fft).real
water2_filt = ifft(water2_fft).real
# Графики
fig = plt.figure(figsize=(12,6), dpi=100)
# plt.plot(water0, 'b-', label='0')
# plt.plot(water1, 'g-', label='1')
# plt.plot(water2, 'r-', label='2')
plt.plot(water0_filt, 'b-o', label='0 filt')
plt.plot(water1_filt, 'g-o', label='1 filt')
plt.plot(water2_filt, 'r-o', label='2 filt')
plt.legend(loc=1)
plt.show()
Частоты, на которых остались ненулевые коэффициенты:
In [22]:
print freq[water0_fft != 0]
print freq[water0_fft != 0]
print freq[water0_fft != 0]
print freq[0], freq[3], freq[9]
С учетом симметриии, получаем, что форму этих кривых определяют в основном три коэффициента: $A_0$, $A_3$, $A_9$
Повторим процедуру, которую мы произвели с водой, на остальных данных
In [23]:
pine_swamp0 = pine_swamp_frames[0].NDVI
pine_swamp1 = pine_swamp_frames[1].NDVI
pine_swamp2 = pine_swamp_frames[2].NDVI
pine_swamp0_fft = fft(pine_swamp0)
pine_swamp1_fft = fft(pine_swamp1)
pine_swamp2_fft = fft(pine_swamp2)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(freq[idx], abs(pine_swamp0_fft[idx]), '-o', label='0')
plt.plot(freq[idx], abs(pine_swamp1_fft[idx]), '-o', label='1')
plt.plot(freq[idx], abs(pine_swamp2_fft[idx]), '-o', label='2')
plt.legend(loc=1)
plt.show()
In [24]:
pine_swamp0_fft = fft(pine_swamp0)
pine_swamp1_fft = fft(pine_swamp1)
pine_swamp2_fft = fft(pine_swamp2)
tol = 20000
pine_swamp0_fft[abs(pine_swamp0_fft)< tol] = 0
pine_swamp1_fft[abs(pine_swamp1_fft)< tol] = 0
pine_swamp2_fft[abs(pine_swamp2_fft)< tol] = 0
pine_swamp0_filt = ifft(pine_swamp0_fft).real
pine_swamp1_filt = ifft(pine_swamp1_fft).real
pine_swamp2_filt = ifft(pine_swamp2_fft).real
# Графики
fig = plt.figure(figsize=(12,6), dpi=100)
# plt.plot(water0, 'b-', label='0')
# plt.plot(water1, 'g-', label='1')
# plt.plot(water2, 'r-', label='2')
plt.plot(pine_swamp0_filt, 'b-o', label='0 filt')
plt.plot(pine_swamp1_filt, 'g-o', label='1 filt')
plt.plot(pine_swamp2_filt, 'r-o', label='2 filt')
plt.legend(loc=1)
plt.show()
print freq[pine_swamp0_fft != 0]
print freq[pine_swamp1_fft != 0]
print freq[pine_swamp1_fft != 0]
print freq[0], freq[3], freq[9], freq[12]
Здесь фигурируют 4 незануленных коэффициента: $A_0$, $A_3$, $A_9$, $A_{12}$. (При пороге в 16000 появляется 5-й коэф. $A_6$, но он приводит к противоестественной "волне" на макушке, т.е. в середине лета)
In [25]:
rostrata0 = rostrata_frames[0].NDVI
rostrata1 = rostrata_frames[1].NDVI
rostrata2 = rostrata_frames[2].NDVI
rostrata0_fft = fft(rostrata0)
rostrata1_fft = fft(rostrata1)
rostrata2_fft = fft(rostrata2)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(freq[idx], abs(rostrata0_fft[idx]), '-o', label='0')
plt.plot(freq[idx], abs(rostrata1_fft[idx]), '-o', label='1')
plt.plot(freq[idx], abs(rostrata2_fft[idx]), '-o', label='2')
plt.legend(loc=1)
plt.show()
In [26]:
rostrata0_fft = fft(rostrata0)
rostrata1_fft = fft(rostrata1)
rostrata2_fft = fft(rostrata2)
tol = 20000
rostrata0_fft[abs(rostrata0_fft)< tol] = 0
rostrata1_fft[abs(rostrata1_fft)< tol] = 0
rostrata2_fft[abs(rostrata2_fft)< tol] = 0
rostrata0_filt = ifft(rostrata0_fft).real
rostrata1_filt = ifft(rostrata1_fft).real
rostrata2_filt = ifft(rostrata2_fft).real
# Графики
fig = plt.figure(figsize=(12,6), dpi=100)
# plt.plot(water0, 'b-', label='0')
# plt.plot(water1, 'g-', label='1')
# plt.plot(water2, 'r-', label='2')
plt.plot(rostrata0_filt, 'b-o', label='0 filt')
plt.plot(rostrata1_filt, 'g-o', label='1 filt')
plt.plot(rostrata2_filt, 'r-o', label='2 filt')
plt.legend(loc=1)
plt.show()
print freq[rostrata0_fft != 0]
print freq[rostrata1_fft != 0]
print freq[rostrata2_fft != 0]
print freq[0], freq[3], freq[9], freq[12]
In [27]:
pine0 = pine_frames[0].NDVI
pine1 = pine_frames[1].NDVI
pine2 = pine_frames[2].NDVI
pine0_fft = fft(pine0)
pine1_fft = fft(pine1)
pine2_fft = fft(pine2)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(freq[idx], abs(pine0_fft[idx]), '-o', label='0')
plt.plot(freq[idx], abs(pine1_fft[idx]), '-o', label='1')
plt.plot(freq[idx], abs(pine2_fft[idx]), '-o', label='2')
plt.legend(loc=1)
plt.show()
In [28]:
pine0_fft = fft(pine0)
pine1_fft = fft(pine1)
pine2_fft = fft(pine2)
tol = 20000
pine0_fft[abs(pine0_fft)< tol] = 0
pine1_fft[abs(pine1_fft)< tol] = 0
pine2_fft[abs(pine2_fft)< tol] = 0
pine0_filt = ifft(pine0_fft).real
pine1_filt = ifft(pine1_fft).real
pine2_filt = ifft(pine2_fft).real
# Графики
fig = plt.figure(figsize=(12,6), dpi=100)
# plt.plot(water0, 'b-', label='0')
# plt.plot(water1, 'g-', label='1')
# plt.plot(water2, 'r-', label='2')
plt.plot(pine0_filt, 'b-o', label='0 filt')
plt.plot(pine1_filt, 'g-o', label='1 filt')
plt.plot(pine2_filt, 'r-o', label='2 filt')
plt.legend(loc=1)
plt.show()
print freq[pine0_fft != 0]
print freq[pine1_fft != 0]
print freq[pine2_fft != 0]
print freq[0], freq[3], freq[9], freq[12]
Как видно из графиков, форму кривых более-менее точно описывают в основном коэффициэты $A_0$, $A_3$, $A_9$, $A_{12}$ (для некоторых кривых достаточно и меньшего количества коэффициетов, но мы выберем те, что описывают все кривые, а не только часть из них).
Построим на одном графике отфильтрованные сигналы с этими коэффициентами (на примере одной точки из каждого класса).
In [29]:
water0_fft = fft(water0)
pine_swamp0_fft = fft(pine_swamp0)
rostrata0_fft = fft(rostrata0)
pine0_fft = fft(pine0)
# Обнулим коэффициенты
size = water0_fft.shape
idx = [0, 3, 9, 12, -3, -9, -12] # Индексы коэффициентов, которые не зануляем
water0_fft_upd = np.zeros(size, dtype=complex)
pine_swamp0_fft_upd = np.zeros(size, dtype=complex)
rostrata0_fft_upd = np.zeros(size, dtype=complex)
pine0_fft_upd = np.zeros(size, dtype=complex)
for t in idx:
water0_fft_upd[t] = water0_fft[t]
pine_swamp0_fft_upd[t] = pine_swamp0_fft[t]
rostrata0_fft_upd[t] = rostrata0_fft[t]
pine0_fft_upd[t] = pine0_fft[t]
# Обратное преобразование
water0_filt = ifft(water0_fft_upd).real
pine_swamp0_filt = ifft(pine_swamp0_fft_upd).real
rostrata0_filt = ifft(rostrata0_fft_upd).real
pine0_filt = ifft(pine0_fft_upd).real
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(water0_filt, '-o', label='water')
plt.plot(pine_swamp0_filt, label='pine_swamp')
plt.plot(rostrata0_filt, '-o', label='rostrata')
plt.plot(pine0_filt, label='pine')
plt.legend(loc=3)
plt.show()
Как видим, на этих графиках прослеживаются основные тенденции, которые были подмечены на "сырых" NDVI:
Водные объекты имеют более низкий NDVI и более короткий "вегетационный период". У всех остальных объектов вегетация начинается и заканчивается приблизительно в одно время.
У сосняков относительно высокий NDVI в зимний период.
Сосновые и осоковые болота очень похожи друг на друга. Наиболее характерное отличие -- в зимний период сосновое болото имеет более высокий NDVI, чем осоковое.
Выше мы выделили четыре коэффициента, которые можно использовать для классификации $A_0$, $A_3$, $A_9$, $A_{12}$. Но в данном случае у нас имеется всего двенадцать различных точек, поэтому нет смысла использовать семь независимых переменных, так как размерность пространства окажется излишне велика. Размерность семь получается в силу того, что коэффициенты представляют собой комплексные числа (1 + 2 + 2 + 2).
Поэтому искусственно снизим размерность до трех -- возмем в качестве объясняющих переменных два коэффициента с наибольшей энергией ($A_0$ и $A_3$).
Построим графики точкек:
In [30]:
result = []
for data in [water0, water1, water2,
pine_swamp0, pine_swamp1, pine_swamp2,
rostrata0, rostrata1, rostrata2,
pine0, pine1, pine2]:
f = fft(data)
result.append([f[0].real, f[3].real, f[3].imag])
result = pd.DataFrame(result, columns=['x1', 'x2', 'x3'])
# Нормализуем для удобства отображения
# result = (result - result.mean())/(result.max() - result.min())
In [31]:
# Тут я пытался построить график в 3d, но он вышел не читаемый
# поэтому я перешел на 2d -- только по первым двум координатам
# from mpl_toolkits.mplot3d import Axes3D
# fig = plt.figure()
# ax = fig.add_subplot(111, projection='3d')
# ax.scatter(result.x1[:3], result.x2[:3], result.x3[:3], marker = '*', s=36) # water
# ax.scatter(result.x1[3:6], result.x2[3:6], result.x3[3:6], marker = '^', s=36) # pine_swamp
# ax.scatter(result.x1[6:9], result.x2[6:9], result.x3[6:9], marker = '>', s=36) # rostr
# ax.scatter(result.x1[9:12], result.x2[9:12], result.x3[9:12], marker = 'v', s=36) # pine
# plt.show()
fig = plt.figure()
plot(result.x1[:3], result.x2[:3], marker = '*') # water
plot(result.x1[3:6], result.x2[3:6], marker = '^') # pine_swamp
plot(result.x1[6:9], result.x2[6:9], marker = '>') # rostr
plot(result.x1[9:12], result.x2[9:12], marker = 'v') # pine
plt.show()
Из графика видно, что точки разных классов лежат хотя и в близких, но все же различных областях пространства. Это позволяет сделать предположение о том, что метод, в основе которого лежат процедуры анализа гармоник, может с успехом использоваться для сегментации ДЗЗ.