Есть временные ряды NDVI, собранные за несколько лет. В данных есть шумы, обусловленные состоянием атмосферы, условиями съемки и т.д. Цель -- найти "модельную" кривую NDVI. Решение: строится зависимость NDVI от номера дня и методом наименьших квадратов подгоняется полином степени $d$.
In [1]:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from tools import get_ndvi_frame # Функция импорта данных ndvi
In [2]:
frame = get_ndvi_frame('data/Landsat/point_1_band3.csv', 'data/Landsat/point_1_band4.csv')
print frame.head()
plot(frame.Day, frame.NDVI, '.')
Out[2]:
Откровенно говоря, результат не впечатляет. Но все равно попробуем подогнать кривую. Для удобства расчетов центрируем наши данные перед подгонкой.
In [3]:
x = frame.Day - frame.Day.mean()
y = frame.NDVI
Начнем с параболы:
In [4]:
coef = np.polyfit(x, y, 2)
print coef
new_x = np.linspace(np.min(x), np.max(x), 100)
predict = np.poly1d(coef)
plot(new_x, predict(new_x), '-', x, y, '.')
Out[4]:
Куб:
In [5]:
coef = np.polyfit(x, y, 3)
new_x = np.linspace(np.min(x), np.max(x), 100)
predict = np.poly1d(coef)
plot(new_x, predict(new_x), '-', x, y, '.')
Out[5]:
Пятая степень
In [6]:
coef = np.polyfit(x, y, 5)
print coef
new_x = np.linspace(np.min(x), np.max(x), 100)
predict = np.poly1d(coef)
plot(new_x, predict(new_x), '-', x, y, '.')
Out[6]:
Ничего хорошего, что не удивительно при таких входных данных. Простой способ не прошел, будем действовать иначе.
Если посмотреть на график, видно, что кривая должна иметь колоколообразную форму. Те точки, которые попали во внутреннюю часть колокола -- скорее всего шумовые и на них не нужно обращать внимания. Поэтому первое, что приходит в голову --- следует построить кривую, проходящую по верхней границе облака точек. Один из способов это сделать --- последовательно пробежать по каждому дню года и определить наибольшее значение NDVI, приходящееся на этот день. Оставив в данных только максимальное значение, мы значительно сократим влияние шумовых точек, и у нас появится возможность построить аппроксимацию верхней границы. Правда, такой подход, по всей видимости, будет несколько завышать значения NDVI в начале и конце вегетационного периода (т.е. удлинять вегетационный период) за счет выбора максимального, а не среднего отклика.
In [7]:
grouped = frame.groupby('Day').max()
print grouped.head()
plot(grouped.index, grouped.NDVI, '.')
Out[7]:
Видно, что группировка и поиск максимума за один день не дает удовлетворительного результата. Придется поступать стандартным способом (например, как строятся ряды NDVI у MODIS): группировать данные за N дней и только потом находить максимум в каждой группе. Для удобства группировки добавим еще одну колонку -- номер группы:
In [8]:
size = 5 # Будем группировать по size дней
frame['GroupNum'] = (frame.Day / size).apply(np.int)
print frame.head()
In [9]:
grouped = frame.groupby('GroupNum').max()
print grouped.head()
plot(grouped.index, grouped.NDVI, '.')
Out[9]:
Это выглядит гораздо лучше, если не считать нехорошие выбросы в начале и конце периода. Попробуем повторить процедуру с подгонкой кривых.
In [10]:
x = grouped.Day - grouped.Day.mean()
y = grouped.NDVI
Как и в прошлый раз, начнем с параболы.
In [11]:
coef = np.polyfit(x, y, 2)
new_x = np.linspace(np.min(x), np.max(x), 100)
predict = np.poly1d(coef)
plot(new_x, predict(new_x), '-', x, y, '.')
Out[11]:
Куб:
In [12]:
coef = np.polyfit(x, y, 3)
new_x = np.linspace(np.min(x), np.max(x), 100)
predict = np.poly1d(coef)
plot(new_x, predict(new_x), '-', x, y, '.')
Out[12]:
Пятая степень:
In [13]:
coef = np.polyfit(x, y, 5)
new_x = np.linspace(np.min(x), np.max(x), 100)
predict = np.poly1d(coef)
plot(new_x, predict(new_x), '-', x, y, '.')
Out[13]: