Ранее рассматривался вопрос о том, что данные ДЗЗ можно классифицировать не только по одномоментному снимку, но и по длинной временной серии. В частности, растительность можно попробовать классифицировать по многолетней серии NDVI: предполагается, что разные классы растительности обладают различными формами кривых NDVI. Эти формы могут выступать в качестве признаков классификации. Форму кривой можно описать, воспользовавашись преобразованием Фурье, которое представляет сигнал в виде числовых коэффциентов. Поэтому полученные коэффициенты могут служить признаками классификации. Более подробно об этом сказано в статье, приведенной выше.
Основная цель данной статьи -- экспериментальная проверка пригодности описанного подхода на реальных данных. Реализация алгоритма выполнена в GRASS GIS с использованием языка Python.
Статья состоит из нескольких обособленных частей:
Использование GRASS в качестве вычислительной платформы остается за рамками данной статьи. Большую часть необходимой информации можно почерпнуть в документации.
В этом разделе приводится код инициализации и некоторые примеры, которые необходимы для понимания программного кода, который будет использован в основной части статьи.
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
Создадим объект, который будет отвечать за связь с GRASS GIS и инициализировать GRASS:
In [2]:
class GRASS():
def __init__(self, gisbase, gisdbase, location, mapset='PERMANENT'):
self.gisbase = os.environ['GISBASE'] = gisbase
self.gisdbase = gisdbase
self.location = location
self.mapset = mapset
# Пока не обновлены пути, импортировать GRASS не удастся
sys.path.append(os.path.join(os.environ['GISBASE'], "etc", "python"))
import grass.script as grass
import grass.script.array as garray
import grass.script.setup as gsetup
self.grass = grass
self.garray = garray
self.gsetup = gsetup
# Активируем GRASS
self.gsetup.init(self.gisbase,
self.gisdbase, self.location, self.mapset)
def get_region_info(self):
gregion = self.grass.region()
return gregion
def get_NDVI_names(self):
"""Возвращает список растров, хранящих значения NDVI,
отсортированный в алфавитном порядке.
"""
raster_list = self.grass.read_command(
'g.list', type='rast',
pattern='MOD13Q1.*250m_16_days_NDVI'
)
raster_list = raster_list.split()
raster_list.sort()
return raster_list
def ndvi_as_array(self):
"""Возвращает растры NDVI в виде трехмерного np.array
"""
names = self.get_NDVI_names()
rast_count = len(names)
region = self.get_region_info()
row_count = region['rows']
col_count = region['cols']
data = np.zeros((row_count, col_count, rast_count))
map = self.garray.array()
for i, name in enumerate(names):
map.read(name)
data[:, :, i] = map.copy()
return data
Перед началом работы нужно обязательно инициализировать GRASS, для этого нужно передать путь к каталогу, в который установлен GRASS, путь к каталогу БД GRASS, а также название MAPSET. Делается это следующим образом:
(Здесь и далее комманды закомментированы, чтобы не выполнялись при импорте, когда документ будет использоваться в качестве скрипта Результат выполнения команд не удалялся).
In [3]:
# grass = GRASS('/opt/grass70-svn', '/home/dima/GIS/GRASSDATA', 'ugra')
После этого будут установлены необходимые переменные и созданным объектом можно будет пользоваться: будут доступны как обычные команды GRASS, так методы, определенные в объекте:
In [4]:
# Команды GRASS
# print grass.grass.read_command('g.region', flags='p')
# # "Самодеятельные" методы
# raster_list = grass.get_NDVI_names()
# print raster_list[:5]
Информация по текущему региону:
In [5]:
# region = grass.get_region_info()
# print region
В нашем случае наиболее важный метод -- ndvi_as_array, при помощи которого можно получить все растры NDVI в виде трехмерного массива (стека растров). В первом измерении хранятся строки растров, во втором -- столбцы, сами растры хранятся в третьем. Считываются данные из текущего региона, с учетом его границ и разрешения. Также учитывается маска.
In [6]:
# data = grass.ndvi_as_array() # Закомментировано, чтобы не вызывать при импорте
Далее полученным массивом можно пользоваться, как обычным numpy.array. Для примера отобразим динамику NDVI для некоторых пикселей:
In [7]:
# fig = plt.figure(figsize=(12,6), dpi=100)
# plt.plot(data[10, 20, :], '-o') # строка 10, колонка 20
# plt.plot(data[40, 40, :], '-o')
# plt.plot(data[400, 140, :], '-o')
# plt.show()