SpatRasterDataset и SpatRasterCollection

До сих пор мы работали с объектами класса SpatRaster.

SpatRaster представляет собой прямоугольную часть мира, разделённую на прямоугольные ячейки одинаковой площади (в единицах заданной системы координат). Для каждой ячейки может быть несколько значений (“слоёв”).

SpatRaster может указывать на один или несколько файлов на диске, в которых хранятся значения ячеек, и/или хранить эти значения в памяти. Эти объекты могут быть созданы с помощью метода rast.

Из объектов SpatRaster можно создавать новые объекты того же класса, а кроме них — SpatRasterDataset и SpatRasterCollection.

SpatRasterDataset — это набор данных, каждый элемент которого представляет собой SpatRaster для одной и той же области пространства (охвата) и системы координат, но, возможно, с разным разрешением. Элементы `SpatRasterDataset`используются для хранения отдельных переменных (например, температуры и осадков) или придания четвёртого измерения (например, высоты, глубины или времени) данным, которые уже имеют три измерения (несколько слоёв).

SpatRasterCollection — это коллекция (список) объектов SpatRaster без ограничений по протяженности или другим геометрическим параметрам. Коллекции применяют для хранения нескольких объектов SpatRaster, чтобы затем объединить их (merge) или создать из них мозаику (mosaic).

SpatRasterDataset создаётся функцией (методом) sds:

r <- rast(system.file("ex/logo.tif", package="terra"))   
x <- sds(r, r/2)
names(r) <- c("first", "second")
r

# Узнаем длину SpatRasterDataset
length(x)

# Извлечём 2-й SpatRaster
x[2]

SpatRasterCollection создаётся функцией sprc:

x <- rast(xmin=-110, xmax=-50, ymin=40, ymax=70, ncols=60, nrows=30)
y <- rast(xmin=-80, xmax=-20, ymax=60, ymin=30)
res(y) <- res(x)
values(x) <- 1:ncell(x)
values(y) <- 1:ncell(y)

z <- sprc(x, y)
z
z[1]

#R

Введение в геопространственный анализ и визуализацию в R

Журнал MDPI Data, наряду с данными, публикует статьи о методах обработки и анализа данных. В нём нам попалась забавная учебная статья

📖 Maesen, P., & Salingros, E. (2024). Introduction to Reproducible Geospatial Analysis and Figures in R: A Tutorial Article. Data, 9(4), 58. https://doi.org/10.3390/data9040058

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

В статье кратко представлены основные понятия растровых и векторных данных, системы координат, также описан базовый рабочий процесс для проведения воспроизводимых исследований в R. Приведены примеры создания различных типов карт (scatter, bubble, choropleth, hexbin и faceted) на основе открытых экологических данных. На этих примерах демонстрируются основные операции с геопространственными векторными данными (чтение, преобразование систем координат, создание геометрий, буферных зон вокруг существующих геометрий, определение пересечений между геометриями и т. п.).

Данные и скрипты есть в Supplementary Materials (скачать https://www.mdpi.com/article/10.3390/data9040058/s1).

📸 Карта средних значений органического углерода почвы (soil organic carbon, SOC) в Германии, измеренных в период 2011–2018 гг.

#R

Пакет mlhrsm для картографирования влажности почвы с высоким пространственным разрешением

Влажность почвы — одна из ключевых переменных в сельском хозяйстве и в экологии. С определением влажности почвы в масштабе поля по данным дистанционного зондирования из космоса существуют большие проблемы: исходные данные имеют низкое пространственное разрешение, а результатам машинного обучения недостаёт точности.

В работе

📖 Peng, Y., Yang, Z., Zhang, Z., & Huang, J. (2024). A Machine Learning-Based High-Resolution Soil Moisture Mapping and Spatial–Temporal Analysis: The mlhrsm Package. Agronomy, 14(3), 421. https://doi.org/10.3390/agronomy14030421

предлагается очередная модель машинного обучения для картографирования влажности почвы. Она основана на алгоритме квантильного случайного леса (quantile random forest) и использует данные наземных датчиков влажности, параметры поверхности земли (растительность, рельеф и почву), а также оценки влажности почвы на поверхности и в прикорневой зоне, полученные по спутниковым данным. В работе используются данные спутников SMAP, Sentinel-1, Landsat, а также данные приборов MODIS. Область исследования: CONUS (contiguous USA), где существуют открытые данные наземных датчиков влажности почвы.

Модель позволяет создавать карты влажности почвы высокого разрешения (от 30 до 500 м, от ежедневных до ежемесячных) и строить оценки неопределенности на участках по территории CONUS на уровнях 0–5 см и 0–1 м.

Точность результатов — примерно такая же, как и у других работ подобного рода. Примеры оценок можно посмотреть в статье.

Привлекает в работе то, что весь расчёт оформлен в виде пакета mlhrsm на языке R с открытым исходным кодом. По сути, статья — это руководство пользователя mlhrsm, где показан расчёт влажности почвы на примере одного поля. Её можно использовать для анализа сильных и слабых сторон подобных моделей, а также как основу для создания собственных моделей.

#почва #сельхоз #R

1️⃣ Скрипт для обработки ежемесячных данных.
2️⃣ Список космических и суборбитальных запусков в первом полугодии 2024 года.

#R

Скрипт для получения данных FIRMS.

#R

geodl

В R есть множество пакетов для работы с пространственными данными. А вот пакетов, где для анализа таких данных используются методы глубокого обучения (deep learning, DL), напротив, совсем мало.

Недавно появился пакет geodl, предоставляющий инструменты для семантической сегментации пространственных данных с помощью DL на основе свёрточной нейронной сети (CNN).

geodl построен на базе пакета torch, который поддерживает реализацию DL с помощью языков R и C++ без необходимости установки среды Python/PyTorch. Это значительно упрощает программную среду, необходимую для реализации DL в R. Растровые данные в geodl обрабатываются с помощью известного пакета terra, который также использует C++. Циклы обучения реализуются с помощью пакета luz.

Подробности о geodl изложены в 📖 препринте.

#R #нейронки

Векторные данные. Введение

Пространственные данные — это данные, имеющие географическую привязку. Они бывают растровыми и векторными. С растровыми данными мы уже познакомились и сейчас займёмся векторными.

Векторные данные представляют собой дискретные объекты с чёткими границами: река, дорога, город. Они состоят из 1) описания геометрии объекта и 2) переменных (атрибутов), связанных с объектом. Например, для города соответствующий векторный объект может содержать полигон, описывающий границы города (геометрию), и атрибуты, такие как “численность населения”, “площадь” и т. п.

Векторные объекты — это точки, линии, полигоны, а также наборы объектов одного типа — мультиточки, мультилинии, мультиполигоны.

Для работы с векторными данными в R широко используются пакеты terra и sf. Первый лучше подходит для операций, где одновременно нужны растровые и векторные данные. Например, чтобы обрезать снимок (растровые данные) по границам административного района (векторные данные). Второй обладает большими возможностями для выполнения операций, специфических для векторных данных. Например, для сглаживания полигонов данных.

Изучать векторные данные начнём с пакета terra. Интерфейс его нам уже знаком и многие функции, привычные по работе с растровыми данными, работают и для данных векторных.

Пакет terra определяет набор классов для представления пространственных данных, имена которых начинаются со Spat. Для векторных данных соответствующим классом является SpatVector.

Начнём с создания объектов SpatVector. Это часто делается при использовании координат, полученных с помощью глобальных навигационных спутниковых систем (GPS/ГЛОНАСС/Beidou). Также это пригодится для создания небольших примеров, иллюстрирующих ту или иную задачу. В большинстве других случаев вы будете считывать готовые объекты SpatVector из файлов или баз данных.

#R

Создание векторных объектов: точки

Зададим координаты (долготу и широту) точек и составим из них матрицу с помощью функции cbind():

longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5, -120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9, 36.2, 39, 41.6, 36.9)
lonlat <- cbind(longitude, latitude)

Для выполнения этих операций пакет terra не нужен. А вот дальше он понадобится, так что не забудьте его загрузить с помощью library(terra).

Обратите внимание, что в отличие от географии, где первой координатой указывается широта (latitude), здесь первой указана долгота (longitude). Так принято во многих компьютерных программах для работы с пространственными данными. Здесь работают традиции не географии, а математики: долгота соответствует координате Х декартовой системы координат, которая обычно в математике указывается первой.

Для создания объектов SpatVector с нуля или из файлов данных служит функция vect пакета terra

pts <- vect(lonlat)

Фактически, vect является методом класса SpatVector. Мы будем часто называть подобные методы функциями, поскольку путаницы это не создаст.

Заглянем внутрь pts:

pts
##  class       : SpatVector
##  geometry    : points
##  dimensions  : 10, 0  (geometries, attributes)
##  extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
##  coord. ref. :

geom(pts)
##       geom part      x    y hole
##  [1,]    1    1 -116.7 45.3    0
##  [2,]    2    1 -120.4 42.6    0
##  [3,]    3    1 -116.7 38.9    0
##  [4,]    4    1 -113.5 42.1    0
##  [5,]    5    1 -115.5 35.7    0
##  [6,]    6    1 -120.8 38.9    0
##  [7,]    7    1 -119.5 36.2    0
##  [8,]    8    1 -113.7 39.0    0
##  [9,]    9    1 -113.7 41.6    0
## [10,]   10    1 -110.7 36.9    0

Видно, что это объект класса SpatVector, который состоит из точек (geometry: points). Всего этих точек 10 и атрибутов у данных нет (dimensions: 10, 0 (geometries, attributes)). Экстент (extent) или охват данных автоматически рассчитан по координатам точек. В свойствах объекта указана система координат (coord. ref.). При создании pts мы её не задали. Давайте сделаем это сейчас.

Задать или узнать текущую систему координат (Coordinate Reference System, CRS) можно с помощью функции crs, знакомой нам по работе с растровыми данными:

crdref <- "+proj=longlat +datum=WGS84"
crs(pts) <- crdref

pts
## class       : SpatVector 
## geometry    : points 
## dimensions  : 10, 0  (geometries, attributes)
## extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs

Такого же результата можно было добиться, задав систему координат при создании pts

pts <- vect(lonlat, crs=crdref)
pts

Теперь добавим к нашему объекту атрибуты (переменные). Сначала создадим таблицу, число строк которой должно совпадать с числом точек (линий, полигонов, …) в векторном объекте:

precip_val <- runif(nrow(lonlat), min=0, max=100)
dt <- data.frame(ID=1:nrow(lonlat), precip=precip_val)

Здесь мы сгенерировали случайные значения осадков (precip_val) в диапазоне от 0 до 100, которых будет столько же, сколько и точек. Затем эти значения поместили в таблицу атрибутов.

Создадим новый векторный объект ptsv с заданными атрибутами

ptsv <- vect(lonlat, atts=dt, crs=crdref)

ptsv
##  class       : SpatVector
##  geometry    : points
##  dimensions  : 10, 2  (geometries, attributes)
##  extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs
##  names       :    ID precip
##  type        : <int>  <num>
##  values      :     1  98.93
##                    2  15.76
##                    3  68.81

У ptsv два атрибута: идентификатор точки (ID) и значение осадков в точке (precip).

#R

Курс “Open-Source Spatial Analytics (R)” [ссылка]

Курс посвящен изучению основ работы в свободной среде программирования R, в первую очередь, для анализа геопространственных данных. Он рассчитан на тех, кто уже имеет некоторые знания о ГИС, работе с картами, картографическими проекциями, векторными и растровыми данными. Опыт работы с R, напротив, не требуется. Предполагается, что вы обучитесь программировать на R по ходу курса.

Справочный материал представлен в виде примеров кода, видео и презентаций. Есть задания для практического обучения. В центральная колонке 📸 представлены модули курса, в левой — примеры и задания.

Курс подготовлен в West Virginia View (🔗 https://wvview.org) — это консорциум государственных, частных и некоммерческих организаций, занимающихся дистанционным зондированием, который является членом AmericaView — общеамериканской сети организаций, способствующей развитию образования в области дистанционного зондирования Земли. Руководитель исследований в West Virginia View — Аарон Максвелл (Aaron Maxwell), доцент кафедры геологии и географии Университета Западной Вирджинии.

Другие курсы West Virginia View:

💻 Methods in Open Science
👨🏻‍💻 GIScience
👨🏼‍💻 Open-Source GIScience
🛰 Remote Sensing
🌍 Digital Cartography
🌐 Client-Side Web GIS
🕸 Geospatial Deep Learning

#R #обучение

Чтение и запись векторных данных

Чтение векторных файлов осуществляет функция vect — та же, что отвечает за создание векторных данных.

Одним из распространенным форматов файлов векторных данных является шейпфайл (shapefile). Это набор из четырёх (или большего числа) файлов с одинаковыми именами, но разными расширениями. Для шейпфайла x в одной папке должны находиться: x.shp, x.shx, x.dbf и x.prj.

Откроем шейпфайл, поставляемый вместе с пакетом terra:

library(terra)
filename <- system.file("ex/lux.shp", package="terra")
## [1] "C:/Users/User/AppData/Local/R/win-library/4.3/terra/ex/lux.shp"

s <- vect(filename)
s
##  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

Функция system.file возвращает полный путь к файлу. Она нужна только для примеров работы с данными, поставляемыми с R. Для собственных файлов используйте функцию vect, указав полный путь к нужному файлу.

vect возвращает объекты SpatVector. Фактически, она создаёт эти объекты с нуля, как мы видели раньше, или из файлов векторных данных различных форматов. В нашем случае построен SpatVector, состоящий из 10 полигонов с 6 атрибутами (переменными).

Для записи векторных служит функция writeVector:

outfile <- "shp_test.shp"
writeVector(s, outfile)

Чтобы перезаписать файл поверх, нужно добавить аргумент overwrite=TRUE

writeVector(s, outfile, overwrite=TRUE)

Для удаления файлов используют функции file.remove или unlink. Будьте осторожны, не спешите!

При удалении шейпфайла нам придётся удалять сразу несколько файлов. В качестве примера удалим shp_test. Сначала мы выделим нужные файлы функцией list.files, указав шаблон имени файла, а затем удалим их при помощи file.remove

ff <- list.files(pattern="^shp_test")
ff
## [1] "shp_test.cpg" "shp_test.dbf" "shp_test.prj" "shp_test.shp" "shp_test.shx"
file.remove(ff)
## logical(0)

TRUE на выходе file.remove показывает, что заданный файл удален.

#R

Преобразование координат объектов `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

Логические операции

Рассмотрим логические операции с векторными данными (SpatVector) в пакете terra.

Начнем с функций crop и intersect — для обрезки и вычисления пересечения векторных данных соответственно:

crop(x, y, ...)
intersect(x, y, ...)

Функции принимают на вход объекты SpatVector, SpatExtent или SpatRaster (у crop первым аргументом могут быть только векторы и растры). С растровым использованием этих функций мы уже сталкивались здесь и здесь.

Построим полигон p1 и лежащий внутри него полигон p2:

p1 <- vect("POLYGON ((0 0, 8 0, 8 9, 0 9, 0 0))")
p2 <- vect("POLYGON ((2 6, 3 6, 3 8, 2 8, 2 6))")
plot(p1, lwd=2)
lines(p2, lwd=2, col="blue")

lines использована для рисования поверх первого plot’a.

Применив к p1 и p2 обрезку и пересечение

cropped <- crop(p1,p2)
isected <- intersect(p1,p2)

мы получим одинаковый результат:

 # class       : SpatVector 
 # geometry    : polygons 
 # dimensions  : 1, 0  (geometries, attributes)
 # extent      : 2, 3, 6, 8  (xmin, xmax, ymin, ymax)
 # coord. ref. :

Разница в работе функций появляется, когда у векторных данных есть атрибуты:

p1[["id1"]] <- 1L
p2[["id4"]] <- 1L

cropped <- crop(p1,p2)
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 1, 1  (geometries, attributes)
# extent      : 2, 3, 6, 8  (xmin, xmax, ymin, ymax)
# coord. ref. :  
# names       :   id1
# type        : <int>
# values      :     1
isected <- intersect(p1,p2)
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 1, 2  (geometries, attributes)
# extent      : 2, 3, 6, 8  (xmin, xmax, ymin, ymax)
# coord. ref. :  
# names       :   id1   id4
# type        : <int> <int>
# values      :     1     1

У crop, в отличие от intersect, геометрия и атрибуты y не передаются на выход.
Обрезать SpatVector можно по прямоугольнику (SpatRaster, SpatExtent) или по другому SpatVector. Если это не полигоны, то используется минимальная выпуклая оболочка.

Объединение векторов SpatVector, а также объектов SpatExtent, осуществляет функция union.

При объединении полигонов нужно учесть один важный момент. Объединим p1 и p2 с помощью union:

united <- union(p1, p2)
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 2, 2  (geometries, attributes)
# extent      : 0, 8, 0, 9  (xmin, xmax, ymin, ymax)
# coord. ref. :  
# names       :   id1   id4
# type        : <int> <int>
# values      :     1    NA
#                   1     1
plot(united[1], lwd=2, col="red")
plot(united[2], add=T, lwd=2, col="blue")

В результате получим вектор, состоящий из двух полигонов: “дырявого” p1 с вырезанным из него p2, и, собственно, p2. Если такой эффект вам не нужен, а нужно просто создать SpatVector из нескольких полигонов, используйте rbind:

merged <- rbind(p1, p2)
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 2, 2  (geometries, attributes)
# extent      : 0, 8, 0, 9  (xmin, xmax, ymin, ymax)
# coord. ref. :  
# names       :   id1   id4
# type        : <int> <int>
# values      :     1    NA
                  NA     1

Объединение линий и точек c помощью union просто объединяет два набора данных без каких-либо геометрических пересечений — так же, как и rbind. На выходе получим атрибуты обоих исходных векторов.

Если x и y имеют разный геометрический тип, то возвращается коллекция SpatVectorCollection.

Функция c() создает их векторов SpatVectorCollection или добавляет объекты в существующую коллекцию:

collected <- c(p1, p2)
# class       : SpatVectorCollection 
# length      : 2 
# geometry    : polygons (1)
               polygons (1)
# names       : ,

#R

Пространственные отношения между геометриями

Пространственные отношения между объектами — это про то, кто кого касается, пересекает, лежит внутри и т. п.

Выяснить пространственные отношения между геометриями векторов помогает функция relate. Она возвращает логическую матрицу, указывающую на наличие или отсутствие определенных пространственных отношений между геометриями x и y:

relate(x, y, relation, ...)

relation — отношение между геометриями: "intersects", "touches", "crosses", "overlaps", "within", "contains", "covers", "coveredby", "disjoint".

Посмотрим как это работает. Создадим новый полигон, лежащий внутри p1 и не пересекающийся с p2:

merged <- rbind(p1, p2)
p3 <- vect("POLYGON ((4 6, 5 6, 5 8, 4 8, 4 6))")
plot(merged, lwd=2)
lines(p3, lwd=2, col="blue")

Найдем, с какими полигонами из merged пересекается (intersects) p3:

relate(merged, p3, "intersects")
      [,1]
[1,]  TRUE
[2,] FALSE

p3 пересекается с первым элементом merged, то есть с p1, и не пересекается со вторым (p2).

is.related(x, y, relation, ...) возвращает логический вектор, указывающий на наличие/отсутствие определенных пространственных отношений между x и любой из геометрий в y.

is.related(merged, p3, "intersects")
[1]  TRUE FALSE

#R

tidyplots

Цель пакета tidyplots (https://tidyplots.org/) — упростить создание графиков для научных публикаций. Он позволяет постепенно добавлять, удалять и корректировать элементы графиков, используя последовательный и интуитивно понятный синтаксис.

Многие пакеты tidy-вселенной так или иначе связаны с Хэдли Уикемом (Hadley Wickham), чего нельзя сказать о tidyplots. Более того, его автор Ян Бродер Энглер (Jan Broder Engler) предложил грамматику графики альтернативную принятой в знаменитом пакете ggplot2, разработанном Уикемом.

Для построения графиков используются функции вида Add, Remove, Adjust, Themes, Split, Output. С Add все понятно, но зачем нужен Remove? А вот зачем: вместо того, чтобы обнулять или делать пустым ненужный элемент, функции Remove говорят: “удалить заголовок легенды” (remove_legend_title) или “удалить ось y” (remove_y_axis). Вроде бы мелочь, но кому-то именно такой способ работы придется по душе. Функции Adjust настраивают элементы графиков, например, размер шрифта или цвет. Themes — это темы, то есть типовые настройки графиков. Split позволяет строить несколько графиков одновременно. Output настраивает вывод графиков в файлы.

Подробнее о работе с пакетом читайте в:

📖 Начало работы с tidyplots

Еще подробнее — в препринте.

Ян Энглер является также автором пакета tidyheatmaps.

Спасибо коллегам за наводку!

#R