Логические данные

Данные, относящиеся к логическому типу logical, имеют всего два значения: TRUE ("истина") и FALSE ("ложь").

Чаще всего, данные этого типа возникают при сравнении переменных:

x = 2; y = 5
z = x > y      # x больше y? 
z              # напечатаем полученное логическое значение
class(z)

Основные логические операции: & (И), | (ИЛИ) и "!" (НЕ):

a = TRUE; b = FALSE # Можно короче: a = T; b = F
a & b          # u И v 
a | b          # u ИЛИ v 
!a             # НЕ-u

Находясь в одном выражении с числами, данные логического типа приводятся к соответствующему числовому типу данных: FALSE представляется как 0, а TRUE — как 1

1 + FALSE
1 + TRUE
class(as.integer(1) + TRUE)

Схема неявного приведения типов в R выглядит так:

logical -> числовые типы -> character

Например:

paste(FALSE, 1.0, "test")

Явное приведение типов выполняется функциями семейства as.*, а проверку можно выполнить с помощью функций is.*.

Элементарные типы данных в R не ограничиваются числовыми, строковыми и логическими. Есть еще категориальные данные (факторы), даты и другие типы. Но о них мы поговорим позже.

#R

Особые величины

Inf, -Inf — бесконечности ∞ и -∞, NaN — не-число (Not-a-Number):

-1 / 0
0 / 0
Inf / NaN
log(Inf)
Inf + NA

NA (Not Available) — пропуск в данных. Например, под числа зарезервировано место, но сами числа не заданы:

(bird_counts <- c(1, 5, NA, 8, NA))

NULL — ничто:

(bird_counts <- c(1, 5, NULL, 8, NULL))

#R

О создании космических аппаратов для проекта «Сфера»

Роскосмос заключил контракт на разработку и изготовление шести космических аппаратов (первый этап развёртывания) системы спутникового ШПД «Скиф», входящей в федеральный проект «Сфера». Защита эскизного проекта планируется до конца 2023 года. Комплектация и начало изготовления спутников — в следующем году. Их попарный запуск — в 2026 и 2027 годах.

Запущенный в прошлом году демонстрационный аппарат «Скиф-Д» завершает этап лётно-экспериментальной отработки, успешно защитив орбитально-частотный ресурс системы.

Кроме того, Роскосмос ведет работу по включению в «Сферу» многоспутниковой системы обзорной съёмки «Грифон» с привлечением внешних инвестиций. Её 136 малых космических аппаратов обеспечат получение данных с территории России каждые 30 часов.

Также в 2023 году Госкорпорацией подписан контракт и начаты работы по формированию облика низкоорбитальной системы прямого спутникового доступа (Direct-to-Cell). После чего начнётся эскизное проектирование системы и создание двух экспериментальных аппаратов. В проекте сделан упор на привлечение частных и малых инновационных компаний и институтов.

Подробности — в материале на сайте!

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

Питер Атвуд (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