Ранее рассматривался вопрос о том, что данные ДЗЗ можно классифицировать не только по одномоментному снимку, но и по длинной временной серии. В частности, растительность можно попробовать классифицировать по многолетней серии 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]
Преобразуем данные в массивы, пригодные для подачи в функции кластеризации
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)
Зададим число кластеров, произведем кластеризацию.
In [6]:
n_clusters = 25
kmeans = KMeans(init='k-means++', n_clusters=n_clusters, n_init=n_clusters)
kmeans.fit(X)
Out[6]:
Метки кластеров и реальное число кластеров:
In [7]:
labels = kmeans.labels_
n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
n_clusters_
Out[7]:
Произведем денормирование и отобразим данные в привычном виде в системе координат амплитуд синусоид.
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]:
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()