Сегментация данных ДЗЗ по временным сериям NDVI на базе преобразования Фурье в GRASS GIS

Ранее рассматривался вопрос о том, что данные ДЗЗ можно классифицировать не только по одномоментному снимку, но и по длинной временной серии. В частности, растительность можно попробовать классифицировать по многолетней серии NDVI: предполагается, что разные классы растительности обладают различными формами кривых NDVI. Эти формы могут выступать в качестве признаков классификации. Форму кривой можно описать, воспользовавашись преобразованием Фурье, которое представляет сигнал в виде числовых коэффциентов. Поэтому полученные коэффициенты могут служить признаками классификации. Более подробно об этом сказано в статье, приведенной выше.

Основная цель данной статьи -- экспериментальная проверка пригодности описанного подхода на реальных данных. Реализация алгоритма выполнена в GRASS GIS с использованием языка Python.

Статья состоит из нескольких обособленных частей:

  1. Описание и исходный код класса, взаимодействующего с GRASS, и который можно использовать в скриптах для автоматизации работы.
  2. Функции анализа.
    1. Выбор коэффициентов разложения.

Объект-обертка вокруг GRASS GIS

Использование 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 (для 7-й версии GRASS)

Создадим объект, который будет отвечать за связь с 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]


projection: 1 (UTM)
zone:       42
datum:      wgs84
ellipsoid:  wgs84
north:      6924647.3356839
south:      6666561.94859838
west:       347504.03290255
east:       617968.50593856
nsres:      502.11164803
ewres:      501.78937483
rows:       514
cols:       539
cells:      277046

['MOD13Q1.A2011001.250m_16_days_NDVI', 'MOD13Q1.A2011017.250m_16_days_NDVI', 'MOD13Q1.A2011033.250m_16_days_NDVI', 'MOD13Q1.A2011049.250m_16_days_NDVI', 'MOD13Q1.A2011065.250m_16_days_NDVI']

Информация по текущему региону:


In [5]:
# region = grass.get_region_info()
# print region


{'rows': 514, 'e': 617968.50593856, 'cells': 277046, 'cols': 539, 'n': 6924647.3356839, 's': 6666561.94859838, 'w': 347504.03290255, 'ewres': 501.78937483, 'nsres': 502.11164803}

В нашем случае наиболее важный метод -- 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()