По сути это касается не только леса, это могут быть водные объекты, сельхозугодия, городская застройка. Вопрос в том, как вычислить любой из этих показателей на определенной единице площади, имея лишь космический снимок.
Прежде всего, скачиваем нужный по дате и территориальному охвату космический снимок с сайта геологической службы США, предварительно зарегистрировавшись на нем.
Войдя в аккаунт, в панели слева выбираем интересующий временной период, например, снимки за последние полгода, или с мая по сентябрь 1990 года. В зависимости от периода выбираем соответствующую версию снимков. Нам нужны летние снимки 2020 года.
Для справки – периоды работы спутников Landsat:
Landsat 1: 23.07.1972 – 06.01.1978
Landsat 2: 22.01.1975 – 22.01.1981
Landsat 3: 05.03.1978 – 31.03.1983
Landsat 4: 16.07.1982 – 1993
Landsat 5: 01.03.1984 – 21.12.2012
Landsat 6: на целевую орбиту не выведен
Landsat 7: 15.04.1999 – н.в.
Landsat 8: 11.02.2013 – н.в.
Скачанная папка со снимком содержит 12 растровых изображений – это каналы, и два текстовых документа. Для дальнейших действий нужен красный (RED) и инфракрасный (Nir) каналы. Если снимок восьмой версии Landsat, нужны четвертый и пятый каналы, если седьмой – третий и четвертый.
Далее вычисляем индекс NDVI (вегетационный индекс, отражающий количество фотосинтетически активной биомассы, проще говоря, зеленой растительности). Во вкладке «Растр» открываем «Калькулятор растров» и высчитываем NDVI по формуле NDVI = (Nir – RED)/(Nir + RED). Смотрим на скрин ниже и обращаем внимание на значения B4 (RED) и B5 (Nir) в имени канала.
Результат должен быть таким.
NDVI варьируется от -1 (отсутствие растительности) до 1 (густая растительность). В нашем случае его диапазон от -0,47 до 0,93.
Схема взята из этой статьи. Статья будет полезна тем, кому нужно более детально разобраться в индексах и их формулах. А мы двигаемся дальше, и для отображения лесной растительности нужно выделить значение от 0,7 до 0,93. Стандартный калькулятор растров в версии 3.10.6 почему-то не желает высчитывать столь дробное значение, поэтому мы обратились к калькулятору модуля Semi-Automatic Classification Plugin (SCP).
Значение 0,9995 обычно рассчитывается как 1, машина в этот раз так рассчитала. Все, что имеет значение 0,9995 – это та растительность, которую мы выделили в NDVI-диапазоне 0,7 – 1. Пока оставляем все как есть и переключаемся на рельеф. Скачиваем с геологической службы США нужные SRTM-снимки.
Добавляем снимки в проект и объединяем: «Растр» > «Прочее» > «Результат/Merge». Результат объединения обрезаем (все в той же вкладке «Растр») и, если необходимо, задаем нужную проекцию.
Наш слой cut_merge имеет минимальное значение -5, а максимальное – 1525. Надо отметить, что в самом SRTM есть незначительная погрешность: наивысшая точка горного Крыма – 1545 метров. Вероятно, если бы снимок был более высокого разрешения, этой неточности удалось бы избежать.
Возвращаемся к лесу. Мы имеем два растровых изображения-результата: первый отражает всю покрытую лесом территорию, второй – всю территорию выше 149 метров. Теперь надо объединить их. Открываем калькулятор растров и задаем следующее выражение.
Если стандартный калькулятор растров начнет капризничать, используйте калькулятор из SCP, о нем говорилось выше.
Работа с растрами закончена. Данный результат стандартными средствами преобразовываем в вектор: «Растр» > «Преобразование» > «Создание полигонов (растр в вектор)».
В получившемся векторном слое открываем таблицу атрибутов и все, что имеет значение 0, выделяем и удаляем. Сохраняем результат. Оставшиеся полигоны со значением 1 – это собственно векторный слой леса на высоте 149+.
Далее в пределах этого векторного леса мы попросту создаем регулярные точки на определенном небольшом расстоянии. Мы выбрали 30 метров. Точек будет очень много, поэтому точечный shape будет долго создаваться. И создастся он в пределах прямоугольника, охватывающего слой. Поэтому необходимо применить следующий инструмент: «Вектор» > «Выборка» > «Выделение по районам».
По завершении алгоритма «Выделение по районам» сохраняем выделение в отдельный слой и добавляем в проект.
Следующий шаг – создание сетки: «Вектор» > «Выборка» > «Create grid (Создание сетки)». Выбираем форму сетки (квадраты, ромбы, шестиугольники), задаем их размер, например, 4000х4000 метров.
Рассчитываем количество точек в каждой ячейке сетки: «Вектор» > «Анализ» > «Подсчитать точки в полигоне». В результате появится еще одна аналогичная сетка, но в ее таблице атрибутов будет новое поле NUMPOINTS с количеством точек в ячейках.
И наконец работа со стилями. В свойствах слоя переходим во вкладку «Стиль» и начинаем творить. Выбираем градуированный знак, поле NUMPOINTS и играемся с градиентами и классификациями.
Результат. Один из множества вариантов. После этого дело только за макетом.
