Кластеризация временных серий NDVI на базе коэффициентов разложения в ряд Фурье

Ранее рассматривался вопрос о том, что данные ДЗЗ можно классифицировать не только по одномоментному снимку, но и по длинной временной серии. В частности, растительность можно попробовать классифицировать по многолетней серии NDVI: предполагается, что разные классы растительности обладают различными формами кривых NDVI. Эти формы могут выступать в качестве признаков классификации. Форму кривой можно описать, воспользовавашись преобразованием Фурье, которое представляет сигнал в виде числовых коэффциентов. Поэтому полученные коэффициенты могут служить признаками классификации.

В статье рассматривался разведочный анализ данных NDVI, из продукта MOD13 (территория в районе Ханты-Мансийска).

Здесь рассматривается вопрос кластеризации этих данных с целью выявления структуры данных и уточнения обучающих классов.

Классификация


In [1]:
import os
import sys

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import matplotlib as mpl

from scipy.fftpack import fft, ifft, fftshift, fftfreq

from FFT_NDVI_Segmentation_GRASS import GRASS

In [2]:
grass = GRASS('/opt/grass70-svn', '/home/dima/GIS/GRASSDATA', 'ugra')
data = grass.ndvi_as_array()

Устройство и назначение следующих функций описано в статье, поэтому не будем на них подробно останавливаться.


In [3]:
def get_max_coef_idx(data_fft, N, return_sym=True):
    """Возвращает индексы у N наиболее значимых гармоник.
    Если return_sym=True, то возвращать все коэффициенты (с учетом симметрии),
    как с положительными, так и отрицательными частотами. (Т.е. 1 + 2(N-1) коэффициентов )
    """
    data_fft_power = abs(data_fft)**2  # энергетический спектр
    
    size = data_fft_power.shape[-1]
    half = 1 + size/2
    
    # индексы N наибольших по модулю коэф.
    idx = argsort(data_fft_power[:, :, :half])  # Учитываем симметрию
    idx = idx[:, :, ::-1] # Нужно по убыванию
    idx = idx[:, :, :N]
    
    # Получили трехмерный массив в котором хранятся
    # номера наиболее значимых гармоник. Но у разных пикселей
    # значимые гармоники могут быть разными.
    # Отсортируем еще раз и найдем N самых важных
    u, counts = np.unique(idx, return_counts=True)  # Считаем элементы
    structure = [('counts', int), ('index', int)]   # Затем объединим их в структуру и отсортируем
    vals = np.array(zip(counts, u), dtype=structure)
    ordered = np.sort(vals, order='counts')[::-1]

    # Отрежем N самых частотных
    idx = ordered[:N]
    res = []
    # Соберем окончательный массив индексов
    for (count, u) in idx:
        res.append(u)
        if u != 0 and return_sym:    # учтем симметрию, если нужно
            res.append(size - u)
            
    return res
    

def filter_ndvi_data(data, N):
    """Возвращает копию данных к которым применены следующие
    действия:
    Преобразование Фурье.
    Зануление всех коэффициентов, кроме N значимых.
    Обратное преобразование Фурье.
    """
    
    data_fft = fft(data, axis=2)
    
    idx = get_max_coef_idx(data_fft, N)
        
    filtered = np.zeros(data_fft.shape, dtype=np.complex)
    filtered[:, :, idx] = data_fft[:, :, idx]
    
    ifiltered = ifft(filtered).real
    
    return ifiltered

Произведем кластеризацию наиболее значимых коэффициентов. Для этого сначала вычленим наиболее значимые коэффициенты.


In [4]:
data_fft = fft(data, axis=2)

N = 4
idx = get_max_coef_idx(data_fft, N, return_sym=False)
print idx
idx.sort()  # После сортировки теряем информацию о значимости
print idx

# Массив самых значимых коэффициентов
coefs = data_fft[:, :, idx]


[3, 0, 9, 6]
[0, 3, 6, 9]

Преобразуем данные в массивы, пригодные для подачи в функции кластеризации


In [5]:
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

rows, cols, features = coefs.shape
samples = rows * cols
print 'Samples:', samples

X1 = np.zeros((samples, features*2 - 1))
X = coefs.reshape((samples, features), order='F')
X1[:, 0] = np.real(X[:, 0])
for i in range(1, N):
    X1[:, 2*i - 1 ] = np.real(X[:, i])
    X1[:, 2*i] = np.imag(X[:, i])

# Нормируем данные
sc = StandardScaler()
X = sc.fit_transform(X1)


Samples: 277046

Зададим число кластеров, произведем кластеризацию.


In [6]:
n_clusters = 25

kmeans = KMeans(init='k-means++', n_clusters=n_clusters, n_init=n_clusters)
kmeans.fit(X)


Out[6]:
KMeans(copy_x=True, init='k-means++', max_iter=300, n_clusters=25, n_init=25,
    n_jobs=1, precompute_distances=True, random_state=None, tol=0.0001,
    verbose=0)

Метки кластеров и реальное число кластеров:


In [7]:
labels = kmeans.labels_

n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)

n_clusters_


Out[7]:
25

Произведем денормирование и отобразим данные в привычном виде в системе координат амплитуд синусоид.


In [8]:
# Денормирование
X = sc.inverse_transform(X)

df = pd.DataFrame({'Class': labels, 
                       'Const': X[:, 0], 
                       'Sin3.imag': X[:, 2],
                       'Sin3.real': X[:, 1],
                       'Sin9.imag': X[:, 6],
                       'Sin9.real': X[:, 5],
                       'Sin6.imag': X[:, 4],
                       'Sin6.real': X[:, 3]
                       })
df['Amp3'] = sqrt(df['Sin3.imag']**2 + df['Sin3.real']**2)
df['Amp9'] = sqrt(df['Sin9.imag']**2 + df['Sin9.real']**2)
df['Amp6'] = sqrt(df['Sin6.imag']**2 + df['Sin6.real']**2)
df.head()


Out[8]:
Class Const Sin3.imag Sin3.real Sin6.imag Sin6.real Sin9.imag Sin9.real Amp3 Amp9 Amp6
0 15 172022 34606.018796 -154991.136251 -37983.216420 32835.320556 -15275.353166 -2099.076564 158807.521400 15418.901932 50208.395769
1 15 176566 40409.758337 -153813.484429 -37977.142784 29372.351059 -13044.617987 -1138.749368 159033.130389 13094.228062 48010.398673
2 9 179187 49757.526935 -148475.440865 -40869.265554 17032.595192 -13717.571823 4900.320874 156591.085400 14566.568621 44276.474181
3 15 185151 35012.175225 -160545.595904 -34122.511964 29670.370648 -11639.798086 579.677582 164319.021352 11654.223508 45218.101653
4 15 178488 36553.700209 -149512.518784 -35471.315376 22762.907198 -12207.098662 5216.719173 153916.101406 13275.067483 42146.935340

In [9]:
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(15, 15), dpi=100)

ax = fig.add_subplot(111, projection='3d')

unique_labels = set(labels)
colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))
for k, col in zip(unique_labels, colors):
    if k == -1:
        # Черный цвет для шумовых точек
        col = 'k'

    class_member_mask = (labels == k)

    xy = df[class_member_mask]
#     plt.plot(xy.Amp3, xy.Const, 'o', markerfacecolor=col,
#              markeredgecolor='k', markersize=6)
    ax.scatter(xy.Amp3, xy.Amp9, xy.Const, color=col)

plt.title('Estimated number of clusters: %d' % n_clusters_)
plt.show()



In [11]:
fig = plt.figure(figsize=(15, 15), dpi=100)

unique_labels = set(labels)
colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))
for k, col in zip(unique_labels, colors):
    if k == -1:
        # Черный цвет для шумовых точек
        col = 'k'

    class_member_mask = (labels == k)

    xy = df[class_member_mask]
    plt.plot(xy.Amp3, xy.Const, 'o', markerfacecolor=col,
             markeredgecolor='k', markersize=6)

plt.title('Estimated number of clusters: %d' % n_clusters_)
plt.show()