Анимация площадки строительства АЭС Аккую (Турция: 36.1444 с.ш., 33.5411 в.д.) по снимкам Sentinel-2 2016–2023 гг. (естественные цвета; снимки 2016–2018 гг. — Level 1C, 2019–2023 гг. — Level 2A).

#снимки

Сегодня мы продолжим рассказ о применении языка R для работы с растровыми данными.

Арифметические и логические действия с растровыми данными

Рассмотрим арифметические и логические действия с растровыми данными:

r <- rast(matrix(1:16, ncol = 4, nrow = 4))

r+1
r^2
r*0.3

r > 5

s <- r+1
s[s > 5] <- NA
s[is.na(s)] <- 5

Каждый раз действие применяется ко всем пикселам растра. В результате получится растр, геометрически равный исходному, но с изменившимися значениями пикселей.

Покажем равенство геометрий исходного и одного из получившихся растров:

compareGeom(r, r+1)

Чтобы выполнять арифметические и логические действия с несколькими растрами, эти растры должны иметь одинаковую геометрию, то есть одинаковое число строк и столбцов, одинаковые размеры пикселей (res), и относится к одной и той же системе координат (Coordinate Reference System, CRS).

Возьмем изображения, полученные 4-м и 5-м каналами прибора Landsat 8 OLI

red <- rast("LC08_044034_20170614_B4.tif")
nir <- rast("LC08_044034_20170614_B5.tif")

# Проверим равенство геометрий
compareGeom(red, nir)

Рассчитаем с их помощью нормализованный разностный вегетационный индекс (NDVI):

ndvi <- (nir - red) / (nir + red)

names(ndvi) <- "LC08_044034_20170614_NDVI"

plot(ndvi)

Загрузим изображение с зеленого канала, и объединим все каналы в изображение с комбинацией каналов “5-4-3” (NIR, Red, Green):

green <- rast("LC08_044034_20170614_B3.tif")

img <- rast(list(nir, red, green))

plotRGB(img, stretch="lin")

#R

Имя файла данных Landsat содержит массу полезной информации. Например, LC08_044034_20170614_B3.tif сообщает нам, что:

* LC08 — снимок Landsat-8 OLI/TIRS (The Operational Land Imager/Thermal Infrared Sensor).
* 20170614 — сделан 14 июня 2017 г.
* 044034 — расположение снимка: Worldwide Reference System (WRS) Path 44, Row 34
* B3 — номер спектрального канала (band 3, зеленый)

#R

Комбинация каналов 5-4-3 (NIR, Red, Green) применяется для изучения состояния растительного покрова, мониторинга дренажа и почвенной мозаики, а также изучения сельскохозяйственных культур. Растительность в этой комбинации имеет оттенки красного, городская застройка — зелено-голубые цвета, почва — от темно до светло коричневого или серого, лед, снег и облака — белые или светло голубые. Насыщенные оттенки красного являются индикаторами здоровой и/или широколиственной растительности. Более светлые оттенки характеризуют травянистую растительность или редколесья/кустарники. Хвойные леса по сравнению с лиственными имеют более темно-красную или даже коричневую окраску. Подробнее — см. здесь.

#комбинация #R

Многослойные растровые данные

Загрузим тестовый снимок КА “Канопус-В-ИК” и переименуем его каналы (слои):

r <- rast("KANOPUS_20190829.tif")
names(r) <- c("blue","green","red","nir")

# class       : SpatRaster 
# dimensions  : 1440, 2052, 4  (nrow, ncol, nlyr)
# resolution  : 2.1, 2.1  (x, y)
# extent      : 624695.4, 629004.6, 6086649, 6089673  (xmin, xmax, ymin, ymax)
# coord. ref. : WGS 84 / UTM zone 44N (EPSG:32644) 
# source      : KANOPUS_20190829.tif 
# names       : blue, green, red, nir

Для выделения подмножества слоев растра используют двойные квадратные скобки [[]] или функцию subset. Аргументами служат индексы или имена интересующих слоев:

# Канал 1 (синий)
r[[1]]
# Снова канал 1 (синий)
r[["blue"]]
# Все каналы, кроме 1 (2:4)
r[[-1]]
# Каналы красный и NIR (3:4)
r[[c("red","nir")]]

На выходе всякий раз получаем новый растровый объект (SpatRaster).

subset() работает аналогично:

subset(r, c(3,2,3,1))
subset(r, 2:3)
subset(r, c("red","blue"))

Первый аргумент — исходный растр, второй — вектор заданных слоев.

Рассчитаем два спектральных индекса — NDVI (Normalized Difference Vegetation Index) и NDWI (Normalized Difference Water Index):

ndvi <- (r[[4]] - r[[3]]) / (r[[4]] + r[[3]])
ndwi <- (r["green"] - r["nir"]) / (r["green"] + r["nir"])

NDVI характеризует количество фотосинтетически активной биомассы. NDWI — это, по сути, перевернутый NDVI: водные объекты имеют высокие значения этого индекса, а растительность — низкие (отрицательные). На снимке водных объектов нет, зато есть дороги и городская застройка, получающие положительные значения NDWI.

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

r <- rast(matrix(1:16, ncol = 4, nrow = 4))
# Составим многослойный растр,
# комбинируя исходный
s <- rast(list(r,2*r,r^2,r+1))
names(s) <- paste("Layer",1:4)
s
# class       : SpatRaster 
# dimensions  : 4, 4, 4  (nrow, ncol, nlyr)
# resolution  : 1, 1  (x, y)
# extent      : 0, 4, 0, 4  (xmin, xmax, ymin, ymax)
# coord. ref. :  
# source(s)   : memory
# names       : Layer 1, Layer 2, Layer 3, Layer 4 
# min values  :       1,       2,       1,       2 
# max values  :      16,      32,     256,      17

s+10
s > 5
s + r

В последнем случае значения каждого пикселя однослойного растра r прибавляются к значениям соответствующих пикселей для каждого слоя s.

Арифметические операции над растровыми данными с разным числом слоев выполняются по обычному правилу: растр с меньшим числом слоев циклически повторяется до тех пор, пока число его слоев не сравняется с числом слоев “более многослойного” растра, после чего выполняется нужная операция. Например:

s1 <- rast(list(r,2*r))
# class       : SpatRaster 
# dimensions  : 4, 4, 2  (nrow, ncol, nlyr)
# resolution  : 1, 1  (x, y)
# extent      : 0, 4, 0, 4  (xmin, xmax, ymin, ymax)
# coord. ref. :  
# source(s)   : memory
# names       : lyr.1, lyr.1 
# min values  :     1,     2 
# max values  :    16,    32

s + s1
# class       : SpatRaster 
# dimensions  : 4, 4, 4  (nrow, ncol, nlyr)
# resolution  : 1, 1  (x, y)
# extent      : 0, 4, 0, 4  (xmin, xmax, ymin, ymax)
# coord. ref. :  
# source(s)   : memory
# names       : Layer 1, Layer 2, Layer 3, Layer 4 
# min values  :       2,       4,       2,       4 
# max values  :      32,      64,     272,      49

Фактически, s (4 слоя) складывается не с s1 (два слоя), а с комбинацией rast(list(s1,s1)), содержащей 4 слоя.

#R

“Канопус-В-ИК” МСС (Многозональная Съемочная Система):

* Канал 1 (синий): 0,46 – 0,52 мкм
* Канал 2 (зеленый): 0,51 – 0,60 мкм
* Канал 3 (красный): 0,63 – 0,69 мкм
* Канал 4 (ближний ИК): 0,75 – 0,84 мкм

Основные характеристики целевой аппаратуры КА “Канопус-В-ИК” — ссылка.

Охват (экстент)

Равенство геометрий двух растров означает, что у них одинаковое число строк и столбцов, одинаковые размеры пикселей в одной и той же системе координат. Наверное, и охват у таких растров тоже одинаков? А вот здесь не все так просто.

Создадим тестовый растр

r <- rast(matrix(1:16, ncol = 4, nrow = 4))
# class       : SpatRaster 
# dimensions  : 4, 4, 1  (nrow, ncol, nlyr)
# resolution  : 1, 1  (x, y)
# extent      : 0, 4, 0, 4  (xmin, xmax, ymin, ymax)
# name        : lyr.1 
# min value   :     1 
# max value   :    16

(Мы показываем только нужные нам сведения про r).

Функция ext возвращает объект класса SpatExtent

class(ext(r))

Давайте немного уменьшим охват r

new_ext <- ext(r) - 0.25
# SpatExtent : 0.25, 3.75, 0.25, 3.75 (xmin, xmax, ymin, ymax)

Заметим, что произошло не просто вычитание числа — произошло именно сокращение охвата: нижние границы увеличились, тогда как верхние уменьшились.

Теперь обрежем r по уменьшенному охвату и сравним геометрии до и после обрезки:

r_crp <- crop(r, new_ext)
compareGeom(r, r_crp)
# [1] TRUE

Геометрии снова одинаковы! Выясняется, что и охват r_crp не уменьшился:

# class       : SpatRaster 
# dimensions  : 4, 4, 1  (nrow, ncol, nlyr)
# resolution  : 1, 1  (x, y)
# extent      : 0, 4, 0, 4  (xmin, xmax, ymin, ymax)

Оказывается, что пока центр ячейки (пиксела) не окажется за пределами нового охвата никакой обрезки растра не происходит. Может сами проверить это, сокращая охват r на значения, меньшие 0.5.

Как только центр ячейки окажется за пределами сократившегося охвата, эта ячейка исчезнет их растра:

new_ext <- ext(r) - 0.51
r_crp <- crop(r, new_ext)
# class       : SpatRaster 
# dimensions  : 2, 2, 1  (nrow, ncol, nlyr)
# resolution  : 1, 1  (x, y)
# extent      : 1, 3, 1, 3  (xmin, xmax, ymin, ymax)
# name        : lyr.1 
# min value   :     6 
# max value   :    11

compareGeom(r, r_crp)
# Ошибка: [compareGeom] extents do not match

Исходный растр 4 х 4, потеряв граничные ячейки, превратился в растр размера 2 х 2.

#R

На всякий случай напомним, что в закрепе есть ссылка на пост Работа с пространственными данными в R. Там собираются все наши публикации по этой теме.

Суда, проходящие через Баб-эль-Мандебский пролив, показаны на снимке Sentinel-1A, сделанном 18 декабря 2022 года (относительная орбита 116, поляризация VV).

#снимки