Визуализация лесных пожаров в Канаде

Питер Атвуд (Peter Atwood) создал впечатляющую анимацию рекордных лесных пожаров в Канаде с мая по октябрь 2023 года, с использованием модели NASA GEOS-FP и данных о пожарах FIRMS. За этот период:

* Наблюдалось более 6 500 пожаров
* Выгорело 18,5 млн. га леса, что составляет около 5% площади лесов Канады
* Эвакуировано около 200 тысяч человек

У Питера Атвуда есть сайт и X с массой интересных карт и видео.

А тем временем, пожароопасный сезон в Канаде продолжается.

#пожары

Сегодня, с утра, познакомимся со структурами данных в R. Подробнее — будет в четверг и в пятницу.

Во вторую половину дня будет про углеродные кредиты и про то, как оценивать биомассу леса, чтобы такие кредиты получить. На примере Verra.

Структуры данных: обзор

Векторы

* одномерные
* содержат данные одного типа

eg <- c("This", "is", "a", "character", "vector")
(eg2 <- c(1, 2, 3, 4, 5, 6, 7))
is.vector(eg)

Матрицы

* двумерные числовые таблицы, со строками и столбцами

# `rnorm` - генератор случайных чисел
# с нормальным распределением
data <- rnorm(8)
mat <- matrix(data, nrow = 2, ncol = 4)
is.vector(mat)
is.matrix(mat)

Массивы

* многомерные числовые данные
* измерения задаются аргументом dim

arr <- array(rnorm(8), dim = c(2, 2, 2))

Списки

* контейнеры для любых комбинаций объектов и структур данных

lst <- list(eg, mat, df, arr, list(eg, eg2))
str(lst) 
class(lst[1])
class(lst[[1]])

Таблицы (Data Frames)

* двумерные таблицы, как в Excel
* каждая колонка может иметь свой тип данных

df <- data.frame(char = c("a", "b"), num = rnorm(2))
str(df)

#R

Углеродные кредиты и компенсации

Спутниковые методы оценки биомассы леса имеют большие погрешности, нередко достигающие 40–50% от средней биомассы леса по району интереса. В то же время, существуют рынки углеродных кредитов, где (вроде бы) дают деньги за предотвращенные или удаленные из атмосферы выбросы. Какими методами оценки биомассы пользуются на этих рынках? Ведь сомнительно, чтобы деньги давали за оценки с точностью 50%…

Сначала — несколько слов об углеродных рынках, кредитах и компенсациях. Естественно, это будет наш взгляд на тему, который не претендует на абсолютную истинность.

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

Введение нормативов на выбросы CO2 заставило предприятия искать пути снижения этих выбросов. Так возникли углеродные рынки, которые превращают выбросы CO2 в товар, устанавливая на него цену. Появились углеродные кредиты (carbon credits) — выпущенные государством квоты на выбросы углекислого газа. Компания, у которой есть неизрасходованные углеродные квоты, может продать их компании, которой таких квот не хватает. Торговля углеродными кредитами происходит в рамках государственных программ торговли углеродными квотами (Emissions Trading Scheme, ETS).

Государственные программы торговли углеродными квотами действуют далеко не везде. Например, в США такая программа существует только в Калифорнии. Поэтому, помимо государственного рынка торговли квотами, появился добровольный углеродный рынок (voluntary carbon market), не имеющий государственного регулирования. Углеродные кредиты, действующие на добровольных рынках, называются "компенсациями углеродных выбросов" (carbon offset) или просто —  углеродными компенсациями.

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

Проект подается в какую-то из программ углеродной компенсации и проверяется на соответствие ее требованиям. Если все в порядке, то рассчитываются объемы CO2, которые удалось удалить из атмосферы (или не допустить их выброса) благодаря проекту. Единица углеродного кредита (компенсации) равна одной тонне сокращенного или удаленного из атмосферы CO2.

Углеродные кредиты, которые генерирует проект, помещается в реестр программы углеродной компенсации, где они дожидаются своего покупателя. Здесь уже — как повезет.

Одни из самых известных программ углеродной компенсации предлагает Verra — НКО со штаб-квартирой в Вашингтоне, работающая с 2006 г. Для каждой программы Verra есть своя методика валидации и верификации, применяемая к проектам. Методика — это довольно обширный документ, в состав которого входит и метод оценки биомассы леса. В следующий раз мы рассмотрим метод оценки надземной биомассы леса VT0005 Tool for measuring above ground live forest biomass using remote sensing, v1.0, который используется во многих методиках Verra.

#климат #лес

Метод оценки надземной биомассы леса в проектах углеродной компенсации

Продолжим разговор, начатый здесь.

Рассмотрим метод оценки надземной биомассы леса VT0005 Tool for measuring above ground live forest biomass using remote sensing, v1.0.

Метод взят из статьи: Asner, G.P., Mascaro, J., Anderson, C. et al. High-fidelity national carbon mapping for resource management and REDD+. Carbon Balance Manage 8, 7 (2013). https://doi.org/10.1186/1750-0680-8-7

Данными ДЗЗ, использующимися в ходе оценки биомассы, являются данные о высоте леса, полученные авиационным лидаром. Asner et al. показывают, что расчеты биомассы по прогнозной модели, опирающейся на подобные данные, дают погрешность около 10% в широком диапазоне экологических условий.

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

Для оценки высоты леса используются данные авиационных лидаров. Заменить авиационные лидары спутниковыми пока нельзя, так как “пятно” (footprint — след сигнала лидара на земной поверхности) авиалидаров, используемых Asner et al. имело диаметр менее метра, а у космического лидара GEDI диаметр “пятна” составляет около 30 метров. Соответственно, поплывут оценки высоты деревьев, а значит и оценки биомассы.

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

Результатом расчетов является число — надземная биомасса леса в границах района интереса. Карту надземной биомассы строить не нужно. Asner et al., кстати, пытаются такую карту построить, и она, естественно, получается не слишком точной. Но это, так сказать, бонус, и в метод VT0005 он не входит.

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

#AGB #лес

Очередной плановый выход в открытый космос с борта МКС завершился сегодня, в 4:30 МСК.

За 7 часов 41 минуту космонавты:

* отключили дополнительный теплообменник от наружных контуров теплового режима модуля “Наука”, осмотрели и сфотографировали место утечки теплоносителя;
* вынесли радиолокатор из модуля “Поиск”, соединили его с адаптером и смонтировали на пассивном устройстве фиксации УФП-2 на “Науке”.
* запустили студенческий наноспутник “Парус-МГТУ”. На спутнике отрабатывается технология развертывания солнечного паруса.

Фото: Роскосмос

Векторы

Вектор — это последовательность элементов одного типа:

# числовой вектор
x <- c(1.5, 6, 8.3, 9, 6, .6, 2e-4)
# символьный вектор
s <- c("s", "t", "r", "i", "n", "g", "another string")
# логический вектор
b <- c(TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE)

Для создания векторов служит функция c(), название которой происходит от английского concatenate (собирать).

Элементы векторов нумеруются, начиная с единицы. Чтобы выбрать элемент вектора, нужно указать его номер в квадратных скобках:

x[1]

Скаляров как таковых в R нет: обычное число представляет собой числовой вектор единичной длины:

a <- 3
a[1]
# но:
a[2]

Возможно вас интересовало, почему перед выводимыми в консоли R результатами стоит [1]? Так вот, это номер элемента вектора, с которого начинается строка вывода:

c("s", "t", "r", "i", "n", "g", "another string")

Двоеточие позволяет создать последовательность элементов с шагом 1

-5:5

Отрицательный индекс в квадратных скобках означает: выбрать все элементы, кроме указанных:

v <- c(1.1, 2.2, 3.3, 4.4, 5.5)
v[1:4]
v[-5]

Векторы в R, как массивы в С, занимают непрерывный блок памяти, поэтому вставлять или удалять элементы в них невозможно. При попытке изменить вектор x в действительности создается новый вектор, который сохраняется с именем исходного (х).

Тип данных, составляющих вектор, и его структуру, как и раньше можно определить с помощью функций class и str:

class(v)
str(v)

Двоеточие : обладает более высоким приоритетом, чем вычитание. Поэтому для создания последовательности чисел от 1 до i-1, последнее число необходимо заключить в скобки

i <- 3 
1:i-1   # Это означает (1:i) - 1, а не 1:(i-1) 
1:(i-1) # так правильно

Функция seq позволяет создавать последовательности с заданным шагом

seq(5,1,by=-.5)

и заданной длины

seq(1,10,length=6)

В seq есть и другие аргументы. Узнать о них можно, как обычно, из справки о функции.

Функция rep позволяет повторить объект заданное число раз

rep(1:3,5)

А вот пример похитрее

rep(1:3,c(5,5,5))
# то же самое:
rep(1:3,rep(5,3))

#R

Векторы
(Продолжение)

Арифметические операции над векторами выполняются поэлементно:

u <- c(1,2,3)
v <- c(4,5,6)
u+v
u*v
u/v

Скалярное произведение векторов записывается так:

u %*% v

Деление на ноль дает в результате Inf (бесконечность):

w <- c(v[1:2],0) # добавляем элемент к фрагменту вектора v
u/w

которая при последующих операциях "поглощает" все конечные значения:

u+u/w

Сложим два вектора разной длины. "Будет ошибка", — скажите вы. А вот и нет:

c(1,2,3) + c(4,5,6,7,8,9,10)

Операция будет выполнена, но R предупредит, что длины векторов-слагаемых различаются.

Сложение и другие подобные операции, требующие равной длины операндов, выполняется по правилу: элементы короткого вектора c(1,2,3) повторяются до тех пор, пока длина этого вектора не сравняется с длиной c(4,5,6,7,8,9,10), после чего выполняется заданная операция. Фактически, складываются векторы:

c(1,2,3,1,2,3,1) + c(4,5,6,7,8,9,10)

Добавление элементов в вектор осуществляется функциями c и append ('добавить'):

vec <- c('a','b')
vec <- c(vec,'c','d')
vec
values <- c('e','f','g')
vec <- append(vec, values)
vec

Удаление элементов из вектора выполняется следующим образом:

a <- sample(1:10)  # генерируем случайные целые числа от 1 до 10
a
remove <- c(3,5,7) # выберем для удаления 3-й, 5-й и 7-й элементы 
a <- a[-remove]    # удалим выбранные элементы
a

Длина вектора, то есть число его элементов, вычисляется функцией length():

length(a)

Указать последний элемент вектора можно так:

a[length(a)]

#R

Векторы
(Окончание)

В R повсеместно используется векторизация, то есть подход к программированию, когда операции выполняются над вектором в целом, а не над отдельными его элементами (скалярами).

a <- c(1,2,3); b <- c(4,5,6)
c1 <- vector() # создаем пустой вектор

Вместо того, чтобы делать так:

for (i in 1:3) {
  c1[i] <- a[i] + b[i]
}

в R поступают так:

c2 <- a + b

Для векторизации расчетов используется логическая индексация:

{r}
a <- c(6,-2,1,8,0,9)
ind_a <- a > 0
ind_a

Логический индекс (ind_a) — вектор, длиной равный исходному (a), элементы которого равны TRUE, если соответствующий элемент исходного вектора удовлетворяет логическому условию (a > 0) и FALSE — в противоположном случае.

Логическая индексация позволяет заменить связку "цикл + условный оператор". Например, чтобы выбрать положительные элементы вектора a не нужно организовывать цикл с проверкой в его теле условия a[i] > 0. Вместо этого поступают так:

a[a > 0] # или a[ind_a]

Примеры использования логических операций:

a <- c(6,-2,1,8,0,9)
a > 0 & a < 9 # логическое И
a < 2 | a > 8 # ИЛИ
# Истинно, если хотя бы один
# из элементов аргумента истинен.
any(a>0)
# Истинно, если все элементы аргумента истинны.
all(a>0)

Если данные содержат пропуски (NA), это может повлиять на результат вычислений. Проверка пропусков реализуется с помощью is.na()



# Данные с пропусками:
a <- c(6,-2,NA,1,8,0,NA,9)
# Их сумма дает:
sum(a)
# Является ли элемент пропуском в данных?
is.na(a)

У многих функций есть аргумент na.rm, управляющий предварительным удалением пропусков

# Cуммирование элементов, 
# с предварительным удалением NA
sum(a, na.rm=T)

Покажем как с помощью векторизации можно легко вычислить определенный интеграл.

Вычислим интеграл (см. рисунок), воспользовавшись методами прямоугольников и трапеций Вспомнить их можно по книге: Турчак Л. И., Плотников П. В. Основы численных методов. – М.: Физматлит, 2003. Для проверки: интеграл равен 7/3 = 2.333(3).

# границы промежутка интегрирования
a <- 1; b <- 2
# число узлов интегрирования
n <- 1000
# координаты узлов сетки
x <- seq(a,b,length.out=n)
# шаг сетки
h <- x[2]-x[1]
# значения подынтегральной функции в узлах сетки
y <- x^2

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

# Нижняя интегральная сумма
sd <- h*sum(y[-length(y)])
sd
# Верхняя интегральная сумма
su <- h*sum(y[-1])
su

Метод трапеций дает более точный результат:

(su+sd)/2

#R

Матрицы и массивы

Матрица (matrix) — это прямоугольная числовая таблица, состоящая из строк и столбцов. Создадим из элементов вектора v1 матрицу размера 2х2:

v1 <- 1:4
m1 <- matrix(v1, nrow=2, ncol=2)
m1
class(m1) 
attributes(m1)

Матрицы можно собирать из блоков, объединяя другие матрицы по строкам (rbind) или по столбцам (cbind)

m2 <- matrix(-2:1, nrow=2, ncol=2)
m3 <- rbind(m1,m2)
m3
m4 <- cbind(m1,m2)
m4

При выборе элементов матрицы указывают индексы строк и столбцов, содержащих нужные элементы:

# Элемент 4-й строки и 2-го столбца матрицы m3
m3[4,2]

Размерность матрицы (и вектора) можно задать с помощью dim:

A <- 1:9
# Фрмируем из вектора матрицу размерности 3 х 3
dim(A) <- c(3,3)
A

Транспонирование матрицы выполняется функцией t():

t(A)

Многомерные массивы в R создаются при помощи функции array. Размерность массива задаётся атрибутом dim:

arr <- array(1:24, dim=c(2,4,3))
arr

#R

На Google Earth Engine появились…

🌲Глобальные карты плотности надземной биомассы Biomass CCI:

ESA CCI Global Forest Above Ground Biomass

Напомним, что данные Biomass CCI текущей 4-й версии имеют пространственное разрешение: 100 м (на экваторе) и содержат плотность биомассы за 2010, 2017, 2018, 2019 и 2020 гг., измеряемую в Мг/га (тонн/га).

🌡Ежедневные приземные метеорологические данные для сельскохозяйственных и агроэкологических исследований:

AgERA5 (ECMWF) dataset

Они основаны на почасовых данных ECMWF ERA5 на уровне поверхности, и содержат показатели:

* Давление водяного пара
* Доля продолжительности выпадения осадков в жидкой фазе
* Доля продолжительности выпадения осадков в твердой фазе
* Общая облачность
* Объем выпавших осадков
* Относительная влажность воздуха на высоте 2 м
* Поток солнечной радиации
* Скорость ветра на высоте 10 м (над поверхностью)
* Температура воздуха на высоте 2 м
* Толщина снежного покрова
* Точка росы
* Эквивалент жидкой воды в снежном покрове

Пространственное разрешение: 0,1° (9600 м).
Временное покрытие: с 1979 года по настоящее время.

#AGB #лес #данные #GEE #погода

Списки

Список (list) похож на вектор, но может содержать объекты разных типов. Создается список функцией list:

lst <- list(n=1:3, s="test", b=c(T,F,F,F,T))
# или
n <- 1:3
s <- "test"
b <- c(T,F,F,F,T)
lst <- list(n,s,b) # состоит из копий векторов n,s,b

Список lst состоит из трех элементов, каждый из которых является вектором: числовым, символьным и логическим соответственно.

Элементы списков, подобно элементам векторов и матриц, могут иметь имена. Узнать или задать эти имена позволяет функция names.

Присвоим элементам списка lst имена num, cha и log:

names(lst) <- c("num","cha","log")

Выбрать элементы списка можно одним из следующих способов:

1. с помощью квадратных скобок [,]. Полученный объект будет в свою очередь являться списком, даже если длина его окажется единичной

lst[1]  # или lst["num"]
class(lst[1])

2. с помощью двойных квадратных скобок [[,]]. Возвращенный объект будет того же типа, каким он был до включения в список

lst[[1]] # или lst[["num"]]
class(lst[[1]])

3. использовать имена элементов списка и знак доллара $:

lst$num
class(lst$num)

Удалить имена элементов можно так: names(lst) <- NULL. Мы пока не будем этого делать.

Добавить элемент в список можно так же, как и в вектор — при помощи функций c() или append():

v <- 5:10             # создадим числовой вектор
l1 <- c(lst,list(v))  # преобразуем его и добавим в список
str(l1)

Если вектор v не преобразовать при помощи list(), то к списку будет добавлен не один элемент-вектор, а шесть элементов-чисел:

l2 <- c(lst,v)
str(l2)

Функции c и append позволяют объединять не только векторы, но и списки:

l3 <- append(list(1:2),list(3:4))

Удаляются элементы списка так же, как элементы вектора:

remove <- c(2,3)      # удалить 2-й и 3-й элементы 
lst <- lst[-remove]

Определим количество элементов в списке lst:

length(lst)

Обратите внимание, что функция length возвращает нам число элементов списка. В нашем случае такой элемент один. Другое дело, что сам этот элемент является числовым вектором. Чтобы найти длину этого последнего нужно применить двойные квадратные скобки

length(lst[[1]])

или преобразовать список в вектор при помощи unlist():

vec_from_lst <- unlist(lst)
class(vec_from_lst)
length(vec_from_lst)

Заметили обращение к элементам списка через знак доллара (lst$num)? Поскольку в именах переменных R можно использовать точку, то обращение к методам объекта или к элементам данных должно осуществляться при помощи другого символа. Этим символом и является $.

Соберем вместе особенности синтаксиса R:

1. Присваивание обозначается стрелкой <- или знаком равенства =.
2. В именах переменных можно использовать точку '.'.
3. Обращение к методам объекта или элементам структуры данных осуществляется через $.

#R