Работа с пространственными данными в R

Начало работы

* R Онлайн
* Установка R и RStudio на локальный компьютер
* Первый скрипт

Основы языка

* Арифметические операции и имена
* Типы данных: числовые, строковые, логические
* Особые величины
* Структуры данных (обзор):
векторы — часть 1, часть 2, часть 3
матрицы и массивы
списки
таблицы
* Управление ходом вычислений
условные операторы
циклы
* Базовая графика
* Функции, функциональное программирование

Пространственные данные

* Растровые данные
* Чтение растровых данных. Метаданные
* Создание растровых данных
* Арифметические и логические действия с растровыми данными
* Многослойные растровые данные
* Охват (экстент)
* Изменение пространственного разрешения
* Замена системы координат
* Слияние и изменение размеров
* Применение функций к слоям данных
* Применение функций к слоям данных-2: app, lapp, tapp и xapp
* Категоризация значений: classify
* Доступ к ячейкам данных
* Фильтрация
* Отображение растровых данных
* SpatRasterDataset и SpatRasterCollection

* Векторные данные
* Создание векторных объектов: точки, линии и полигоны
* Чтение и запись данных
* Преобразование координат
* Работа с элементами SpatVector
* Логические операции
* Пространственные отношения между объектами

* Взаимодействие растров и векторов
* Обрезка и маскирование растров по векторным данным
* Извлечение значений растров

...

Основные пакеты R для работы с пространственными данными
Справка по функциям пакета terra
Ссылки по работе с пространственными данными в R

#R

Спутниковую группировку наблюдения Земли “Грифон” включат в проект “Сфера”

ТАСС сообщил, что проект "Грифон" включат в программу создания многоспутниковой группировки "Сфера".

Проект будет реализовываться на базе специализированного отраслевого оператора с привлечением внешних инвестиций, сообщили в Роскосмосе.

Программа "Сфера" предполагает создание на орбите пяти спутниковых группировок, предоставляющих телекоммуникационные услуги, и еще пяти группировок дистанционного зондирования Земли. Одной из таких группировок станет “Грифон”.

"Грифон" — проект спутниковой группировки наблюдения Земли, которая будет включать 136 космических аппаратов на базе платформы CubeSat. Она позволит получать данные каждые 30 часов с территории России и не реже 38 часов по всему миру, — сообщает ТАСС. Вероятно, речь идет о том, что территория России будет полностью покрыта данными спутниковой съемки за 30 часов.

Ранее, Роскосмос сообщал, что инициатором создания группировки “Грифон” является головной институт госкорпорации — ЦНИИмаш, “в партнерстве с ведущими ВУЗами страны”.

ComNews сообщает, что группировка "Грифон" будет состоять из космических аппаратов, ведущих обзорную съемку. Она должна дополнять группировку сверхвысокодетального мониторинга "Беркут".

Источник рисунка

#россия

Числовые данные

В R три числовых типа данных:

* numeric — вещественные двойной точности (= double);
* integer — целые;
* complex — комплексные.

Числа по умолчанию создаются вещественными

a <- 1 # или a <- 1.0
class(a)
str(a)

Функция class отвечает на вопрос, к какому типу данных (классу) относится объект. Функция str возвращает структуру объекта: тип данных, к которому относится объект и его значение.

Чтобы создать целое число, его запись нужно завершить буквой L.

b <- 10L
class(b)

Целое число можно преобразовать из вещественного с помощью функции as.integer:

b <- as.integer(a)

Для проверки принадлежности к тому или иному типу данных используют функции is.*

is.integer(b)
is.integer(a)
is.numeric(a)

Комплексное число получается добавлением к целому или вещественному мнимой части, обозначенной символом i:

c <- sqrt(-1)     # так не пойдёт!
str(c)
c <- sqrt(-1+0i)  # правильный способ
class(c)

Неявное приведение числовых типов данных выполняется по схеме:

integer -> numeric -> complex

Например, складывая целое и вещественное числа, мы получим в результате вещественное число, а добавляя к сумме комплексное число получим комплексный же результат.

Явное приведение типов данных выполняется при помощи функций as.*:

as.numeric(1.5+1i) # при потере мнимой части выдается предупреждение

#R

Символьные (строковые) данные

Данные символьного типа character представляют собой строки, заключенные в одинарные или двойные кавычки

s1 <- "Text string."
s2 <- 'Another string.'
s12 <- paste(s1,s2)
str(s12)

Соединение строк (конкатенация) выполняется функцией paste.

Кавычки одного типа могут находиться внутри кавычек другого типа:

"First 'attempt'"
'Second "attempt"'

Если тип кавычек совпадает, то внутренние кавычки нужно экранировать при помощи '\':

"Third \"attempt\""

Выделить подстроку позволяет функция substr:

substr("abcdefg",3,5) # выделить подстроку с 3-го по 5-й элемент

Функции sub и gsub заменяют часть строки, соответствующую заданному образцу. sub заменяет только первое вхождение, а gsub выполняет глобальную замену подстроки

# sub("что_заменять", "на_что_заменить", "где_искать")
sub("a","A","rama")  # sub = substitute
gsub("a","A","rama") # gsub = global substitute
sub("a","A",c("mama","rama","banana"))

#R

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

Данные, относящиеся к логическому типу 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