Здесь содержится код вспомогательных функций и пояснения к их работе:
Проблемы. Здесь нехорошо то, что из этого файла будут импортироваться код функций. При этом весь код, поясняющий работу этих функций будет выполняться в момент импорта. Нужно придумать способ спрятать поясняющий код, чтобы при импорте он не вызывался.
In [1]:
import datetime
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
In [2]:
url = 'data/Landsat/point_1_band3.csv'
frame3 = pd.io.parsers.read_csv(url, sep=';', na_values='null')
frame3.head()
Out[2]:
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]:
Повторим процедуру с 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]:
Расчитаем 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]:
In [6]:
frame['NDVI'] = (frame.B4 - frame.B3)/(frame.B4 + frame.B3)
frame.head()
Out[6]:
Отортируем данные по году и дню, построим график 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]:
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
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)