Преобразование координат объектов `SpatVector`

Для преобразования координат объектов SpatVector и SpatRaster из одной системы координат в другую в пакете terra используется функция project.

Входные параметры project():

1. x: объект SpatVector или SpatRaster, который нужно преобразовать.
2. y: новая система координат (CRS) в виде строки (WKT, PROJ4 или EPSG-кода) или объекта SpatRaster, у которого будет взята система координат.
3. method: метод интерполяции, используемый для преобразования растровых данных. По умолчанию используется "bilinear".
4. res: пространственное разрешение итогового растра (если x — объект SpatRaster).
5. size: размер итогового растра (если x является SpatRaster).
6. filename: имя файла для сохранения результата.
7. ...: Дополнительные аргументы, передаваемые другим методам.

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

Создадим вектор (SpatVector) в системе координат WGS84:

library(terra)

# Создадим вектор из WKT.
v <- vect("POLYGON ((0 -5, 10 0, 10 -10, 0 -5))", crs="EPSG:4326")

plot(v, border='blue', col='yellow', lwd=3, main = "Исходный полигон в WGS84")

Преобразуем SpatVector в новую систему координат EPSG:3857:

# Проектируем в EPSG:3857
v_projected <- project(v, "EPSG:3857")

# Проверяем старую и новую системы координат
cat(crs(v))
cat(crs(v_projected))

plot(v_projected, col = "red", main = "Перепроецированный полигон в EPSG:3857")

Преобразуем SpatVector в систему координат другого объекта SpatVector:

# Создадим новый вектор в UTM Zone 33N (EPSG:32633)
v2 <- vect(cbind(x = c(1, 2, 3), y = c(1, 2, 3)), crs = "EPSG:32633")

# Перепроецирование v в систему координат v2
v1_projected <- project(v, v2)

# Проверяем новую систему координат
cat(crs(v1_projected))

plot(v1_projected, col = "red", main = "Перепроецированный полигон в UTM Zone 33N")

#R

rsi — загрузка данных из STAC и расчет спектральных индексов [ссылка]

Пакет rsi (от repeated spatial infelicities) предоставляет пользователю:

- Интерфейс к проекту Awesome Spectral Indices project, который содержит список спектральных индексов в виде таблицы tibble.
- Метод эффективного вычисления этих спектральных индексов.
- Метод загрузки данных с любого сервера STAC, с дополнительными настройками для загрузки популярных данных Landsat, Sentinel-1 и Sentinel-2 с бесплатных и публичных серверов STAC.
- Метод объединения нескольких растров, содержащих различные наборы данных, в единый растровый стек.

Функция spectral_indices() возвращает таблицу спектральных индексов.

Функция get_stac_data() позволяет загружать изображения из любого доступного каталога STAC. Например, можно загрузить композит каналов Landsat с маской облачности:

aoi <- sf::st_point(c(-74.912131, 44.080410))
aoi <- sf::st_set_crs(sf::st_sfc(aoi), 4326)
aoi <- sf::st_buffer(sf::st_transform(aoi, 5070), 1000)

landsat_image <- get_stac_data(
  aoi,
  start_date = "2022-06-01",
  end_date = "2022-06-30",
  pixel_x_size = 30,
  pixel_y_size = 30,
  asset_names = c("red", "blue", "green"),
  stac_source = "https://planetarycomputer.microsoft.com/api/stac/v1/",
  collection = "landsat-c2-l2",
  mask_band = "qa_pixel",
  mask_function = landsat_mask_function,
  output_filename = tempfile(fileext = ".tif"),
  item_filter_function = landsat_platform_filter,
  platforms = c("landsat-9", "landsat-8")
)

Для популярных данных, например для снимков Landsat, есть отдельные функции, где большинство параметров настроено по умолчанию:

landsat_image <- get_landsat_imagery(
  aoi,
  start_date = "2022-06-01",
  end_date = "2022-06-30",
  output_filename = tempfile(fileext = ".tif")
)

По умолчанию, данные загружаются из Microsoft's Planetary Computer API.

Теперь на основе полученных каналов снимков Landsat рассчитаем спектральные индексы при помощи calculate_indices():

indices <- calculate_indices(
  landsat_image,
  available_indices,
  output_filename = tempfile(fileext = ".tif")
)

Наконец, в rsi есть утилита для эффективного объединения растров, содержащих различные данные об одном и том же месте, в VRT, что позволяет программам типа GDAL рассматривать эти отдельные источники данных как единый файл.

Например, мы можем объединить наши снимки Landsat с полученными индексами:

raster_stack <- stack_rasters(
  c(landsat_image, indices),
  tempfile(fileext = ".vrt")
)

#R #индексы

Работа с элементами SpatVector

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

Рассмотрим несколько примеров.

1️⃣ Определим число элементов векторных данных. Сначала создадим тестовый SpatVector из данных, поставляемых вместе с пакетом:

library(terra)

v <- vect(system.file("ex/lux.shp", package="terra"))

v
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 12, 6  (geometries, attributes)
# extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
# source      : lux.shp
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP
# type        : <num>    <chr> <num>    <chr> <num> <int>
#   values      :     1 Diekirch     1 Clervaux   312 18081
# 1 Diekirch     2 Diekirch   218 32543
# 1 Diekirch     3  Redange   259 18664

Вектор** v содержит 12 элементов и 6 атрибутов (переменных):

names(v)
# [1] "ID_1"   "NAME_1" "ID_2"   "NAME_2" "AREA"   "POP"

Число элементов v можно подсчитать функциями:

length(v)
nrow(v)

2️⃣ Добавление атрибута в вектор. Создавать векторные данные с заданным набором атрибутов мы умеем. Теперь добавим атрибут к уже имеющимся данным.

Добавим идентификаторы, равные номеру элемента в векторе. Сделаем это двумя способами:

v[["ID_new_1"]] <- 1:nrow(v)
v$ID_new_2 <- seq.int(nrow(v))

3️⃣ Получение координат элементов (геометрии). Координаты элементов векторов без атрибутов возвращает функция geom:

geom(v)
#      geom part        x        y hole
# [1,]    1    1 6.026519 50.17767    0
# [2,]    1    1 6.031361 50.16563    0
# [3,]    1    1 6.035646 50.16410    0
# [4,]    1    1 6.042747 50.16157    0
# [5,]    1    1 6.043894 50.16116    0
# ...

На выходе получается матрица значений координат. Или вектор (просто vector), или список, или таблица — в зависимости от настроек функции, которых очень много.

4️⃣ Конвейер функций. Конвейерная обработка функций в R (|>) встроена в язык, начиная с версии R 4.1.0. Конвейер принимает вывод одной функции и передает его в другую функцию в качестве аргумента. Иногда это делает процесс обработки данных более наглядным.

Например, вместо

v <- vect(system.file("ex/lux.shp", package="terra"))

мы могли бы записать

v <- system.file("ex/lux.shp", package="terra") |>
     vect()

* Разумеется, речь идет о перегрузке функций, точнее о перегрузке методов классов Spat* пакета terra.
** Для краткости, здесь мы называем векторные данные просто векторами.

#R

Расширения ggplot2

У популярного графического пакета ggplot2 так много расширений, что в них легко потеряться. К счастью есть веб-ресурсы, помогающие выбрать нужное расширение.

🔗ggplot2 extensions (https://exts.ggplot2.tidyverse.org) — список расширений ggplot2, который можно дополнять, и галерея примеров применения расширений.

🖥 Awesome ggplot2 — обновляемый список учебников, пакетов и других ресурсов, связанных с ggplot2.

📚 ggplot2 extended (https://ggplot2-extended-book.com/) — электронная книга от Antti Rask, посвященная расширениям ggplot2. Пока в ней в основном заготовки будущих разделов, но содержание многообещающее. Спасибо коллегам за наводку!

#R