Ранее рассматривался вопрос о том, что данные ДЗЗ можно классифицировать не только по одномоментному снимку, но и по длинной временной серии. В частности, растительность можно попробовать классифицировать по многолетней серии NDVI: предполагается, что разные классы растительности обладают различными формами кривых NDVI. Эти формы могут выступать в качестве признаков классификации. Форму кривой можно описать, воспользовавашись преобразованием Фурье, которое представляет сигнал в виде числовых коэффциентов. Поэтому полученные коэффициенты могут служить признаками классификации. Более подробно об этом сказано в статье, приведенной выше.
В данной статье рассматривается разведочный анализ на данных NDVI, из продукта MOD13 (территория в районе Ханты-Мансийска).
In [1]:
# coding=utf-8
import os
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
from FFT_NDVI_Segmentation_GRASS import GRASS
Перед началом работы нужно обязательно инициализировать GRASS, для этого нужно передать путь к каталогу, в который установлен GRASS, путь к каталогу БД GRASS, а также название MAPSET.
In [2]:
grass = GRASS('/opt/grass70-svn', '/home/dima/GIS/GRASSDATA', 'ugra')
Информация по текущему региону:
In [3]:
data = grass.ndvi_as_array()
In [4]:
from scipy.fftpack import fft, ifft, fftshift, fftfreq
Как было показно ранее после преобразования Фурье можно анализировать не все коэффициенты, а лишь те, которые описывают основную структуру сигнала. В этом подразделе рассматривается вопрос о том, какие коэффициенты должны быть включены в анализ, а какие могут быть опущены без потери точности.
Основная идея такова: Оставляются те коэффициенты, которые несут наибольшую энергию, а коэффициенты с энергией меньше порога обнуляются. Вопрос о том, как можно выбрать этот порог, решается следующим образом:
Самое интересное происходит в пятом пункте -- выбор функции оценки погрешности.
Построим для удобства работы вспомогательную функцию, которая будет производить шаги 2-- 4. Будем рассматривать пиксели растров по отдельности (т.е. произведем одномерное преобразование Фурье). Таким образом после преобразования Фурье у нас будет трехмерных массив комплексных чисел -- первые два измерения отвечают за широту/долготу пикселя, а третье измерение -- коэффициенты разложения.
In [5]:
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 [6]:
N = 15 # Большое, чтобы были видны характерные детали
filtered = filter_ndvi_data(data, N)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(data[10, 20, :], '-bo', label='Data')
plt.plot(filtered[10, 20, :], '-ro', label='Filtered')
plt.legend()
plt.show()
Можно взять массив data и прогнать его через функцию filter_ndvi_data, а потом сравнить полученный сигнал с исходным. Разность исходного и восстановленного сигналов может служить фукнцией качества. Однако очевидно, что:
Поэтому подберем $N$ на основе процедуры перекрестной проверки. В обучающее множество положим отсчеты с четными индексами, в тестовое -- с нечетными. (Такой подход к разбиению обусловлен упрощением кода и вычислений -- для преобразования Фурье нужно, чтобы отсчеты были равноотстоящими, а изъятие данных случайным образом значительно затруднит вычисления, поэтому разбиение произведем систематически).
После разбиения на обучающем множестве прогоним функцию filter_ndvi_data и получим сглаженные кривые NDVI для каждого пикселя. Затем произведем линейную интерполяцию полученных кривых и проинтерполированные значения со значениями, лежащих в тестовом множестве. Величина погрешности на тестовом множестве и будет оценкой качества.
In [7]:
# Функция разбиения на тестовое и обучающее множество
def split_data(data):
size = data.shape[-1]
train = data[:, :, 0:size:2]
test = data[:, :, 1:size:2]
# Если число отсчетов в data четное,
# то в test один остчет "торчит" наружу
# (за пределы области изменения обучающего множества)
# поэтому удалим последний отсчет в теством множестве
if len(train[0, 0, :]) == len(test[0, 0, :]):
test = test[:, :, :-1]
return (train.copy(), test.copy())
Покажем на примере одного пикселя, как будем считать ошибку фильтрации:
In [8]:
train, test = split_data(data)
N = 6
size = data.shape[-1]
filtered = filter_ndvi_data(train, N)
# Индексы отсчетов обучающего множества: 0, 2, ..., N
t_train = np.arange(0, size, 2)
# Индексы отсчетов тестового множества: 1, 3, ..., N-1
t_test = t_train[:-1].copy()
t_test += 1
# Проинтерполируем отфильтрованные результаты в точках, попавших в тествоем множество
recieved = np.interp(t_test, t_train, filtered[20, 20, :])
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(t_train, filtered[20, 20, :], '-go', label='Filtered')
plt.plot(t_train, train[20, 20, :], '-bo', label='Train')
plt.plot(t_test, test[20, 20, :], '-ro', label='Test')
plt.plot(t_test, recieved, '^', label='Interpolated filtered')
plt.legend()
plt.show()
# Считаем ошибку
from sklearn.metrics import mean_squared_error
err = mean_squared_error(test[20, 20, :], recieved)
print err
Оформим приведенный код в виде функции:
In [9]:
def estimate_error(train, test, N):
filtered = filter_ndvi_data(train, N)
# Индексы отсчетов обучающего множества: 0, 2, ..., N
t_train = np.arange(0, size, 2)
# Индексы отсчетов тестового множества: 1, 3, ..., N-1
t_test = t_train[:-1].copy()
t_test += 1
err = 0
dim_x, dim_y, dim_ndvi = train.shape
recieved = np.zeros(test.shape)
# TODO: этот цикл очень неэффективный. Исправить на использование numpy
for i in range(dim_x):
for j in range(dim_y):
recieved[i, j, :] = np.interp(t_test, t_train, filtered[i, j, :])
err += mean_squared_error(test[i, j, :], recieved[i, j])
return err
Находим "нужное" N. Код ниже закомментирован, поскольку выполняется довольно долго.
In [10]:
# size = data.shape[-1]
# train, test = split_data(data)
# errors = np.zeros(size/2 -1)
# errors[:] = inf
# for N in range(1, size/2):
# cur_err = estimate_error(train, test, N)
# print N, cur_err
# errors[N-1] = cur_err
# fig = plt.figure(figsize=(12,6), dpi=100)
# plt.plot(range(2, size/2), errors[1:], '-o', label='Errors') # ошибка при N=1 слишком велика
# plt.legend()
# plt.show()
Из графика видно, что минимум ошибки достигается при N=4. При этом падение ошибки идет резко при N=1 и N=2, после N=2 (N=3, N=4) падение ошибки незначительное.
Как видно из фрагмента вывода ошибок, при классификации нужно будет использовать не менее двух коэффициентов. Оптимум получим при четырех, но он незначительно лучше, чем при двух. Поэтому мы можем с чистой совестью использовать 2 коэффициента. Тем не менее сохраним N=4 наиболее значимых коэффициентов в виде массивов (а вдруг пригодятся):
In [11]:
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]
Сохраним полученные массивы как растры GRASS.
Коэффициент $A_0$ (константа) является действительным, поэтому обработаем его отдельно:
In [12]:
band = grass.garray.array()
band[...] = coefs[:, :, 0].real
band.write('Const', overwrite=True)
Out[12]:
Сохраним оставшиеся коэффициенты:
In [13]:
for i in range(1, 4):
band = grass.garray.array()
band[...] = coefs[:, :, i].real
band.write('Sin' + str(idx[i])+'.real', overwrite=True)
band[...] = coefs[:, :, i].imag
band.write('Sin' + str(idx[i])+'.image', overwrite=True)
Прочитаем растры с тренировочными данными. В наборе имеется три растра train1, train2 и train3, которые соотвествуют обучающим классам на трех уровнях иерархии. В таблице приведены значения, которыми закодированы классы:
| № класса в train1 | название класса в train1 | № класса в train2 | название класса в train2 | № класса в train3 | название класса в train3 | Сложный для распознавания? | Коментарий | Редкий? |
|---|---|---|---|---|---|---|---|---|
| 1 | Ледники постоянные | 1 | Ледники постоянные | 1 | Ледники постоянные | |||
| 2 | Вода открытая | 2 | Вода «синяя» | 2 | Вода «синяя» | |||
| 2 | Вода открытая | 3 | Вода «черная» | 3 | Вода «черная» | |||
| 3 | Лишенные растительности участки | 4 | Пески | 4 | Пески | |||
| 3 | Лишенные растительности участки | 5 | Застройка | 5 | Застройка | |||
| 4 | Пойменные луга | 6 | Пойменные луга | 6 | Пойменные луга | |||
| 5 | Болота | 7 | Олиготрофные | 10 | Грядово-мочажинные | |||
| 5 | Болота | 7 | Олиготрофные | 11 | Рямы | |||
| 5 | Болота | 7 | Олиготрофные | 12 | Плоскобугристые | 1 | Близки с классом лишайниковых тундр | |
| 5 | Болота | 7 | Олиготрофные | 13 | Грядово-озерковые и ГТК | |||
| 5 | Болота | 8 | Мезотрофные | 20 | Вахтово-сфагновые | 1 | Смешиваются с ГМК (5-7-10) | Да |
| 5 | Болота | 8 | Мезотрофные | 21 | Гипновые топи (лазиокарпа) | 1 | Плохо отделяются (смешиваются с несколькими классами) | |
| 5 | Болота | 8 | Мезотрофные | 22 | Хасыреи | 1 | Плохо отделяются (смешиваются с несколькими классами) | |
| 5 | Болота | 9 | Мезотрофные | 31 | Разнотравно-гипновые богаого грунтового питания | Плохо отделяются (смешиваются с несколькими классами) | Да | |
| 5 | Болота | 9 | Эутрофные | 30 | Согры | 1 | Смешиваются с темнохвойными лесами | |
| 6 | Леса | 10 | Лиственные | 40 | Лиственные | |||
| 6 | Леса | 11 | Смешанные | 41 | Смешанные | |||
| 6 | Леса | 12 | Хвойные (не беломошные) | 42 | Темнохвойные | |||
| 6 | Леса | 12 | Хвойные (не беломошные) | 44 | Сосново-кустарничково-моховые (редкостойные) | 1 | Смешиваются с рямами | |
| 6 | Леса | 13 | Беломошные | 43 | Сосновые и лиственничные беломошные | |||
| 6 | Леса | 14 | Гари (возраст примерно 10 и более лет) | 45 | Гари | |||
| 7 | Тундры | 15 | Тундры горные | 50 | Тундры горные | |||
| 7 | Тундры | 16 | Тундры равнинные | 51 | Тундры равнинные | |||
| 8 | Пашня | 17 | Пашня | 70 | Пашня |
Прочитаем обучающие растры и проверим, какие классы присутствуют в данных. (Класс номер 0 -- отсутствие информации).
In [14]:
train1 = grass.garray.array()
train1.read('train1')
train2 = grass.garray.array()
train2.read('train2')
train3 = grass.garray.array()
train3.read('train3')
print 'Уникальные классы в train1:', np.unique(train1)
print 'Уникальные классы в train2:', np.unique(train2)
print 'Уникальные классы в train3:', np.unique(train3)
Индексы пикселей, для которых известны номера классов:
In [15]:
rows1, cols1 = np.nonzero(train1)
rows2, cols2 = np.nonzero(train2)
rows3, cols3 = np.nonzero(train1)
Рассмотрим сначала классы на первом уровне иерархии. Возьмем коэффициенты разложения Фурье и сравним их для разных классов.
In [16]:
data_fft = fft(data, axis=2)
N = 4
idx = get_max_coef_idx(data_fft, N, return_sym=False)
print idx
coefs1 = data_fft[:, :, idx]
Создадим фрейм данных.
In [17]:
labels = {2: 'Water', 3: 'Non-vegetation', 4: 'Grassland', 5: 'Swamp', 6: 'Forest'}
cats = pd.Series([labels[cat] for cat in train1[rows1, cols1]], dtype="category")
df1 = pd.DataFrame({'Class': train1[rows1, cols1],
'Names': cats,
'Const': coefs1[rows1, cols1, 1].real,
'Sin3.imag': coefs1[rows1, cols1, 0].imag,
'Sin3.real': coefs1[rows1, cols1, 0].real,
'Sin9.imag': coefs1[rows1, cols1, 2].imag,
'Sin9.real': coefs1[rows1, cols1, 2].real,
'Sin6.imag': coefs1[rows1, cols1, 3].imag,
'Sin6.real': coefs1[rows1, cols1, 3].real,
})
df1.head()
Out[17]:
In [18]:
df1.hist(column='Class', bins=50)
Out[18]:
Видно, что классы представлены неравномерно. В частности для третьего класса (участков, лишенных растительности) данных очень мало. Большее всего примеров на лесные территории.
Построим "ящики с усами".
In [19]:
df1.boxplot(column='Const', by='Names')
Out[19]:
Из диаграммы, построенной по одному коэффициенту (константе, которая может интерпретироватсья как среднее значение NDVI), видно:
Посмотрим, что происходит при использовании синусоид в качестве объясняющих переменных.
In [20]:
df1.boxplot(column='Sin3.real', by='Names')
Out[20]:
In [21]:
df1.boxplot(column='Sin3.imag', by='Names')
Out[21]:
Сами по себе коэффициенты мало что привнесли к тому, что мы уже знаем. Посмотрим, что дадут амплитуды синусоид. Для этого добавим в таблицу еще одну колонку -- апмлитуд.
In [22]:
df1['Amp3'] = sqrt(df1['Sin3.imag']**2 + df1['Sin3.real']**2)
df1['Amp9'] = sqrt(df1['Sin9.imag']**2 + df1['Sin9.real']**2)
df1['Amp6'] = sqrt(df1['Sin6.imag']**2 + df1['Sin6.real']**2)
df1.head()
Out[22]:
In [23]:
df1.boxplot(column='Amp3', by='Names')
Out[23]:
По диаграмме амплитуд получаем достаточно предсказуемые результаты -- амплитуды синусоид на территориях, покрытых растительностью, значительно выше, чем амплитуды синусоид у водных объектов и участков, лишенных растительности.
Построим диаграмму распределения классов в зависимости от амплитуды синусоиды и свободгого члена: по оси абсцисс отложим значения константы из разложения Фурье, а по оси ординат -- амплитуды синуса с "тремя горбами". Эта диаграмма относительно легко интерпретируется в "физических" величинах. В первом приближении константа равна среднему значению NDVI за анализируемый период, а амплитуда синуса -- амплитуде изменения NDVI в течении сезона (мы анализируем трехлетний период, поэтому и рассматриваем на диаграмме синус с тремя периодами, тем более, как было показано выше, этот синус имеет наибольшую предсказательную способность среди остальных коэффициентов).
In [24]:
water = df1[df1['Class'] == 2]
forest = df1[df1['Class'] == 6]
non_veget = df1[df1['Class'] == 3]
grass_land = df1[df1['Class'] == 4]
swamp = df1[df1['Class'] == 5]
fig = plt.figure(figsize=(8, 8), dpi=100)
plt.plot(water['Amp3'], water['Const'], 'r+', label='Water')
plt.plot(forest['Amp3'], forest['Const'], 'g+', label='Forest')
plt.plot(non_veget['Amp3'], non_veget['Const'], 'b*', label='Non vegetation')
plt.plot(grass_land['Amp3'], grass_land['Const'], 'r*', label='Grasslands')
plt.plot(swamp['Amp3'], swamp['Const'], 'b+', label='Swamp')
plt.legend()
plt.ylabel('Const')
plt.xlabel('Amp3')
plt.show()
Дополнительно к тому, что мы уже узнали ранее, из диаграммы видим:
Если посмотреть диаграммы по другим синусам, то каких-то новых результатов увидеть не удается.
In [25]:
labels = {1: 'Glaciers',
2: 'Blue water',
3: 'Black water',
4: 'Sands',
5: 'Build',
6: 'Grass',
7: 'Oligotr.',
8: 'Mesotr.',
9: 'Eutrofnye',
10: 'Hardwood',
11: 'Mixed f.',
12: 'Con.(not wh-moss)',
13: 'Wh-moss Conifer',
14: 'Burned',
15: 'Mount.tundra',
16: 'Pl.tundra',
17: 'Arable land'}
train2 = grass.garray.array()
train2.read('train2')
cats = pd.Series([labels[cat] for cat in train2[rows1, cols1]], dtype="category")
df2 = pd.DataFrame({'Class': train2[rows1, cols1],
'Names': cats,
'Const': coefs1[rows1, cols1, 1].real,
'Sin3.imag': coefs1[rows1, cols1, 0].imag,
'Sin3.real': coefs1[rows1, cols1, 0].real,
'Sin9.imag': coefs1[rows1, cols1, 2].imag,
'Sin9.real': coefs1[rows1, cols1, 2].real,
'Sin6.imag': coefs1[rows1, cols1, 3].imag,
'Sin6.real': coefs1[rows1, cols1, 3].real,
})
df2['Amp3'] = sqrt(df2['Sin3.imag']**2 + df2['Sin3.real']**2)
df2['Amp9'] = sqrt(df2['Sin9.imag']**2 + df2['Sin9.real']**2)
df2['Amp6'] = sqrt(df2['Sin6.imag']**2 + df2['Sin6.real']**2)
df2.head()
Out[25]:
In [26]:
df2.hist(column='Class', bins=50)
Out[26]:
Как и ранее видим, что количество примеров в разных классах сильно варьируется.
In [27]:
df2.boxplot(column='Const', by='Names', figsize=(15,10))
Out[27]:
In [28]:
df2.boxplot(column='Amp3', by='Names', figsize=(15,10))
Out[28]:
In [29]:
blue_w = df2[df2['Class'] == 2]
black_w = df2[df2['Class'] == 3]
build = df2[df2['Class'] == 5]
grassland = df2[df2['Class'] == 6]
ologotr = df2[df2['Class'] == 7]
mesotr = df2[df2['Class'] == 8]
eutr = df2[df2['Class'] == 9]
hardw = df2[df2['Class'] == 10]
mixed = df2[df2['Class'] == 11]
non_wh_mos_con = df2[df2['Class'] == 12]
wh_moss_c = df2[df2['Class'] == 13]
burned = df2[df2['Class'] == 14]
fig = plt.figure(figsize=(8, 8), dpi=100)
plt.plot(blue_w['Amp3'], blue_w['Const'], 'b+', label='Blue Water')
plt.plot(black_w['Amp3'], black_w['Const'], 'k+', label='Black Water')
plt.plot(build['Amp3'], build['Const'], 'y^', label='Buildins')
plt.plot(grassland['Amp3'], grassland['Const'], 'g*', label='Grass')
plt.plot(ologotr['Amp3'], ologotr['Const'], 'rv', label='Oligotrophic')
plt.plot(mesotr['Amp3'], mesotr['Const'], 'b*', label='Mesotrophic')
plt.plot(eutr['Amp3'], eutr['Const'], 'c^', label='Eutrofic')
plt.plot(hardw['Amp3'], hardw['Const'], 'ro', label='Hardwood')
plt.plot(mixed['Amp3'], mixed['Const'], 'g+', label='Mixed forest')
plt.plot(non_wh_mos_con['Amp3'], non_wh_mos_con['Const'], 'k^', label='Conifer (not white moss)')
plt.plot(wh_moss_c['Amp3'], wh_moss_c['Const'], 'mo', label='Conifer (white moss)')
plt.plot(burned['Amp3'], burned['Const'], 'y+', label='Burned')
plt.ylabel('Const')
plt.xlabel('Amp3')
plt.legend(loc=2)
plt.show()
Для того, чтобы разглядеть детали ниже на диаграмме показаны болота и леса (за исключением хвойных лесов).
In [30]:
fig = plt.figure(figsize=(8, 8), dpi=100)
plt.plot(ologotr['Amp3'], ologotr['Const'], 'rv', label='Oligotrophic')
plt.plot(mesotr['Amp3'], mesotr['Const'], 'b*', label='Mesotrophic')
plt.plot(eutr['Amp3'], eutr['Const'], 'c^', label='Eutrofic')
plt.plot(hardw['Amp3'], hardw['Const'], 'ro', label='Hardwood')
plt.plot(mixed['Amp3'], mixed['Const'], 'g+', label='Mixed forest')
plt.plot(burned['Amp3'], burned['Const'], 'y+', label='Burned')
plt.legend(loc=0)
plt.ylabel('Const')
plt.xlabel('Amp3')
plt.show()
На третьем уровне иерархии нет смысла рассматривать все классы, поскольку часть из них не изменилась. Изменению подверглись лишь лесные территории и болота. Поэтому проведем анализ только этих классов, причем остановимся на тех, которые помечены в документации как плохо отделяющиеся от других классов.
В описании обучающих данных два класса под номером 43. Скорее всего ошибка и класс "Сосновые и лиственничные беломошные леса" должен идти под номером 44.
In [31]:
labels = {
10: 'Ridge-hollow',
11: 'Ryam',
12: 'Plain-hillock',
13: 'Ridge-lake',
20: 'Shift-sphagnum',
21: 'Sedge marshes',
22: 'Hasyrei',
31: 'Herb-rich',
30: 'Sogra',
40: 'Hardwood',
41: 'Mixed',
42: 'Dark coniferous',
44: 'Pine-shrub-moss',
43: 'Pine and larch wh-moss',
45: 'Burned'
}
train3 = grass.garray.array()
train3.read('train3')
cats = pd.Series([labels[cat] if cat in labels else '-' for cat in train3[rows1, cols1]], dtype="category")
df3 = pd.DataFrame({'Class': train3[rows1, cols1],
'Names': cats,
'Const': coefs1[rows1, cols1, 1].real,
'Sin3.imag': coefs1[rows1, cols1, 0].imag,
'Sin3.real': coefs1[rows1, cols1, 0].real,
'Sin9.imag': coefs1[rows1, cols1, 2].imag,
'Sin9.real': coefs1[rows1, cols1, 2].real,
'Sin6.imag': coefs1[rows1, cols1, 3].imag,
'Sin6.real': coefs1[rows1, cols1, 3].real,
})
df3 = df3[df3['Names'] != '-']
df3['Amp3'] = sqrt(df3['Sin3.imag']**2 + df3['Sin3.real']**2)
df3['Amp9'] = sqrt(df3['Sin9.imag']**2 + df3['Sin9.real']**2)
df3['Amp6'] = sqrt(df3['Sin6.imag']**2 + df3['Sin6.real']**2)
df3.head()
Out[31]:
In [32]:
df3.hist(column='Class', bins=50)
Out[32]:
In [33]:
labels = {
10: 'Ridge-hollow',
11: 'Ryam',
12: 'Plain-hillock',
13: 'Ridge-lake',
20: 'Shift-sphagnum',
21: 'Sedge marshes',
22: 'Hasyrei',
31: 'Herb-rich',
30: 'Sogra',
40: 'Hardwood',
41: 'Mixed',
42: 'Dark coniferous',
44: 'Pine-shrub-moss',
43: 'Pine and larch wh-moss',
45: 'Burned'
}
r_h = df3[df3['Class'] == 10]
ryam = df3[df3['Class'] == 11]
p_h = df3[df3['Class'] == 12]
r_l = df3[df3['Class'] == 13]
s_s = df3[df3['Class'] == 20]
s_m = df3[df3['Class'] == 21]
haserey = df3[df3['Class'] == 22]
sogra = df3[df3['Class'] == 30]
h_r = df3[df3['Class'] == 31]
hardw = df3[df3['Class'] == 40]
mixed = df3[df3['Class'] == 41]
dark = df3[df3['Class'] == 42]
p_s_m = df3[df3['Class'] == 43]
p_l_m = df3[df3['Class'] == 44] # Проверить номер класса
burned = df3[df3['Class'] == 45]
fig = plt.figure(figsize=(13, 13), dpi=100)
plt.plot(r_h['Amp3'], r_h['Const'], 'b+', label='Ridge-hollow (10)')
plt.plot(r_l['Amp3'], r_l['Const'], 'k+', label='Ridge-lake (13)')
plt.plot(s_s['Amp3'], s_s['Const'], 'g+', label='Shift-sphagnum (20)')
plt.plot(s_m['Amp3'], s_m['Const'], 'c+', label='Sedge marshes (21)')
plt.plot(h_r['Amp3'], h_r['Const'], 'mv', label='Herb-rich (31)')
plt.plot(hardw['Amp3'], hardw['Const'], 'y+', label='Hardwood (40)')
plt.plot(mixed['Amp3'], mixed['Const'], 'r^', label='Mixed forest (41)')
plt.plot(dark['Amp3'], dark['Const'], 'b^', label='Dark coniferous (42)')
plt.plot(p_s_m['Amp3'], p_s_m['Const'], 'g^', label='Pine-shrub-moss (43)')
plt.plot(p_l_m['Amp3'], p_l_m['Const'], 'ko', label='?Pine and larch wh-moss? (44)')
plt.plot(burned['Amp3'], burned['Const'], 'c^', label='Burned (45)')
# Нет в обучающией выборке:
# plt.plot(ryam['Amp3'], ryam['Const'], 'go', label='Ryam (11)')
# plt.plot(p_h['Amp3'], p_h['Const'], 'yo', label='Plain-hillock (12)')
# plt.plot(haserey['Amp3'], haserey['Const'], 'yo', label='Haserey (22)')
# plt.plot(sogra['Amp3'], sogra['Const'], 'cv', label='Sogra (30)')
plt.ylabel('Const')
plt.xlabel('Amp3')
plt.legend(loc=1)
plt.show()
Рассмотрим следующие классы, которые помечены как сложные для распознавания. При этом были выброшены из анализа классы 12, 22, 30, поскольку они не присутствуют в тренировочных данных
| № класса в train1 | название класса в train1 | № класса в train2 | название класса в train2 | № класса в train3 | название класса в train3 | Коментарий |
|---|---|---|---|---|---|---|
| 5 | Болота | 8 | Мезотрофные | 20 | Вахтово-сфагновые | Смешиваются с ГМК (5-7-10) |
| 5 | Болота | 8 | Мезотрофные | 21 | Гипновые топи (лазиокарпа) | Плохо отделяются (смешиваются с несколькими классами) |
| 5 | Болота | 9 | Мезотрофные | 31 | Разнотравно-гипновые богатого грунтового питания | Плохо отделяются (смешиваются с несколькими классами) |
| 6 | Леса | 12 | Хвойные (не беломошные) | 43 | Сосново-кустарничково-моховые (редкостойные) | Смешиваются с рямами |
In [34]:
fig = plt.figure(figsize=(8, 8), dpi=100)
plt.plot(blue_w['Amp3'], blue_w['Const'], 'b+', label='Blue Water')
plt.plot(black_w['Amp3'], black_w['Const'], 'k+', label='Black Water')
plt.plot(build['Amp3'], build['Const'], 'y^', label='Buildins')
plt.plot(grassland['Amp3'], grassland['Const'], 'g*', label='Grass')
plt.plot(ologotr['Amp3'], ologotr['Const'], 'rv', label='Oligotrophic')
plt.plot(mesotr['Amp3'], mesotr['Const'], 'b*', label='Mesotrophic')
plt.plot(eutr['Amp3'], eutr['Const'], 'c^', label='Eutrofic')
plt.plot(hardw['Amp3'], hardw['Const'], 'ro', label='Hardwood')
plt.plot(mixed['Amp3'], mixed['Const'], 'g+', label='Mixed forest')
plt.plot(non_wh_mos_con['Amp3'], non_wh_mos_con['Const'], 'k^', label='Conifer (not white moss)')
plt.plot(wh_moss_c['Amp3'], wh_moss_c['Const'], 'mo', label='Conifer (white moss)')
plt.plot(burned['Amp3'], burned['Const'], 'y+', label='Burned')
plt.plot([0, 180000], [-35000, 235000], 'k-')
plt.ylabel('Const')
plt.xlabel('Amp3')
plt.legend(loc=2)
plt.show()
Если посмотреть на диаграммы точек (например, тут), то бросается в глаза, что часть классов резко "обрезана" по нижней границе -- существует какая-то прямая линия, за которую не точки не выходят (кроме синей воды).
Эта линия интерпертируется следующим образом:
Найдем приблизительной уравнение прямой и вычтем его из данных. Это позволит "нормализировать" уровень NDVI для классов.
Уравнение зависимости константы от амлитуды:
In [35]:
def amp(c):
x1, x2 = 0, 180000
y1, y2 = -35000, 235000
return (c-x1)*(y2-y1)/(x2-x1) + y1
Собственно сам график.
In [36]:
df3['Const_a'] = df3['Const'] - amp(df3['Amp3'])
grassland = df3[df3['Class'] == 6]
r_h = df3[df3['Class'] == 10]
ryam = df3[df3['Class'] == 11]
p_h = df3[df3['Class'] == 12]
r_l = df3[df3['Class'] == 13]
s_s = df3[df3['Class'] == 20]
s_m = df3[df3['Class'] == 21]
haserey = df3[df3['Class'] == 22]
sogra = df3[df3['Class'] == 30]
h_r = df3[df3['Class'] == 31]
hardw = df3[df3['Class'] == 40]
mixed = df3[df3['Class'] == 41]
dark = df3[df3['Class'] == 42]
p_s_m = df3[df3['Class'] == 44]
p_l_m = df3[df3['Class'] == 43]
burned = df3[df3['Class'] == 45]
fig = plt.figure(figsize=(10, 10), dpi=100)
plt.plot(r_h['Amp3'], r_h['Const_a'], 'b+', label='Ridge-hollow (10)')
plt.plot(r_l['Amp3'], r_l['Const_a'], 'k+', label='Ridge-lake (13)')
plt.plot(s_s['Amp3'], s_s['Const_a'], 'g+', label='Shift-sphagnum (20)')
plt.plot(s_m['Amp3'], s_m['Const_a'], 'c+', label='Sedge marshes (21)')
plt.plot(h_r['Amp3'], h_r['Const_a'], 'mv', label='Herb-rich (31)')
plt.plot(hardw['Amp3'], hardw['Const_a'], 'y+', label='Hardwood (40)')
plt.plot(mixed['Amp3'], mixed['Const_a'], 'r^', label='Mixed forest (41)')
plt.plot(dark['Amp3'], dark['Const_a'], 'b^', label='Dark coniferous (42)')
plt.plot(p_s_m['Amp3'], p_s_m['Const_a'], 'g^', label='Pine-shrub-moss (43)')
plt.plot(p_l_m['Amp3'], p_l_m['Const_a'], 'ko', label='?Pine and larch wh-moss? (44)')
plt.plot(burned['Amp3'], burned['Const_a'], 'c^', label='Burned (45)')
plt.plot(grassland['Amp3'], grassland['Const_a'], 'g*', label='Grass')
# Нет в обучающией выборке:
# plt.plot(ryam['Amp3'], ryam['Const'], 'go', label='Ryam (11)')
# plt.plot(p_h['Amp3'], p_h['Const'], 'yo', label='Plain-hillock (12)')
# plt.plot(haserey['Amp3'], haserey['Const'], 'yo', label='Haserey (22)')
# plt.plot(sogra['Amp3'], sogra['Const'], 'cv', label='Sogra (30)')
plt.ylabel('Const (adjst)')
plt.xlabel('Amp3')
plt.legend(loc=1)
plt.show()
Такая нормализованная диаграмма позволяет увидеть картину несколько более детально, в первую очередь она позволяет сравнивать уровни минимальные NDVI разных классов (видимо, уровни NDVI в зимний период?). После нормализации горизонтальные линии соотвествуют примерно одинаковым минимальным значениям NDVI, так, на нижнем уровне диаграммы расположилось горизонтальное облако точек пойменных лугов оно лежит ниже всех остальных классов, вероятно потому, что в зимний период луга полностью покрыты снегом. Часть классов приобрела четко выраженную горизонтальную направленность, например, вахтово-сфагновые болота, гипновые топи, грядово-озерковые болота. Беломошные леса также вытянулись горизонтально.
Тот факт, что, к примеру, гипновые топи лежат на диаграмме выше, чем пойменные луга позволяет предположить, что они покрываются снегом не полностью, и на них остается что-то, повышает уровень NDVI (ветки деревьев/кустов? хвойные?). То, что основной кластер грядово-мочажинных болот лежит на диаграмме еще выше, также должно иметь объяснение (еще больше веток?).
In [36]: