Здесь содержится код вспомогательных функций и пояснения к их работе:

  1. Функция get_ndvi_frame и пояснения к ее реализации.
  2. Функция get_modis_ndvi_frame.

Проблемы. Здесь нехорошо то, что из этого файла будут импортироваться код функций. При этом весь код, поясняющий работу этих функций будет выполняться в момент импорта. Нужно придумать способ спрятать поясняющий код, чтобы при импорте он не вызывался.

Инициализация

Импортируем необходимые библиотеки и создадим основыне функции, которые понадобятся в дальнейшем:


In [1]:
import datetime

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt

Чтение данных Landsat и расчет NDVI

Пояснения к реализации

Прочитаем табличные данные, и рассчитаем значения NDVI. Для примера возьмем точку номер 1.


In [2]:
url = 'data/Landsat/point_1_band3.csv'
frame3 = pd.io.parsers.read_csv(url, sep=';', na_values='null')
frame3.head()


Out[2]:
X Y Name Sensor Year Day Value
0 572677.37 7016668.9 toar_LT51590152009293KIS00_3 LT5 2009 293 NaN
1 572677.37 7016668.9 toar_LT51590161987073XXX01_3 LT5 1987 73 0.63195
2 572677.37 7016668.9 toar_LT51620151985264KIS00_3 LT5 1985 264 NaN
3 572677.37 7016668.9 toar_LT51600151989133KIS00_3 LT5 1989 133 NaN
4 572677.37 7016668.9 toar_LT51600151986141XXX03_3 LT5 1986 141 NaN

In [3]:
# Отрежем номер канала в названии растра и используем название как 
# индекс строки (будет удобно склеивать несколько файлов в один)
index = frame3.Name.apply(lambda x: x[:-2])

# Удалим лишние колонки
frame3 = frame3.drop(['X', 'Y', 'Sensor', 'Name'], axis=1)

frame3.index = index
frame3.head()


Out[3]:
Year Day Value
Name
toar_LT51590152009293KIS00 2009 293 NaN
toar_LT51590161987073XXX01 1987 73 0.63195
toar_LT51620151985264KIS00 1985 264 NaN
toar_LT51600151989133KIS00 1989 133 NaN
toar_LT51600151986141XXX03 1986 141 NaN

Повторим процедуру с 4-м каналом:


In [4]:
url = 'data/Landsat/point_1_band4.csv'
frame4 = pd.io.parsers.read_csv(url, sep=';', na_values='null')
index = frame4.Name.apply(lambda x: x[:-2])
frame4 = frame4.drop(['X', 'Y', 'Sensor', 'Name'], axis=1)
frame4.index = index
frame4.head()


Out[4]:
Year Day Value
Name
toar_LT51610151987103XXX02 1987 103 NaN
toar_LT51600162007183IKR01 2007 183 0.367236
toar_LT51600161990088KIS00 1990 88 0.711419
toar_LT51610161994186KIS00 1994 186 0.311629
toar_LT51580161986271XXX02 1986 271 NaN

Расчитаем NDVI, для этого возьмем данные из 4-го канала и поместим их в таблицу с данными по 3-му каналу. Затем удалим все строки, в которых нет данных по отражательной способности.


In [5]:
val4 = frame4.Value
frame = pd.concat([frame3, val4], axis=1, )
frame.columns = ['Year', 'Day', 'B3', 'B4']    # Переименуем столбцы

frame = frame.dropna()
frame.head()


Out[5]:
Year Day B3 B4
toar_LT41590161988100XXX03 1988 100 0.614420 0.677441
toar_LT41590161988180XXX03 1988 180 0.078024 0.277332
toar_LT41590161988196XXX03 1988 196 0.085931 0.330386
toar_LT41590161988212XXX03 1988 212 0.354359 0.451560
toar_LT41590161988244XXX03 1988 244 0.433639 0.547898

In [6]:
frame['NDVI'] = (frame.B4 - frame.B3)/(frame.B4 + frame.B3)
frame.head()


Out[6]:
Year Day B3 B4 NDVI
toar_LT41590161988100XXX03 1988 100 0.614420 0.677441 0.048782
toar_LT41590161988180XXX03 1988 180 0.078024 0.277332 0.560867
toar_LT41590161988196XXX03 1988 196 0.085931 0.330386 0.587186
toar_LT41590161988212XXX03 1988 212 0.354359 0.451560 0.120610
toar_LT41590161988244XXX03 1988 244 0.433639 0.547898 0.116409

Отортируем данные по году и дню, построим график NDVI


In [7]:
frame = frame.sort(columns=['Year', 'Day'])
# plt.plot(frame.NDVI)

Обработка дат

Предыдущий график показывает динамику NDVI совершенно неверно: он не учитывает дату съемки, поскольку по оси абсцисс отложены не дни, а номера строчек отсортированной таблицы с данными. Для того, чтобы график был похож на правду, нужно по оси абсцисс откладывать дни. Для этого переведем даты съемки из формата "(год, день)" в формат "номер дня", при этом номер дня будет отсчитываться от какой-либо фиксированной даты. Для примера возьмем дату 1 января 1 года (можно выбрать любую другую) и переведем все даты съемки в новый формат.

Для удобства напишем вспомогательную функцию day_number(year, day) которая по году и дню выведет номер дня, отсчитывая его от 1 января 1 года


In [8]:
def day_number(year, day):
    # Номер первого января заданного года
    d = datetime.date(year, 1, 1)
    d = d.toordinal()
    return d + day - 1

# Проверка
# print '1-е января первого года:', day_number(1, 1)
# print '3-е января второго года:', day_number(2, 3)
# print '1-е февраля 1984 года:', day_number(1984, 32)

# Эта функция служит для обработки frame:
def day_number(row):
    d = datetime.date(int(row['Year']), 1, 1)
    d = d.toordinal()
    return d + row['Day'] - 1

Добавим колонку с номером дня.


In [9]:
frame['DayNum'] = frame.apply(day_number, axis=1)
frame.head()


Out[9]:
Year Day B3 B4 NDVI DayNum
toar_LT51610161984127KIS00 1984 127 0.653173 0.618082 -0.027603 724402
toar_LT51590161984129KIS00 1984 129 0.577255 0.540409 -0.032967 724404
toar_LT51600161984184XXX05 1984 184 0.075905 0.285286 0.579696 724459
toar_LT51610161985065KIS00 1985 65 0.474825 0.523001 0.048281 724706
toar_LT51590161985067KIS00 1985 67 0.588754 0.646601 0.046826 724708

Код функции


In [10]:
def get_ndvi_frame(band3_name, band4_name):
    
    def day_number(row):
        d = datetime.date(int(row['Year']), 1, 1)
        d = d.toordinal()
        return d + row['Day'] - 1
    
    frame3 = pd.io.parsers.read_csv(band3_name, sep=';', na_values='null')
    
    index = frame3.Name.apply(lambda x: x[:-2])
    frame3 = frame3.drop(['X', 'Y', 'Sensor', 'Name'], axis=1)
    frame3.index = index
    
    frame4 = pd.io.parsers.read_csv(band4_name, sep=';', na_values='null')
    index = frame4.Name.apply(lambda x: x[:-2])
    frame4 = frame4.drop(['X', 'Y', 'Sensor', 'Name'], axis=1)
    frame4.index = index
    
    val4 = frame4.Value
    frame = pd.concat([frame3, val4], axis=1, )
    frame.columns = ['Year', 'Day', 'B3', 'B4']    # Переименуем столбцы
    frame = frame.dropna()
    
    frame['NDVI'] = (frame.B4 - frame.B3)/(frame.B4 + frame.B3)

    frame['DayNum'] = frame.apply(day_number, axis=1)
    frame = frame.sort(columns = ['DayNum'])
    return frame

Чтение данных MODIS

Аналогично чтению данных Landsat определим функцию, которая читает данные MODIS и возвращает результат в виде фрейма даных.


In [11]:
def get_modis_ndvi_frame(datafile):
    
    def day_number(row):
        d = datetime.date(int(row['Year']), 1, 1)
        d = d.toordinal()
        return d + row['Day'] - 1
    
    frame = pd.io.parsers.read_csv(datafile, sep=';', na_values='null')
    frame = frame.drop(['X', 'Y', 'Sensor', 'Name'], axis=1)
    
    frame.columns = ['Year', 'Day', 'NDVI']    # Переименуем столбцы
    frame = frame.dropna()

    frame['DayNum'] = frame.apply(day_number, axis=1)
    frame = frame.sort(columns = ['DayNum'])
    return frame

Проверка работы функции (отключена, чтобы не производить лишней работы при импортах):


In [12]:
# datafile = 'data/MODIS/point_5_NDVI.csv'
# frame = get_modis_ndvi_frame(datafile)
# plot(frame.DayNum, frame.NDVI)