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

В 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

Таблицы

Таблицы (data frames), как и матрицы, хранят двумерные данные. Но, в отличие от матриц, в разных колонках таблицы могут находиться данные разных типов. По внутреннему устройству, таблица — это список колонок равной длины. Если вы видели таблицу в Excel, то вам уже знакомы таблицы в R.

Создаются таблицы с помощью функции data.frame:

df <- data.frame(a=1:5, b=6:10, c=letters[1:5])
str(df)

Выбор элементов таблицы рассмотрим на примерах:

df[1,]            # 1-я строка
df[,2]            # 2-ая колонка
df[4,3]           # элемент 4-й строки 3-й колонки
df[4,c(1,3)]      # элемент 4-й строки из 1-ой и 3-й колонок
df$c              # колонка с именем "c" 
df[,"b"]          # колонка с именем "b" (2-я)

Часто для выбора нужного фрагмента таблицы удобно использовать логические условия:

# строка, в которой элемент 3-й колонки равен "d"
df[df$c == 'd', ] # 4-я строка
# строки, в которых элементы 1-й колонки больше или равны 3
df[df$a >= 3,]    # строки 3-5

Таблицы можно объединять по строкам или по столбцам:

df2 <- data.frame(a=6:10, b=11:15, c=letters[6:10])
rbind(df,df2)
cbind(df,df2)

Еще один способ объединения таблиц дает в результате список:

l3 <- c(df,df2)
str(l3)

Для объединения таблиц, имеющих общие колонки, используется функция merge. Как правило, с ее помощью объединяют таблицы с общими ключами:

first <- data.frame(id=1:5,
                    COL=LETTERS[1:5],
                    col=letters[1:5])
secnd <- data.frame(id=1:5, col=letters[1:5])
# объединяем таблицы по колонке id
total <- merge(first,secnd,by="id")
# объединяем таблицы по колонкам id и col
total <- merge(first,secnd,by=c("id","col")) 

Размеры таблицы определяются следующим образом:

# длиной таблицы является длина списка, 
# состоящего из колонок
length(df)
ncol(df)   # число колонок
nrow(df)   # число строк

Для работы с именами колонок, кроме names(), существует colnames(). Функция rownames задает имена строк таблицы:

rownames(df) <- LETTERS[1:nrow(df)]

Теперь выбрать элемент таблицы df можно так:

df["A","b"] # df[1,2]

#R

По итогам совещания президента РФ Владимира Путина в РКК "Энергия" РИА Новости сообщает, что предприятия "Роскосмоса" занимаются организацией серийного производства спутников, благодаря чему к 2036 году предполагается развернуть на различных орбитах более двух тысяч новых космических аппаратов разных типов. Для сравнения, на июнь этого года орбитальная группировка России насчитывала 225 спутников.

Ранее гендиректор "Роскомоса" Юрий Борисов говорил, что минимальная необходимая спутниковая группировка России к 2030 году должна состоять из 1000–2000 аппаратов. Для этого уже к 2025 году нужно производить по 250 спутников в год, а к 2030 — один аппарат в день.

До 2036 года далеко, а вот 2025 год уже близко. Скоро все и узнаем.

#россия

Обожаю такие заголовки

https://www.riatomsk.ru/article/20231024/tusur-bespilotniki-taksaciya-lesa/

— нелепо, смешно и с арифметическими ошибками.

А говорит эксперт примерно следующее:

В лесном массиве, площадью один квадратный километр, бригада таксологов будет ходить со скоростью 4 километра в час — это три рабочих дня, и один день нужен на обработку данных. Для беспилотника, со средней скоростью полета 30 километров в час, понадобится 3 часа для облета — это один рабочий день, и еще один день на обработку. Кроме того, бригада таксологов составляет 4–5 человек, а для беспилотника потребуются только двое. В итоге, применение беспилотников позволит в 4 раза сократить затраты человеко-часов на проведение таксации леса.

#лес

Заседание секции “Наблюдения и исследования Земли из космоса” Совета РАН по космосу

25 октября 2023 года состоялось заседание секции “Наблюдения и исследования Земли из космоса” Совета РАН по космосу.

На заседании:

* заслушали и обсудили доклад директора ФГБУ НИЦ “Планета” Тасенко Сергея Викторовича на тему: “Национальная спутниковая группировка для обеспечения решения задач гидрометеорологии и мониторинга окружающей среды”.
* обсудили и утвердили план работы, регламент работы и тематики проблем секции.

Кроме упомянутой, в Совете РАН по космосу есть и другие секции:

* Космическая биология и физиология
* Космическое материаловедение
* Физика Солнца
* Космическая геодинамика
* Солнечно-земные связи

а еще экспертные комиссии:

* Фундаментальные проблемы космических перелетов
* Проблемы научно-образовательных спутников
* Фундаментальные проблемы физики микрогравитации

и экспертная рабочая группа по космическим угрозам.

Все они что-то делают. Но что именно, можно узнать лишь по косвенным признакам — сайт Совета РАН не обновлялся с декабря 2020 года. А хотелось бы знать подробности. Координация исследований в области наблюдения Земли из космоса — это, знаете ли, важно.

Данные ICESat-2 для оценки высоты леса

Лидар ATLAS на борту спутника ICESat-2 — лидар с подсчетом фотонов (photon counting LIDAR). Он испускает лазерные импульсы с высокой частотой и детектирует отдельные фотоны, отраженные объектами на земной поверхности. Измерение времени полета фотона, совместно с информацией о положении и ориентации прибора, позволяет рассчитать координаты отражающего объекта земной поверхности, включая его высоту.

Пространственное разрешение прибора ATLAS, то есть размер “пятна” (footprint) лазерного луча на поверхности Земли составляет около 14 м. Частота генерации лазерных импульсов подобрана так, чтобы обеспечить смещение “пятна” каждого следующего импульса на расстояние около 70 см по мере пролета спутника. В результате происходит непрерывное получение информации вдоль трека спутника на земной поверхности для полосы шириной 14 м. Расстояние между двумя соседними треками измерений составляет около 3 км на экваторе.

Система ATLAS/ICESat-2, так же как и ее предшественник GLAS/ICESat, разработана в основном для мониторинга ледового покрова, но может быть использована и для оценки высоты леса на больших территориях, в том числе для бореальных лесов (Neuenschwander et al., 2020a). Результаты применения данных ATLAS/ICESat-2 для оценки высоты лесов России приведены в (Барталев и др., 2022) и (Шинкаренко и др., 2023).

Для оценки высоты леса используется информационный продукт ATL08, построенный по данным ATLAS/ICESat-2 и содержащий информацию о значениях высоты растительности. В указанном продукте полученное лидаром облако точек разделяется на четыре класса: подстилающая поверхность, листва/кроны (canopy), верхний полог листвы/крон (top-of-canopy) и шум. Далее, непрерывный трек, соответствующий полосе земной поверхности шириной 14 м, разбивается на сегменты длиной 100 м. Каждый стометровый сегмент состоит из пяти последовательных сегментов, длиной 20 м, взятых из продукта ATL03. Пространственное расположение таких сегментов характеризуется координатами их центров. Таким образом, данные ATL08 представлены в виде сегментов шириной 14 м и длиной 20 м (Neuenschwander et al., 2022), что позволяет обеспечить достаточно высокую пространственную детальность измерений высоты.

Элемент набора данных ATL08 соответствует одному сегменту и содержит параметры, характеризующие распределение относительной высоты растительности на данном участке: значения минимума, максимума, среднего, стандартного отклонения и различные перцентили высот растительности, а также количество измерений и показатель степени их неопределенности (uncertainty). Высота отразивших фотоны элементов растительности (кроны/листва) оценивается относительно пространственно-интерполированных точек, отнесенных к классу подстилающей поверхности. В качестве параметра, характеризующего высоту леса, используется величина h_canopy — 98 перцентиль высот растительности.


Барталев С.А., Богодухов М.А., Жарко В.О., Сидоренков В.М. Исследование возможностей использования данных ICESat-2 для оценки высоты лесов России. Современные проблемы дистанционного зондирования Земли из космоса. 2022. Т. 19. № 4. С. 195-206. https://doi.org/10.21046/2070-7401-2022-19-4-195-206

Шинкаренко С.С., Барталев С.А., Богодухов М.А., Жарко В.О. Оценка возможностей определения высоты и проективного покрытия защитных лесных насаждений по данным ICESat-2. Современные проблемы дистанционного зондирования Земли из космоса. 2023. Т. 20. №2. С. 174–183. http://dx.doi.org/10.21046/2070-7401-2023-20-2-174-183

Neuenschwander, A., Guenther, E., White, J. C., Duncanson, L., & Montesano, P. (2020). Validation of ICESat-2 terrain and canopy heights in boreal forests. Remote Sensing of Environment, 251, 112110. https://doi.org/10.1016/j.rse.2020.112110

Neuenschwander A., Pitts K., Jelley B. et al. Ice, Cloud, and Land Elevation Satellite 2 (ICESat-2): Algorithm Theoretical Basis Document (ATBD) for Land-Vegetation Along-Track Products (ATL08). Version 6. 2022. 144 p. https://doi.org/10.5067/8ANPSL1NN7YS

#лидар #лес