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

Зададим координаты (долготу и широту) точек и составим из них матрицу с помощью функции 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

Introduction to Environmental Data Science

Книга Introduction to Environmental Data Science (https://bookdown.org/igisc/EnvDataSci/) Джерри Дэвиса (Jerry Davis) посвящена анализу пространственных данных, связанных с исследованиями окружающей среды, на языке программирования R.

В книге есть краткое введение в R, описана работа с векторными и растровыми данными, моделирование, анализ спутниковых снимков и временных рядов.

Environmental Data Science book (https://edsbook.org/welcome) — аналогичный ресурс на Python.

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

#R #python