class: center, middle, inverse, title-slide .title[ # Пространственные данные ] .subtitle[ ## Визуализация и анализ географических данных на языке R ] .author[ ### Тимофей Самсонов ] .institute[ ### МГУ имени Ломоносова, Географический факультет ] .date[ ### 2023-10-31 ] --- # Пространственные данные R отличается хорошо развитой экосистемой пакетов для работы с пространственными данными. Умеет работать с ними автономно и путем взаимодействия с внешними библиотеками (GDAL, GEOS, PROJ, Leaflet) и ГИС-пакетами (QGIS, SAGA, GRASS, Whitebox). <img src="img/spackages.png" width="659" /> --- # Дополнительные книги __Spatial Data Science with applications in R__ (Pebesma, Bivand): <iframe src="https://r-spatial.org/book/" width="100%" height="400px" data-external="1"></iframe> https://r-spatial.org/book/ --- # Дополнительные книги __Geocomputation with R__ (Lovelace, Nowosad, Muenchow): <iframe src="https://geocompr.robinlovelace.net" width="100%" height="400px" data-external="1"></iframe> https://geocompr.robinlovelace.net --- ## Используемые пакеты ```r library(sf) library(stars) library(mapview) library(dplyr) library(readr) ``` --- ## Модели пространственных данных __Пространственный объект__ (feature) — это цифровая модель материального или абстрактного объекта реального или виртуального мира с указанием его идентификатора, координатных и атрибутивных данных. <br> __Пространственные данные__ (spatial data) — это данные о пространственных объектах и их наборах. <br> .pull-left[ .red[_Векторная модель_] пространственных данных представляет отдельные объекты путем координатного описания их границ, а также привязанных к ним характеристик — атрибутов. ] .pull-right[ .blue[_Растровая модель_] пространственных данных представляет географическое пространство в виде непрерывного покрытия матрицей ячеек, к каждой из которых привязаны характеристики. ] --- class: center, middle # Векторные данные --- ## Simple Features __Simple Features__ (официально _Simple Features Access_) — это стандарт [OGC 06-103](http://www.opengeospatial.org/standards/sfa), который определяет общую модель хранения и доступа к векторным объектам в географических информационных системах. - Все геометрии состоят из точек. - Точки являются координатами в 2-, 3- или 4-мерном пространстве. - Все точки в геометрии имеют одинаковую размерность. В дополнение к обязательным измерениям `\(X\)` и `\(Y\)` _возможны_ два дополнительных: - `\(Z\)`, обозначающее высоту - `\(M\)`, обозначающее некоторую меру, связанную с точкой — например, время Таким образом, существует 4 варианта геометрии: `\(XY\)`, `\(XYZ\)`, `\(XYM\)` и `\(XYZM\)`. В географических координатах `\(X\)` соответствует долготе, `\(Y\)` соответствует широте. --- ## Simple Features Стандарт включает в себя 17 типов геометрий. Наиболее употребительны следующие. <br> .blue[__Простые геометрии__] Тип | Описание ----|-------------------------------------------------------------------------------- `POINT` | нуль-мерная геометрия, содержащая одну точку `LINESTRING` | последовательность точек, соединенных прямыми, несамопересекающимися отрезками; одномерная геометрия `POLYGON` | геометрия с положительной площадью (двумерная); последовательность точек, отрезки между которыми формируют замкнутое кольцо без самопересечений; первое кольцо является внешним, ноль и более остальных колец представляют дырки внутри полигона --- ## Simple Features Стандарт включает в себя 17 типов геометрий. Наиболее употребительны следующие. <br> .red[__Мультигеометрии__] Тип | Описание ----|-------------------------------------------------------------------------------- `MULTIPOINT` | множество точек; геометрия типа `MULTIPOINT` называется _простой_ если ни одна пара точек в `MULTIPOINT` не совпадает `MULTILINESTRING` | множество линий `MULTIPOLYGON` | множество полигонов `GEOMETRYCOLLECTION` | множество геометрий произвольного типа за исключением `GEOMETRYCOLLECTION` Оставшиеся виды геометрий _Simple Features_ включают: `CIRCULARSTRING`, `COMPOUNDCURVE`, `CURVEPOLYGON`, `MULTICURVE`, `MULTISURFACE`, `CURVE`, `SURFACE`, `POLYHEDRALSURFACE`, `TIN`, `TRIANGLE`. --- ## Simple Features <!-- --> --- ## Форматы представления .blue[__Well-Known Text (WKT)__] — текстовый формат (удобен для визуализации) ``` ## POINT (0.5 0.5) ## LINESTRING (0 1, 0.5 1.5, 1.2 1.2, 2 1.3, 3 2) ## POLYGON ((0.5 0.5, 2 0, 3 2, 1.5 4, 0 3, 0.5 0.5), (1 1, 0.8 2, 2 2.2, 1.4 1.1, 1 1)) ## MULTIPOINT ((0.5 0.5), (1 3), (2 1), (0.2 2), (2 3), (1.5 1.5)) ## MULTILINESTRING ((0.5 1.5, 1.2 1.2, 2 1.3), (0 1.5, 0.5 2, 1.2 1.7), (2 1.8, 3 2.5)) ## MULTIPOLYGON (((0.5 0.5, 2 0, 3 2, 1.5 4, 0 3, 0.5 0.5), (1 1, 0.8 2, 2 2.2, 1.4 1.1, 1 1)), ((3 3.3, 3.5 3.1, 4 3, 4 3.7, 3.7 3.96, 3.2 4, 3 3.3), (3.2 3.4, 3.8 3.2, 3.8 3.7, 3.3 3.8, 3.2 3.4)), ((3 1.2, 2.5 0.2, 3.5 0.2, 3.5 1.2, 3 1.2)), ((0 1, 0.1 0.8, 0.2 0.5, 0.1 0.3, 0 0.7, 0 1))) ## GEOMETRYCOLLECTION (POLYGON ((0.5 0.5, 2 0, 3 2, 1.5 4, 0 3, 0.5 0.5), (1 1, 0.8 2, 2 2.2, 1.4 1.1, 1 1)), MULTIPOINT ((3.5 -0.5), (4 2), (5 0), (3.2 1), (5 2), (4.5 0.5)), MULTILINESTRING ((3 3.5, 3.7 3.2, 4.5 3.3), (2.5 3.5, 3 4, 3.7 3.7), (4.5 3.8, 5.5 4.5))) ``` .red[__Well-Known Binary (WKB)__] — бинарный формат (предпочтителен для хранения) ``` ## 01 01 00 00 00 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 e0 3f ## 01 02 00 00 00 05 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 f0 3f 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 f8 3f 33 33 33 33 33 33 f3 3f 33 33 33 33 33 33 f3 3f 00 00 00 00 00 00 00 40 cd cc cc cc cc cc f4 3f 00 00 00 00 00 00 08 40 00 00 00 00 00 00 00 40 ## 01 03 00 00 00 02 00 00 00 06 00 00 00 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 00 40 00 00 00 00 00 00 00 00 00 00 00 00 00 00 08 40 00 00 00 00 00 00 00 40 00 00 00 00 00 00 f8 3f 00 00 00 00 00 00 10 40 00 00 00 00 00 00 00 00 00 00 00 00 00 00 08 40 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 e0 3f 05 00 00 00 00 00 00 00 00 00 f0 3f 00 00 00 00 00 00 f0 3f 9a 99 99 99 99 99 e9 3f 00 00 00 00 00 00 00 40 00 00 00 00 00 00 00 40 9a 99 99 99 99 99 01 40 66 66 66 66 66 66 f6 3f 9a 99 99 99 99 99 f1 3f 00 00 00 00 00 00 f0 3f 00 00 00 00 00 00 f0 3f ## 01 04 00 00 00 06 00 00 00 01 01 00 00 00 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 e0 3f 01 01 00 00 00 00 00 00 00 00 00 f0 3f 00 00 00 00 00 00 08 40 01 01 00 00 00 00 00 00 00 00 00 00 40 00 00 00 00 00 00 f0 3f 01 01 00 00 00 9a 99 99 99 99 99 c9 3f 00 00 00 00 00 00 00 40 01 01 00 00 00 00 00 00 00 00 00 00 40 00 00 00 00 00 00 08 40 01 01 00 00 00 00 00 00 00 00 00 f8 3f 00 00 00 00 00 00 f8 3f ## 01 05 00 00 00 03 00 00 00 01 02 00 00 00 03 00 00 00 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 f8 3f 33 33 33 33 33 33 f3 3f 33 33 33 33 33 33 f3 3f 00 00 00 00 00 00 00 40 cd cc cc cc cc cc f4 3f 01 02 00 00 00 03 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 f8 3f 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 00 40 33 33 33 33 33 33 f3 3f 33 33 33 33 33 33 fb 3f 01 02 00 00 00 02 00 00 00 00 00 00 00 00 00 00 40 cd cc cc cc cc cc fc 3f 00 00 00 00 00 00 08 40 00 00 00 00 00 00 04 40 ## 01 06 00 00 00 04 00 00 00 01 03 00 00 00 02 00 00 00 06 00 00 00 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 00 40 00 00 00 00 00 00 00 00 00 00 00 00 00 00 08 40 00 00 00 00 00 00 00 40 00 00 00 00 00 00 f8 3f 00 00 00 00 00 00 10 40 00 00 00 00 00 00 00 00 00 00 00 00 00 00 08 40 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 e0 3f 05 00 00 00 00 00 00 00 00 00 f0 3f 00 00 00 00 00 00 f0 3f 9a 99 99 99 99 99 e9 3f 00 00 00 00 00 00 00 40 00 00 00 00 00 00 00 40 9a 99 99 99 99 99 01 40 66 66 66 66 66 66 f6 3f 9a 99 99 99 99 99 f1 3f 00 00 00 00 00 00 f0 3f 00 00 00 00 00 00 f0 3f 01 03 00 00 00 02 00 00 00 07 00 00 00 00 00 00 00 00 00 08 40 66 66 66 66 66 66 0a 40 00 00 00 00 00 00 0c 40 cd cc cc cc cc cc 08 40 00 00 00 00 00 00 10 40 00 00 00 00 00 00 08 40 00 00 00 00 00 00 10 40 9a 99 99 99 99 99 0d 40 9a 99 99 99 99 99 0d 40 ae 47 e1 7a 14 ae 0f 40 9a 99 99 99 99 99 09 40 00 00 00 00 00 00 10 40 00 00 00 00 00 00 08 40 66 66 66 66 66 66 0a 40 05 00 00 00 9a 99 99 99 99 99 09 40 33 33 33 33 33 33 0b 40 66 66 66 66 66 66 0e 40 9a 99 99 99 99 99 09 40 66 66 66 66 66 66 0e 40 9a 99 99 99 99 99 0d 40 66 66 66 66 66 66 0a 40 66 66 66 66 66 66 0e 40 9a 99 99 99 99 99 09 40 33 33 33 33 33 33 0b 40 01 03 00 00 00 01 00 00 00 05 00 00 00 00 00 00 00 00 00 08 40 33 33 33 33 33 33 f3 3f 00 00 00 00 00 00 04 40 9a 99 99 99 99 99 c9 3f 00 00 00 00 00 00 0c 40 9a 99 99 99 99 99 c9 3f 00 00 00 00 00 00 0c 40 33 33 33 33 33 33 f3 3f 00 00 00 00 00 00 08 40 33 33 33 33 33 33 f3 3f 01 03 00 00 00 01 00 00 00 06 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 f0 3f 9a 99 99 99 99 99 b9 3f 9a 99 99 99 99 99 e9 3f 9a 99 99 99 99 99 c9 3f 00 00 00 00 00 00 e0 3f 9a 99 99 99 99 99 b9 3f 33 33 33 33 33 33 d3 3f 00 00 00 00 00 00 00 00 66 66 66 66 66 66 e6 3f 00 00 00 00 00 00 00 00 00 00 00 00 00 00 f0 3f ## 01 07 00 00 00 03 00 00 00 01 03 00 00 00 02 00 00 00 06 00 00 00 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 00 40 00 00 00 00 00 00 00 00 00 00 00 00 00 00 08 40 00 00 00 00 00 00 00 40 00 00 00 00 00 00 f8 3f 00 00 00 00 00 00 10 40 00 00 00 00 00 00 00 00 00 00 00 00 00 00 08 40 00 00 00 00 00 00 e0 3f 00 00 00 00 00 00 e0 3f 05 00 00 00 00 00 00 00 00 00 f0 3f 00 00 00 00 00 00 f0 3f 9a 99 99 99 99 99 e9 3f 00 00 00 00 00 00 00 40 00 00 00 00 00 00 00 40 9a 99 99 99 99 99 01 40 66 66 66 66 66 66 f6 3f 9a 99 99 99 99 99 f1 3f 00 00 00 00 00 00 f0 3f 00 00 00 00 00 00 f0 3f 01 04 00 00 00 06 00 00 00 01 01 00 00 00 00 00 00 00 00 00 0c 40 00 00 00 00 00 00 e0 bf 01 01 00 00 00 00 00 00 00 00 00 10 40 00 00 00 00 00 00 00 40 01 01 00 00 00 00 00 00 00 00 00 14 40 00 00 00 00 00 00 00 00 01 01 00 00 00 9a 99 99 99 99 99 09 40 00 00 00 00 00 00 f0 3f 01 01 00 00 00 00 00 00 00 00 00 14 40 00 00 00 00 00 00 00 40 01 01 00 00 00 00 00 00 00 00 00 12 40 00 00 00 00 00 00 e0 3f 01 05 00 00 00 03 00 00 00 01 02 00 00 00 03 00 00 00 00 00 00 00 00 00 08 40 00 00 00 00 00 00 0c 40 9a 99 99 99 99 99 0d 40 9a 99 99 99 99 99 09 40 00 00 00 00 00 00 12 40 66 66 66 66 66 66 0a 40 01 02 00 00 00 03 00 00 00 00 00 00 00 00 00 04 40 00 00 00 00 00 00 0c 40 00 00 00 00 00 00 08 40 00 00 00 00 00 00 10 40 9a 99 99 99 99 99 0d 40 9a 99 99 99 99 99 0d 40 01 02 00 00 00 02 00 00 00 00 00 00 00 00 00 12 40 66 66 66 66 66 66 0e 40 00 00 00 00 00 00 16 40 00 00 00 00 00 00 12 40 ``` --- ## Базовые библиотеки .left-40[ В R существует высоко развитая инфраструктура для работы с векторными данными, которая обеспечивается пакетом [__sf__](https://cran.r-project.org/web/packages/sf/index.html). <br> __sf__ опирается на библиотеки [PROJ](https://proj.org), [GDAL](https://gdal.org), [GEOS](https://trac.osgeo.org/geos/) и [S2](https://s2geometry.io), которые устанавливаются вместе с ним. ] .right-60[ <!-- --> ] --- ## Чтение Для чтения данных средствами __sf__ необходимо использовать функцию `st_read()`: ```r countries = st_read('../r-geo-course/data/ne/countries.gpkg') ## Reading layer `admin_0_map_units' from data source ## `/Users/tsamsonov/GitHub/r-geo-course/data/ne/countries.gpkg' ## using driver `GPKG' ## Simple feature collection with 183 features and 72 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513 ## Geodetic CRS: WGS 84 ``` - Коллекция из 183 пространственных объектов с 72 атрибутами - Тип геометрии `MULTIPOLYGON`, размерность геометрии `\(XY\)` - Ограничивающий прямоугольник (разброс координат) по осям `\(X\)` и `\(Y\)` имеет диапазон `\([-180, 180] \times [-90, 83.64513]\)` - Проекция (CRS — coordinate reference system) имеет название _WGS 84_. --- ## Чтение Подгрузим также для работы данные по другим типам объектов: ```r oceans = st_read('../r-geo-course/data/ne/oceans.gpkg', quiet = TRUE) rivers = st_read('../r-geo-course/data/ne/rivers.gpkg', quiet = TRUE) lakes = st_read('../r-geo-course/data/ne/lakes.gpkg', quiet = TRUE) cities = st_read('../r-geo-course/data/ne/cities.gpkg', quiet = TRUE) ``` > Параметр `quiet = TRUE` отключает вывод информации о загруженных данных. --- ## Внутренняя структура Представление пространственных объектов типа Simple Features реализовано в виде иерархии из трех классов объектов: 1. `sf` (simple features) — объект класса `data.frame`, представляющий множество пространственных объектов со списком-колонкой для хранения геометрии. 1. `sfc` (simple features geometry column) — список-колонка в объекте `sf`, представляющий множество геометрий пространственных объектов. 1. `sfg` (simple feature geometry) — геометрия пространственного объекта внутри списка `sfc`. Поскольку Simple Features реализованы в виде обычных фреймов данных, _любая операция, применимая к фрейму данных, будет также применима к объекту типа_ `sf`: ```r class(countries) ## [1] "sf" "data.frame" ``` --- ## Внутренняя структура Геометрия пространственных объектов хранится в одном из столбцов. В данном случае он имеет название `geometry`: ```r head(countries['geometry']) ## Simple feature collection with 6 features and 0 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -73.41544 ymin: -55.25 xmax: 75.15803 ymax: 42.68825 ## Geodetic CRS: WGS 84 ## geometry ## 1 MULTIPOLYGON (((61.21082 35... ## 2 MULTIPOLYGON (((23.90415 -1... ## 3 MULTIPOLYGON (((21.02004 40... ## 4 MULTIPOLYGON (((51.57952 24... ## 5 MULTIPOLYGON (((-66.95992 -... ## 6 MULTIPOLYGON (((43.58275 41... ``` --- ## Внутренняя структура Геометрический столбец можно извлечь, применив функцию `st_geometry()`. Полученный объект будет иметь тип __sfc__ (Simple Feature Geometry Column): ```r outlines = st_geometry(countries) class(outlines) ## [1] "sfc_MULTIPOLYGON" "sfc" head(outlines) ## Geometry set for 6 features ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -73.41544 ymin: -55.25 xmax: 75.15803 ymax: 42.68825 ## Geodetic CRS: WGS 84 ## First 5 geometries: ``` В данном случае объекты имеют класс `sfc_MULTIPOLYGON`, который является расширением класса `sfc` (simple feature geometry column). --- ## Внутренняя структура .pull-left[ Поскольку объект класса `sfc` представляет собой список, любой элемент (отдельный объект) можно извлечь по его порядковому номеру: ```r obj = outlines[[8]] class(obj) ``` [1] "XY" "MULTIPOLYGON" "sfg" <br> Геометрия 8-го объекта имеет класс `sfg`, реализованный в виде мультиполигонов (`MULTIPOLYGON`) с плоскими координатами (`XY`) ] .pull-right[ Чтобы добраться до координат, необходимо развернуть иерархию списков, из которых состоит `sfg`: ```r outlines[[8]][[1]] ## [[1]] ## [,1] [,2] ## [1,] 68.9350 -48.6250 ## [2,] 69.5800 -48.9400 ## [3,] 70.5250 -49.0650 ## [4,] 70.5600 -49.2550 ## [5,] 70.2800 -49.7100 ## [6,] 68.7450 -49.7750 ## [7,] 68.7200 -49.2425 ## [8,] 68.8675 -48.8300 ## [9,] 68.9350 -48.6250 ``` ] --- ## Визуализация (базовая графика) По умолчанию выполняется графический обзор атрибутов объектов ```r plot(countries) ``` <!-- --> --- ## Визуализация (базовая графика) Если задача стоит нарисовать границы объектов, то нужно отображать объект __sfc__: ```r plot(st_geometry(countries)) ``` <!-- --> --- ## Визуализация (базовая графика) Для быстрого построения тематических карт по выбранному показателю необходимо при вызове функции `plot()` указать соответствующий атрибут фрейма данных: ```r plot(countries['sovereignt'], key.pos = NULL) # Здесь легенда не нужна ``` <!-- --> --- ## Визуализация (базовая графика) Для отображения координатной сетки надо указать параметр `graticule = TRUE`, а подписей координат — `axes = TRUE`: ```r plot(countries['gdp_md_est'], graticule = TRUE, axes = TRUE) ``` <!-- --> --- ## Визуализация (базовая графика) Для совмещения нескольких слоев на одной карте необходимо при втором и последующих вызовах функции `plot()` указать параметр `add = TRUE`. .pull-left[ .code-small[ ```r cities_large = cities |> filter(scalerank == 0, ! name %in% c('Washington, D.C.', 'Paris', 'Riyadh', 'Rome', 'São Paulo', 'Kolkata')) plot(st_geometry(countries), lwd = 0.5, border = 'gray', xlim = c(-180,180)) plot(oceans, col = 'steelblue1', border = 'steelblue', add = TRUE) plot(lakes, col = 'steelblue1', border = 'steelblue', add = TRUE) plot(rivers, col = 'steelblue', add = TRUE) plot(cities_large, col = 'black', pch = 19, cex = 0.25, add = TRUE) text(cities_large$longitude, cities_large$latitude, label = cities_large$name, cex = 0.5, pos = 2, offset = 0.25) ``` ] ] .pull-right[ <!-- --> ] --- ## Визуализация (интерактивная) Для интерактивного исследования объектов можно использовать пакет `mapview`: ```r mapview(countries) ``` <img src="img/mapview1.png" width="481" /> --- ## Визуализация (интерактивная) Отображение показателей цветом осуществляется через параметр `col.regions`: ```r nconts = length(unique(countries$continent)) mapview(countries, zcol = 'continent', col.regions = RColorBrewer::brewer.pal(nconts, 'Set1')) ``` <img src="img/mapview2.png" width="481" /> --- ## Визуализация (интерактивная) Несколько слоев складываются через `+` как в __ggplot2__: .code-small[ ```r mapview(countries, zcol = 'continent', * col.regions = RColorBrewer::brewer.pal(nconts, 'Set1')) + mapview(cities_large, col.regions = 'black', label = 'name', cex = 3) } |> leafem::addStaticLabels(cities_large, label = cities_large$name, offset = c(0.1, 0), style = list("color" = "black", "font-weight" = "bold")) ``` ] <img src="img/mapview3.png" width="481" /> --- ## Атрибутивные операции К объектам `sf` применимы стандартные преобразования фреймов данных. .pull-left[ __Италия__: ```r italy = countries |> filter(sovereignt == 'Italy') plot(st_geometry(italy)) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-28-1.png" width="250px" /> ] .pull-right[ __Cтраны более 100 млн. жителей__: ```r largest = countries |> select(pop_est) |> filter(pop_est > 100000000) plot(outlines, col = 'lightgrey') plot(largest, col = 'red', add = TRUE) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-29-1.png" width="400px" /> ] --- ## Атрибутивные операции .pull-left[ __Валовой региональный продукт по континентам:__ ```r continents = countries %>% filter(., st_is_valid(.)) |> group_by(continent) |> summarise(gdp = sum(gdp_md_est)) plot(continents['gdp']) ``` ] .pull-right[ <img src="09_SpatialData_files/figure-html/unnamed-chunk-30-1.png" width="400px" /> ] --- ## Геометрические объекты (sfg) Для создания геометрических объектов класса `sfg` существует ряд функций: Функция | Тип пространственного объекта --------------------------|----------------------------------- `st_point()` | _POINT_ `st_linestring()` | _LINESTRING_ `st_polygon()` | _POLYGON_ `st_multipoint()` | _MULTIPOINT_ `st_multilinestring()` | _MULTILINESTRING_ `st_multipolygon()` | _MULTIPOLYGON_ `st_geometrycollection()` | _GEOMETRYCOLLECTION_ <br><br><br><br><br><br><br><br><br><br><br> В зависимости от типа создаваемого объекта, данные функции принимают координаты в виде вектора, матрицы или списка. --- ## Точки Проще всего создаются отдельные __точки__ (_POINT_): ```r st_point(c(0, 2)) # XY POINT ## POINT (0 2) st_point(c(0, 2, -1)) # XYZ POINT ## POINT Z (0 2 -1) st_point(c(0, 2, 5), dim = 'XYM') # XYM POINT ## POINT M (0 2 5) st_point(c(0, 2, -1, 5)) # XYZM POINT ## POINT ZM (0 2 -1 5) ``` Дополнительный параметр `dim=` служит для уточнения типа геометрии точек и по сути нужен только тогда, когда необходимо создать редко используемые точки типа _XYM_. Во всех остальных случаях (_XY_, _XYZ_, _XYZM_) размерность геометрии распознается по умолчанию. --- ## Мультиточки и линии При создании __мультиточек__ (`MULTIPOINT`) и __линий__ (`LINESTRING`) необходимо подавать на вход функции уже матрицу координат: ```r coords = matrix(c( 0, 2, 1, 3, 3, 1, 5, 0 ), ncol = 2, byrow = TRUE) mp = st_multipoint(coords) # XY MULTIPOINT print(mp) ## MULTIPOINT ((0 2), (1 3), (3 1), (5 0)) ls = st_linestring(coords) # XY LINESTRING print(ls) ## LINESTRING (0 2, 1 3, 3 1, 5 0) ``` --- ## Мультиточки и линии .pull-left[ В первом случае геометрия состоит из отдельных точек. Во втором случае те же самые точки соединены линией: ```r plot(ls) plot(mp, col = 'red', pch = 19, add = TRUE) ``` > .small[Создание трех-(_XYZ_, _XYM_) и четырехмерных (_ZYXM_) мультиточек и линий выполняется аналогично, но матрица должна содержать не 2, а, соответственно 3 или 4 столбца, и при необходимости параметр `dim = 'XYM'`.] ] .pull-right[ <!-- --> ] --- ## Полигоны Полигоны создаются на основе списков матриц координат. .pull-left[ ```r # Координаты границы coords = matrix(c( 1, 0, 0, 2, 2, 3, 4, 2, 3, 0.5, 1, 0 ), ncol = 2, byrow = TRUE) # Простой полигон без дырки pol = st_polygon(list(coords)) print(pol) ## POLYGON ((1 0, 0 2, 2 3, 4 2, 3 0.5, 1 0)) ``` ] .pull-right[ ```r plot(pol, col = 'lightblue') ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-35-1.png" width="400px" /> ] --- ## Полигоны Вторая и последующие матрицы координат в списке определяют границы дырок. .pull-left[ ```r # Координаты дырки hole = matrix(c( 2, 1, 3, 1.5, 3, 2, 2, 2, 1.5, 1.5, 2, 1 ), ncol = 2, byrow = TRUE) # Полигон с дыркой pol2 = st_polygon(list(coords, hole)) print(pol2) ## POLYGON ((1 0, 0 2, 2 3, 4 2, 3 0.5, 1 0), (2 1, 3 1.5, 3 2, 2 2, 1.5 1.5, 2 1)) ``` ] .pull-right[ ```r plot(pol2, col = 'lightblue') ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-37-1.png" width="400px" /> ] --- ## Мультиполигоны .pull-left[ ```r coords1 = matrix(c( 0.5, 0, 0, 1, 1, 1.5, 2, 1, 1.5, 0.25, 0.5, 0 ), ncol = 2, byrow = TRUE) coords2 = matrix(c( 3, 1, 2.5, 2, 3.5, 2.5, 4, 2, 4, 1.25, 3, 1 ), ncol = 2, byrow = TRUE) ``` ] .pull-right[ Мультиполигоны (`MULTIPOLYGON`) и мультилинии (`MULTILINESTRING`) требуются, когда один географический объект состоит из нескольких геометрических объектов. В мультиполигонах добавляется еще один уровень вложенности списков. ```r mpol = st_multipolygon( list( list(coords1), list(coords2) ) ) print(mpol) ## MULTIPOLYGON (((0.5 0, 0 1, 1 1.5, 2 1, 1.5 0.25, 0.5 0)), ((3 1, 2.5 2, 3.5 2.5, 4 2, 4 1.25, 3 1))) ``` ] --- ## Мультиполигоны ```r plot(pol, col = 'grey') # Обычный полигон (серый) plot(mpol, col = 'pink', add = TRUE) # Мультиполигон (розовый) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-40-1.png" width="500px" /> --- ## Мультиполигоны Пример из практики: остров на озере как часть суши. .pull-left[ ```r coords4 = matrix(c( 2.2, 1.2, 2.8, 1.5, 2.8, 1.8, 2.2, 1.8, 2.0, 1.6, 2.2, 1.2 ), ncol = 2, byrow = TRUE) island = st_polygon(list(coords4)) mpol2 = st_multipolygon( list(pol2, island) ) print(mpol2) ## MULTIPOLYGON (((1 0, 0 2, 2 3, 4 2, 3 0.5, 1 0), (2 1, 3 1.5, 3 2, 2 2, 1.5 1.5, 2 1)), ((2.2 1.2, 2.8 1.5, 2.8 1.8, 2.2 1.8, 2 1.6, 2.2 1.2))) ``` ] .pull-right[ ```r plot(mpol2, col = 'olivedrab') ``` <!-- --> ] --- ## Мультилинии Мультилиния, в отличие от мультиполигона, не требует дополнительного списка верхнего уровня, поскольку линии не могут содержать дыр. Например, можно собрать мультилинию из двух частей, соответствующих участкам реки до и после озера: .pull-left[ ```r coords1 = matrix(c( -3, 0, -1, 2, 0, 2 ), ncol = 2, byrow = TRUE) coords2 = matrix(c( 4, 2, 5, 3, 6, 5 ), ncol = 2, byrow = TRUE) mline = st_multilinestring(list(coords1, coords2)) print(mline) ## MULTILINESTRING ((-3 0, -1 2, 0 2), (4 2, 5 3, 6 5)) ``` ] .pull-right[ ```r plot(mline, lwd = 3, col = 'blue') plot(pol2, col = 'lightblue', add = TRUE) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-44-1.png" height="250px" /> ] --- ## Геометрические коллекции Геометрическая коллекция (`GEOMETRYCOLLECTION`) позволяет хранить вместе любые виды геометрий. Коллекции редко создаются напрямую, чаще они получаются в результате выполнения геометрических операций типа оверлея. ```r col = st_geometrycollection(list(ls, mp, mline, pol2)) print(col) ## GEOMETRYCOLLECTION (LINESTRING (0 2, 1 3, 3 1, 5 0), MULTIPOINT ((0 2), (1 3), (3 1), (5 0)), MULTILINESTRING ((-3 0, -1 2, 0 2), (4 2, 5 3, 6 5)), POLYGON ((1 0, 0 2, 2 3, 4 2, 3 0.5, 1 0), (2 1, 3 1.5, 3 2, 2 2, 1.5 1.5, 2 1))) ``` .pull-left[ ```r plot(col) ``` ] .pull-right[ <img src="09_SpatialData_files/figure-html/unnamed-chunk-46-1.png" height="200px" /> ] --- ## Списки геометрических объектов (sfc) Списки геометрических объектов (класс `sfc`) используются в таблицах пространственных объектов в качестве столбца, который хранит геометрию объектов. ```r moscow.sfg = st_point(c(37.615, 55.752)) irkutsk.sfg = st_point(c(104.296, 52.298)) petro.sfg = st_point(c(158.651, 53.044)) cities.sfc = st_sfc(moscow.sfg, irkutsk.sfg, petro.sfg) print(cities.sfc) ## Geometry set for 3 features ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: 37.615 ymin: 52.298 xmax: 158.651 ymax: 55.752 ## CRS: NA ## POINT (37.615 55.752) ## POINT (104.296 52.298) ## POINT (158.651 53.044) ``` --- ## Списки геометрических объектов (sfc) При создании списка геометрий для него может быть определена система координат: ```r st_crs(cities.sfc) = st_crs(4326) # WGS84 print(cities.sfc) ## Geometry set for 3 features ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: 37.615 ymin: 52.298 xmax: 158.651 ymax: 55.752 ## Geodetic CRS: WGS 84 ``` > Для списка геометрий может быть определена только одна система координат --- ## Списки геометрических объектов (sfc) Куда легли созданные точки? ```r plot(cities.sfc, pch = 19) countries |> filter(sovereignt == 'Russia') |> st_geometry() |> plot(add = TRUE) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-49-1.png" height="250px" /> --- ## Пространственные объекты (`sf`) Пространственные объекты (класс `sf`) организуются в виде фрейма данных, один из столбцов которого имеет класс `sfc`. ```r city.attr = data.frame( name = c('Москва', 'Иркутск', 'Петропавловск-Камчатский'), established = c(1147, 1661, 1740), population = c(12500, 620, 180) ) cites.sf = st_sf(city.attr, geometry = cities.sfc) print(cites.sf) ## Simple feature collection with 3 features and 3 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: 37.615 ymin: 52.298 xmax: 158.651 ymax: 55.752 ## Geodetic CRS: WGS 84 ## name established population geometry ## 1 Москва 1147 12500 POINT (37.615 55.752) ## 2 Иркутск 1661 620 POINT (104.296 52.298) ## 3 Петропавловск-Камчатский 1740 180 POINT (158.651 53.044) ``` --- ## Точки по координатам Для решения этой задачи можно воспользоваться функцией `st_as_sf()`. Рассмотрим задачу на примере файла координат станций из базы метеорологических данных [__ВНИИГМИ-МЦД__](http://meteo.ru/): ```r stations = read_fwf( '../r-geo-course/data/vniigmi/stations.txt', col_positions = fwf_widths(diff(c(1, 7, 42, 47, 53, 59, 67, 71)), col_names = c('id', 'name', 'lat', 't1', 'lon', 't2', 'z')), locale = locale(encoding = 'CP1251') ) print(stations) ## # A tibble: 1,124 × 7 ## id name lat t1 lon t2 z ## <dbl> <chr> <dbl> <chr> <dbl> <chr> <dbl> ## 1 20046 Им.Э.Т.Кренкеля,ГМО 80.6 с.ш. 58 в.д. 21 ## 2 20069 Остров_Визе 79.5 с.ш. 77.0 в.д. 10 ## 3 20087 Голомянный 79.6 с.ш. 90.6 в.д. 7 ## 4 20107 Баренцбург 78.1 с.ш. 14.2 в.д. 73 ## 5 20289 Русский 77.2 с.ш. 96.4 в.д. 9 ## 6 20292 Им.Е.К.Федорова,ГМО 77.7 с.ш. 104. в.д. 12 ## 7 20353 мыс_Желания 77.0 с.ш. 68.6 в.д. 9 ## 8 20476 Стерлегова 75.4 с.ш. 88.9 в.д. 10 ## 9 20667 Им.М.В.Попова 73.3 с.ш. 70.0 в.д. 4 ## 10 20674 Остров_Диксон 73.5 с.ш. 80.4 в.д. 42 ## # ℹ 1,114 more rows ``` --- ## Точки по координатам Теперь создадим пространственные точки на основе этой таблицы, взяв координаты из столбцов _lat_ и _lon_ соответственно и указав код системы координат: .pull-left[ ```r sf_stations = st_as_sf( stations, coords = c("lon", "lat"), crs = 4326 ) plot( st_geometry(sf_stations), pch = 19, col = 'red', cex = 0.25 ) plot( st_geometry(countries), border = 'grey', add = TRUE ) box() ``` ] .pull-right[ <img src="09_SpatialData_files/figure-html/unnamed-chunk-52-1.png" height="450px" /> ] --- ## Преобразование типа геометрии Для преобразования типов геометрии существует функция `st_cast()`. Функция принимает объекты классов `sfg`, `sfc` или `sf`, а также название типа геометрии, к которому необходимо привести входные объекты. ```r italy.borders = st_cast(italy, 'MULTILINESTRING') class(st_geometry(italy.borders)) ## [1] "sfc_MULTILINESTRING" "sfc" italy.regions = st_cast(italy.borders, 'MULTIPOLYGON') class(st_geometry(italy.regions)) ## [1] "sfc_MULTIPOLYGON" "sfc" italy.points = st_cast(italy.borders, 'POINT') class(st_geometry(italy.points)) ## [1] "sfc_POINT" "sfc" ``` --- ## Преобразование типа геометрии Для преобразования типов геометрии существует функция `st_cast()`. Функция принимает объекты классов `sfg`, `sfc` или `sf`, а также название типа геометрии, к которому необходимо привести входные объекты. .pull-left[ ```r plot( st_geometry(italy.regions), lwd = 0.5 ) plot( italy.points, pch = 20, add = TRUE ) ``` ] .pull-right[ <img src="09_SpatialData_files/figure-html/unnamed-chunk-54-1.png" height="400px" /> ] --- ## Полигонизация и разбиение линий __Полигонизация__ — это процесс преобразования линии или мультилинии в полигон(ы). Может быть полезной если необходимо восстановить площади объектов на основе их границ. На вход нужна _одна_ мультилиния: .pull-left[ ```r # Создадим три линии coords1 = rbind(c(0, 0), c(0, 6)) line1 = st_linestring(coords1) coords2 = rbind(c(-1,1), c(5,1)) line2 = st_linestring(coords2) coords3 = rbind(c(-1,5), c(4,0)) line3 = st_linestring(coords3) # Создадим мультилинию mls = st_multilinestring( list(line1, line2, line3) ) ``` ] .pull-right[ ```r plot(mls) # мультилиния plot(st_cast(mls, 'MULTIPOINT'), pch = 20, cex = 2, add = TRUE) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-56-1.png" height="250px" /> ] --- ## Полигонизация и разбиение линий .pull-left[ Для того чтобы мультилинию превратить в полигоны, необходимо чтобы в местах пересечений и самопересечений стояли узла. Добавить узлы можно с помощью функции `st_node()`: ```r # Неудачная попытка st_polygonize(mls) ## GEOMETRYCOLLECTION EMPTY # Добавление пересечений mls2 = st_node(mls) poly2 = st_polygonize(mls2) print(poly2) ## GEOMETRYCOLLECTION (POLYGON ((0 1, 0 4, 3 1, 0 1))) ``` ] .pull-right[ ```r plot(mls2) plot(poly2, col = 'grey', add = TRUE) plot(st_cast(mls2, 'MULTIPOINT'), pch = 20, cex = 2, add = TRUE) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-58-1.png" height="300px" /> ] --- ## Геометрические атрибуты К описательным характеристикам геометрии относятся ограничивающий прямоугольник, периметр (для линий и полигонов), площадь (для полигонов), центроид и список координат. ```r st_bbox(italy) # Органичивающий прямоугольник ## xmin ymin xmax ymax ## 6.749955 36.619987 18.480247 47.115393 st_area(italy) # Площадь ## 314577521836 [m^2] st_length(italy.borders) # Периметр ## 5319369 [m] st_centroid(italy) |> # Центроид (может быть не внутри для невыпуклых фигур) st_geometry() ## Geometry set for 1 feature ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: 12.2687 ymin: 42.67074 xmax: 12.2687 ymax: 42.67074 ## Geodetic CRS: WGS 84 ``` --- ## Геометрические атрибуты К описательным характеристикам геометрии относятся ограничивающий прямоугольник, периметр (для линий и полигонов), площадь (для полигонов), центроид и список координат. ```r st_point_on_surface(italy) |> # Точка гарантированно внутри st_geometry() ## Geometry set for 1 feature ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: 12.63118 ymin: 42.55822 xmax: 12.63118 ymax: 42.55822 ## Geodetic CRS: WGS 84 st_coordinates(italy) |> head(4) # Список координат ## X Y L1 L2 L3 ## [1,] 10.44270 46.89355 1 1 1 ## [2,] 11.04856 46.75136 1 1 1 ## [3,] 11.16483 46.94158 1 1 1 ## [4,] 12.15309 47.11539 1 1 1 ``` --- ## Геометрические атрибуты Пример — ограничивающий прямоугольник и центроиды для Индонезии. .pull-left[ .code-small[ ```r indonesia = countries |> filter(sovereignt == 'Indonesia') box = st_as_sfc(st_bbox(indonesia)) plot(indonesia |> st_geometry(), col = 'lightgrey') plot(box, border = 'red', add = TRUE) plot(st_centroid(indonesia), col = 'darkgreen', pch = 19, add = TRUE) plot(st_point_on_surface(indonesia), col = 'violetred', pch = 19, add = TRUE) ``` ] ] .pull-left[ <br> <img src="09_SpatialData_files/figure-html/unnamed-chunk-61-1.png" width="600px" /> ] --- ## Экспорт Для экспорта векторных пространственных данных можно воспользоваться функцией `st_write()`, которая определит формат выходного файла по указанному вами расширению: ```r st_write(cites.sf, 'data/mycities.shp') # Шейп-файл st_write(cites.sf, 'data/topo.gpkg', 'mycities') # Geopackage ``` --- class: center, middle # Растровые данные --- ## Растр и его геометрия __Растр__ представляет из себя матрицу значений. Каждой ячейке матрицы соответствует прямоугольная пространственная область фиксированного размера, которая называется _пикселом_. <br><br> __Геометрия растра__ определяет, где именно располагаются в пространстве пикселы растра и может быть описана путем указания следующих компонент <br><br> Параметр | Назначение ---------|--------- `NCOLS` | Количество столбцов `NROWS` | Количество строк `XLLCENTER` | Координата `\(X\)` центра левой нижней ячейки растра `YLLCENTER` | Координата `\(Y\)` центра левой нижней ячейки растра `CELLSIZE` | Размер ячейки --- ## Пакеты для работы с растрами В настоящее время для работы с растровыми данными в R используются два пакета: [__stars__](https://r-spatial.github.io/stars/) и [__terra__](https://cran.r-project.org/web/packages/terra/index.html). - __terra__ является наследником пакета [raster](https://cran.r-project.org/web/packages/terra/index.html), который исторически был основным средством работы с растровыми данными и обладает широким спектром функций растрового анализа. - __stars__ — относительно новый, разработан с целью поддержки многомерных данных и более тесного взаимодействия с пакетом `sf`. В целом можно сказать, что пакеты __terra__ и __stars__ частично пересекаются по функциональности, но скорее дополняют друг друга, нежели дублируют. > В этой и ближайших лекциях мы будем работать с растрами в формате __stars__, поскольку он концептуально близок к пакету __sf__. --- ## Чтение Для чтения растров любой размерности можно использовать функцию `read_stars()`: ```r dem = read_stars('../r-geo-course/data/world/gebco.tif') # Цифровая модель рельефа img = read_stars('../r-geo-course/data/world/BlueMarbleJuly.tif') # Цветной космический снимок (RGB) img ## stars object with 3 dimensions and 1 attribute ## attribute(s): ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## BlueMarbleJuly.tif 1 13 33 63.13569 75 255 ## dimension(s): ## from to offset delta refsys point x/y ## x 1 720 -180 0.5 WGS 84 FALSE [x] ## y 1 360 90 -0.5 WGS 84 FALSE [y] ## band 1 3 NA NA NA NA ``` --- ## Внутренняя структура Для работы с данными типа stars необходимо понимать их внутреннюю структуру. Для начала можно взглянуть на нее посредством стандартной функции `str()`: ```r str(img) ## List of 1 ## $ BlueMarbleJuly.tif: num [1:720, 1:360, 1:3] 6 4 7 7 7 7 7 8 8 8 ... ## - attr(*, "dimensions")=List of 3 ## ..$ x :List of 7 ## .. ..$ from : num 1 ## .. ..$ to : num 720 ## .. ..$ offset: num -180 ## .. ..$ delta : num 0.5 ## .. ..$ refsys:List of 2 ## .. .. ..$ input: chr "WGS 84" ## .. .. ..$ wkt : chr "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic Sys"| __truncated__ ## .. .. ..- attr(*, "class")= chr "crs" ## .. ..$ point : logi FALSE ## .. ..$ values: NULL ## .. ..- attr(*, "class")= chr "dimension" ## ..$ y :List of 7 ## .. ..$ from : num 1 ## .. ..$ to : num 360 ## .. ..$ offset: num 90 ## .. ..$ delta : num -0.5 ## .. ..$ refsys:List of 2 ## .. .. ..$ input: chr "WGS 84" ## .. .. ..$ wkt : chr "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic Sys"| __truncated__ ## .. .. ..- attr(*, "class")= chr "crs" ## .. ..$ point : logi FALSE ## .. ..$ values: NULL ## .. ..- attr(*, "class")= chr "dimension" ## ..$ band:List of 7 ## .. ..$ from : num 1 ## .. ..$ to : int 3 ## .. ..$ offset: num NA ## .. ..$ delta : num NA ## .. ..$ refsys: chr NA ## .. ..$ point : logi NA ## .. ..$ values: NULL ## .. ..- attr(*, "class")= chr "dimension" ## ..- attr(*, "raster")=List of 4 ## .. ..$ affine : num [1:2] 0 0 ## .. ..$ dimensions : chr [1:2] "x" "y" ## .. ..$ curvilinear: logi FALSE ## .. ..$ blocksizes : int [1:3, 1:2] 720 720 720 1 1 1 ## .. .. ..- attr(*, "dimnames")=List of 2 ## .. .. .. ..$ : NULL ## .. .. .. ..$ : chr [1:2] "x" "y" ## .. ..- attr(*, "class")= chr "stars_raster" ## ..- attr(*, "class")= chr "dimensions" ## - attr(*, "class")= chr "stars" ``` --- ## Внутренняя структура Данный трёхканальный растр представляет собой список из единственного элемента с названием `BlueMarbleJuly.tif` — это имя было присвоено автоматически при чтении растра. Каждый такой элемент соответствует _переменной_ данных. <br><br> В данном случае переменная одна — это интенсивность цвета. Хранится она в виде трехмерного массива (`array`) размерностью `\(720 \times 360 \times 3\)`: ```r str(img[[1]]) ## num [1:720, 1:360, 1:3] 6 4 7 7 7 7 7 8 8 8 ... img[[1]][100, 200, 2] ## [1] 14 ``` --- ## Измерения (dimensions) Каждой оси этого массива внутри `stars` соответствует измерение (`dimension`), которое определяет параметры отображения индексов массива на соответствующую систему координат (пространственную, временную, спектральную и т.д.). <br><br> Список измерений и его элементы можно получить так: ```r attr(img, 'dimensions') ## from to offset delta refsys point x/y ## x 1 720 -180 0.5 WGS 84 FALSE [x] ## y 1 360 90 -0.5 WGS 84 FALSE [y] ## band 1 3 NA NA NA NA attr(img, 'dimensions')[['x']] ## $from ## [1] 1 ## ## $to ## [1] 720 ## ## $offset ## [1] -180 ## ## $delta ## [1] 0.5 ## ## $refsys ## Coordinate Reference System: ## User input: WGS 84 ## wkt: ## GEOGCRS["WGS 84", ## ENSEMBLE["World Geodetic System 1984 ensemble", ## MEMBER["World Geodetic System 1984 (Transit)"], ## MEMBER["World Geodetic System 1984 (G730)"], ## MEMBER["World Geodetic System 1984 (G873)"], ## MEMBER["World Geodetic System 1984 (G1150)"], ## MEMBER["World Geodetic System 1984 (G1674)"], ## MEMBER["World Geodetic System 1984 (G1762)"], ## MEMBER["World Geodetic System 1984 (G2139)"], ## ELLIPSOID["WGS 84",6378137,298.257223563, ## LENGTHUNIT["metre",1]], ## ENSEMBLEACCURACY[2.0]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433]], ## CS[ellipsoidal,2], ## AXIS["geodetic latitude (Lat)",north, ## ORDER[1], ## ANGLEUNIT["degree",0.0174532925199433]], ## AXIS["geodetic longitude (Lon)",east, ## ORDER[2], ## ANGLEUNIT["degree",0.0174532925199433]], ## USAGE[ ## SCOPE["Horizontal component of 3D system."], ## AREA["World."], ## BBOX[-90,-180,90,180]], ## ID["EPSG",4326]] ## ## $point ## [1] FALSE ## ## $values ## NULL ## ## attr(,"class") ## [1] "dimension" ``` --- ## Измерения (dimensions) Атрибут `dimensions` имеет свой собственный атрибут `raster`, который необходим для того чтобы определить какие именно измерения растра являются пространственными, а также установить преобразования, которые будут над ними производиться при анализе или визулизации: ```r img |> attr('dimensions') |> attr('raster') |> str() ## List of 4 ## $ affine : num [1:2] 0 0 ## $ dimensions : chr [1:2] "x" "y" ## $ curvilinear: logi FALSE ## $ blocksizes : int [1:3, 1:2] 720 720 720 1 1 1 ## ..- attr(*, "dimnames")=List of 2 ## .. ..$ : NULL ## .. ..$ : chr [1:2] "x" "y" ## - attr(*, "class")= chr "stars_raster" ``` --- ## Визуализация растров (статичная) .pull-left[ Для визуализации одноканальных растров используется функция `plot()`. В простейшем виде ей достаточно просто передать визуализируемый растр: ```r plot(dem) ``` ] .pull-right[ <img src="09_SpatialData_files/figure-html/unnamed-chunk-68-1.png" height="400px" /> ] --- ## Визуализация растров (статичная) .pull-left[ Задав шкалу и соответствующие параметры, можно сделать визуализацию более информативной: .code-small[ ```r brks = c(-12000, -5000, -2500, -1000, -200, 0, 200, 500, 1000, 2000, 4000, 8000) clrs = c("steelblue4", "steelblue3", "steelblue2", "steelblue1", "lightskyblue1", "darkseagreen", "lightgoldenrod1", "darkgoldenrod1", "darkorange", "coral2", "firebrick3") plot(dem, breaks = brks, col = clrs) ``` ] ] .pull-right[ <img src="09_SpatialData_files/figure-html/unnamed-chunk-69-1.png" height="400px" /> ] --- ## Визуализация растров (статичная) Для синтезирования цветного изображения на основе многоканального растра необходимо объект `stars` предварительно подать в функцию `st_rgb()`: ```r plot(st_rgb(img)) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-70-1.png" height="300px" /> --- ## Визуализация растров (статичная) Для синтеза в заданной последовательности каналов достаточно их перемешать в последнем индексе: .pull-left[ .code-small[ ```r st_rgb(img[,,,c(3, 1, 2)]) |> plot() ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-71-1.png" height="300px" /> ] ] .pull-right[ .code-small[ ```r st_rgb(img[,,,c(2, 3, 1)]) |> plot() ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-72-1.png" height="300px" /> ] ] --- ## Визуализация растров (статичная) Вы можете совмещать на картах несколько растровых и векторных слоев точно так же как и при совмещении векторных данных (указав параметр `add = TRUE` при вызове функции `plot()`): ```r plot(st_rgb(img), reset = FALSE) plot(outlines, border = rgb(1,1,1,0.5), lwd = 0.5, add = TRUE) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-73-1.png" height="300px" /> --- ## Визуализация растров (интерактивная) Объекты типа `stars` могут быть визуализированы аналогично векторным на интерактивных картах `mapview`: ```r mapview(dem, at = brks, col = clrs) ``` <img src="img/mapview4.png" width="482" /> --- ## Обрезка __Ограничивающим прямоугольником:__ ```r box = st_bbox(c(xmin=-80, xmax=-10, ymax=85, ymin=58), crs=st_crs(4326)) *dem_greenland = dem[box] plot(dem_greenland) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-76-1.png" height="300px" /> --- ## Обрезка __Произвольным контуром:__ ```r country = countries |> filter(name == 'Afghanistan') *dem_country = dem[country] plot(dem_country) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-77-1.png" height="300px" /> --- ## Индексирование Ортогональная структура объектов типа stars позволяет выполнять по ним различные срезы, отсекая ненужные данные. Для этого используется оператор квадратной скобки `[`, который работает следующим образом: - первый аргумент выбирает атрибут - второй и последующий аргументы выбирают измерения. Таким образом, при работе с растрами, которые содержат один атрибут, вам необходимо указать 4 индекса: `[var, x, y, band]`, где `var` - это название или порядковый номер атрибута, а `x, y, band` — порядковые номера двух координатных и одного семантического измерения. --- ## Индексирование .pull-left[ __Выбрать 1 канал:__ ```r ch1 = img[,,,1] ch1 ## stars object with 3 dimensions and 1 attribute ## attribute(s): ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## BlueMarbleJuly.tif 1 5 14 51.1141 47 255 ## dimension(s): ## from to offset delta refsys point x/y ## x 1 720 -180 0.5 WGS 84 FALSE [x] ## y 1 360 90 -0.5 WGS 84 FALSE [y] ## band 1 1 NA NA NA NA ``` ] .pull-right[ ```r plot(ch1) ``` <!-- --> ] --- ## Индексирование .pull-left[ __Выбрать диапазон ячеек растра__ ```r frag = img[, 320:470, 100:255, ] frag ## stars object with 3 dimensions and 1 attribute ## attribute(s): ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## BlueMarbleJuly.tif 1 20 50 65.96065 106 221 ## dimension(s): ## from to offset delta refsys point x/y ## x 320 470 -180 0.5 WGS 84 FALSE [x] ## y 100 255 90 -0.5 WGS 84 FALSE [y] ## band 1 3 NA NA NA NA ``` ] .pull-right[ ```r plot(st_rgb(frag)) ``` <img src="09_SpatialData_files/figure-html/unnamed-chunk-81-1.png" height="400px" /> ] --- ## Манипуляции Объекты типа `stars` поддерживают манипуляции, аналогичные тем, что могут применяться к векторным данным. Посмотрим это на примере данных по высоте земной поверхности с учетом и без покровного оледенения: ```r etopo = read_stars(c( '../r-geo-course/data/etopo1_bed.tif', '../r-geo-course/data/etopo1_ice.tif' )) etopo ## stars object with 2 dimensions and 2 attributes ## attribute(s): ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## bed -10632 -4287 -2451.5 -2113.199 86 6159 ## ice -10632 -4287 -2451.5 -1892.726 215 6159 ## dimension(s): ## from to offset delta refsys point x/y ## x 1 720 -180 0.5 WGS 84 FALSE [x] ## y 1 360 90 -0.5 WGS 84 FALSE [y] ``` --- ## Манипуляции _Переименуем_ каналы в `ice` и `bed` и посчитаем толщину покровного оледеления как их разность через _мутирование_: .pull-left[ ```r etopo = etopo |> rename(bed = 1, ice = 2) |> mutate(depth = ice - bed) plot(etopo['depth'], col = cm.colors(5), breaks = c(0, 500, 1000, 2000, 3000, 4000), main = 'Мощность покровного оледенения', reset = FALSE) plot(oceans, col = 'steelblue', add = TRUE) ``` ] .pull-right[ <img src="09_SpatialData_files/figure-html/unnamed-chunk-83-1.png" height="400" /> ] --- ## Манипуляции _Фильтрация_ происходит по измерениям, но применяется не к индексам ячеек, а к соответствующим величинам измерений: .pull-left[ ```r greenland = etopo |> filter(x > -80, x < -10, y > 58, y < 85) plot(greenland) ``` ] .pull-right[ <img src="09_SpatialData_files/figure-html/unnamed-chunk-84-1.png" height="400" /> ] --- ## Манипуляции _Выбор_ переменной позволяет оставить только ее: ```r icedepth = etopo |> select(depth) icedepth ## stars object with 2 dimensions and 1 attribute ## attribute(s): ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## depth -198 0 0 220.4736 0 4286 ## dimension(s): ## from to offset delta refsys point x/y ## x 1 720 -180 0.5 WGS 84 FALSE [x] ## y 1 360 90 -0.5 WGS 84 FALSE [y] ``` --- ## Экспорт Чтобы экспортировать (сохранить в файл) любой растр, можно воспользоваться функцией `write_stars()`, указав имя выходного файла: ```r write_stars(greenland, 'data/world/greenland.tif') ``` --- class: center, middle --- ## Пространственная привязка __Пространственная привязка__ (_spatial reference_ или _georeference_) — важнейшая составляющая пространственных данных, которая говорит о том, как правильно интерпретировать координаты объектов. Пространственная привязка в простейшем случае включает несколько фундаментальных компонент: 1. _Эллипсоид вращения_ — тело, по отношению к которому вычисляются геодезические координаты точек (широты и долготы) 1. _Исходные геодезические даты (датум)_ — параметры положения эллипсоида в теле Земли 1. _Географическая система координат_ — включает датум, положение начального меридиана и единицы измерения широт и долгот 1. _Проекция_ — математический способ перехода от географических координат на эллипсоиде к плоским прямоугольным координатам карты. 1. _Плоская прямоугольная система координат_ — включает проекцию, ее параметры и единицы измерения координат. --- ## Пространственная привязка Если точки имеют также координаты `\(Z\)`, то для их правильной интерпретации необходимы дополнительные компоненты пространственной привязки: 1. _Система счета высот_ (геодезические, нормальные, ортометрические) - определяет метод вычисления высот и глубин (координата `\(Z/H\)`) 1. _Модель геоида, квазигеоида или эллипсоида_ — определяет поверхность, относительно которой вычисляются высоты. 1. _Вертикальная система координат_ — фактическая реализация системы счета высот относительно конкретной поверхности относимости с заданным положением нулевого уровня. Например, в России это _Балтийская система нормальных высот_ с нулем в г. Кронштадт. Аналогичным образом требуется введение системы счета дополнительных координат `\(M\)`, если они используются в представлении координат. --- ## Форматы описания пространственной привязки Существует три распространенных способа задания (хранения) пространственной привязки: - _WKT (Well-Known Text)_ — представление в виде иерархического списка. Это _наиболее полный_ формат описания пространственной привязки, который рекомендуется к использованию для избежания неоднозначностей. - _PROJ.4 String_ — представление в виде строки. - _EPSG (European Petroleum Survey Group)_ — представление в виде числового кода. Для поиска проекций в перечисленных форматах представления удобно воспользоваться порталом [spatialreference.org](spatialreference.org). --- ## PROJ.4 String __PROJ.4 String__ — строковый формат представления информации о пространственной привязки, введенный в библиотеке [__PROJ__](http://proj.org), начиная с версии 4. Основные параметры: ``` +datum Datum name (see `proj -ld`) +ellps Ellipsoid name (see `proj -le`) +lat_0 Latitude of origin +lat_1 Latitude of first standard parallel +lat_2 Latitude of second standard parallel +lat_ts Latitude of true scale +lon_0 Central meridian +proj Projection name (see `proj -l`) +units meters, US survey feet, etc. +vunits vertical units. +x_0 False easting +y_0 False northing +zone UTM zone ``` --- ## PROJ.4 String Примеры записи координат в формате PROJ.4: - Географические координаты в _WGS84_ (без проекции): ``` ## +proj=longlat +datum=WGS84 +no_defs ``` - Координаты в проекции _Web Mercator_ (проекция Google Maps, Яндекс.Карт и т.д.): ``` ## +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs ``` - Координаты в _конической равнопромежуточной проекции_: ``` ## +proj=eqdc +lat_0=0 +lon_0=0 +lat_1=60 +lat_2=60 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs ``` - Координаты в проекции _UTM, зона 37_: ``` ## +proj=utm +zone=37 +datum=WGS84 +units=m +no_defs ``` --- ## WKT (Well-Known Text) __WKT__ предполагает представление компонент пространственной привязки к виде иерархического списка. Например, _полярная стереографическая проекция_: .code-vsmall[ ``` PROJCS["WGS 84 / EPSG Russia Polar Stereographic", GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]], PROJECTION["Polar_Stereographic"], PARAMETER["latitude_of_origin",90], PARAMETER["central_meridian",105], PARAMETER["scale_factor",0.994], PARAMETER["false_easting",2000000], PARAMETER["false_northing",2000000], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AXIS["X",EAST], AXIS["Y",NORTH], AUTHORITY["EPSG","5940"]] ``` ] --- ## EPSG (European Petroleum Survey Group) __EPSG (European Petroleum Survey Group)__ — европейская рабочая группа нефтегазовой области, которая ведет реестр систем координат с уникальными цифровыми кодами вида `EPSG:xxxxxx`. <br> Коды EPSG оказались удобны, поэтому используются повсеместно для быстрой инициализации проекций со стандартными параметрами. Например, ранее рассмотренные проекции имеют следующие коды EPSG: - _WGS84_: `EPSG:4326` - _Web Mercator_: `EPSG:3857` - _UTM_: `EPSG:326..` , например для UTM 37N: `EPSG:32637` --- ## Преобразование координат Преобразование координат включает три различных операции: 1. __Трансформирование__ — пересчет географических координат с одного датума на другой 1. __Проецирование__ — от географических координат к плоским прямоугольным 1. __Обратное проецирование__ — от плоских координат к географическим. Например, чтобы пересчитать координаты _UTM_ в проекцию _Гаусса-Крюгера_, необходимо: 1. Обратно проецировать координаты в географические _WGS84_ 1. Трансформировать географические координаты c _WGS84_ в _ГСК-2011_ 1. Проецировать координаты _ГСК-2011_ в проекцию _Гаусса-Крюгера_ _Несоответствие датумов часто является причиной того, что данные из разных наборов плохо совмещаются друг с другом_ --- ## Операции с пространственной привязкой Работа с пространственной привязкой данных в R состоит в основном из четырех операций: - чтение информации о системе координат - создание информации о системе координат - замена информации о системе координат - изменение системы координат (проецирование) Первые три операции (чтение, создание, замена) осуществляются функцией `st_crs()`. --- ## Чтение системы координат Чтобы прочитать информацию о проекции, достаточно передать в качестве параметра объект типа `sf` или `stars`: .code-small[ ```r st_crs(countries) # Координатная система ## Coordinate Reference System: ## User input: WGS 84 ## wkt: ## GEOGCRS["WGS 84", ## ENSEMBLE["World Geodetic System 1984 ensemble", ## MEMBER["World Geodetic System 1984 (Transit)"], ## MEMBER["World Geodetic System 1984 (G730)"], ## MEMBER["World Geodetic System 1984 (G873)"], ## MEMBER["World Geodetic System 1984 (G1150)"], ## MEMBER["World Geodetic System 1984 (G1674)"], ## MEMBER["World Geodetic System 1984 (G1762)"], ## MEMBER["World Geodetic System 1984 (G2139)"], ## ELLIPSOID["WGS 84",6378137,298.257223563, ## LENGTHUNIT["metre",1]], ## ENSEMBLEACCURACY[2.0]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433]], ## CS[ellipsoidal,2], ## AXIS["geodetic latitude (Lat)",north, ## ORDER[1], ## ANGLEUNIT["degree",0.0174532925199433]], ## AXIS["geodetic longitude (Lon)",east, ## ORDER[2], ## ANGLEUNIT["degree",0.0174532925199433]], ## USAGE[ ## SCOPE["Horizontal component of 3D system."], ## AREA["World."], ## BBOX[-90,-180,90,180]], ## ID["EPSG",4326]] ``` ] --- ## Создание системы координат Эта же функция позволяет создать новую координатную систему, путем передачи ей кода _EPSG_ или строки _PROJ.4_: .code-small[ ```r st_crs(3857) # Проекция Меркатора для карт мира ## Coordinate Reference System: ## User input: EPSG:3857 ## wkt: ## PROJCRS["WGS 84 / Pseudo-Mercator", ## BASEGEOGCRS["WGS 84", ## ENSEMBLE["World Geodetic System 1984 ensemble", ## MEMBER["World Geodetic System 1984 (Transit)"], ## MEMBER["World Geodetic System 1984 (G730)"], ## MEMBER["World Geodetic System 1984 (G873)"], ## MEMBER["World Geodetic System 1984 (G1150)"], ## MEMBER["World Geodetic System 1984 (G1674)"], ## MEMBER["World Geodetic System 1984 (G1762)"], ## MEMBER["World Geodetic System 1984 (G2139)"], ## ELLIPSOID["WGS 84",6378137,298.257223563, ## LENGTHUNIT["metre",1]], ## ENSEMBLEACCURACY[2.0]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433]], ## ID["EPSG",4326]], ## CONVERSION["Popular Visualisation Pseudo-Mercator", ## METHOD["Popular Visualisation Pseudo Mercator", ## ID["EPSG",1024]], ## PARAMETER["Latitude of natural origin",0, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8801]], ## PARAMETER["Longitude of natural origin",0, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8802]], ## PARAMETER["False easting",0, ## LENGTHUNIT["metre",1], ## ID["EPSG",8806]], ## PARAMETER["False northing",0, ## LENGTHUNIT["metre",1], ## ID["EPSG",8807]]], ## CS[Cartesian,2], ## AXIS["easting (X)",east, ## ORDER[1], ## LENGTHUNIT["metre",1]], ## AXIS["northing (Y)",north, ## ORDER[2], ## LENGTHUNIT["metre",1]], ## USAGE[ ## SCOPE["Web mapping and visualisation."], ## AREA["World between 85.06°S and 85.06°N."], ## BBOX[-85.06,-180,85.06,180]], ## ID["EPSG",3857]] ``` ] --- ## Создание системы координат Эта же функция позволяет создать новую координатную систему, путем передачи ей кода _EPSG_ или строки _PROJ.4_: .code-small[ ```r st_crs('+proj=robin') # Проекция Робинсона для карт мира ## Coordinate Reference System: ## User input: +proj=robin ## wkt: ## PROJCRS["unknown", ## BASEGEOGCRS["unknown", ## DATUM["World Geodetic System 1984", ## ELLIPSOID["WGS 84",6378137,298.257223563, ## LENGTHUNIT["metre",1]], ## ID["EPSG",6326]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8901]]], ## CONVERSION["unknown", ## METHOD["Robinson"], ## PARAMETER["Longitude of natural origin",0, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8802]], ## PARAMETER["False easting",0, ## LENGTHUNIT["metre",1], ## ID["EPSG",8806]], ## PARAMETER["False northing",0, ## LENGTHUNIT["metre",1], ## ID["EPSG",8807]]], ## CS[Cartesian,2], ## AXIS["(E)",east, ## ORDER[1], ## LENGTHUNIT["metre",1, ## ID["EPSG",9001]]], ## AXIS["(N)",north, ## ORDER[2], ## LENGTHUNIT["metre",1, ## ID["EPSG",9001]]]] ``` ] --- ## Создание системы координат Эта же функция позволяет создать новую координатную систему, путем передачи ей кода _EPSG_ или строки _PROJ.4_: .code-small[ ```r # Проекция UTM, зона 37. st_crs('+proj=utm +zone=37 +datum=WGS84 +units=m') ## Coordinate Reference System: ## User input: +proj=utm +zone=37 +datum=WGS84 +units=m ## wkt: ## PROJCRS["unknown", ## BASEGEOGCRS["unknown", ## DATUM["World Geodetic System 1984", ## ELLIPSOID["WGS 84",6378137,298.257223563, ## LENGTHUNIT["metre",1]], ## ID["EPSG",6326]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8901]]], ## CONVERSION["UTM zone 37N", ## METHOD["Transverse Mercator", ## ID["EPSG",9807]], ## PARAMETER["Latitude of natural origin",0, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8801]], ## PARAMETER["Longitude of natural origin",39, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8802]], ## PARAMETER["Scale factor at natural origin",0.9996, ## SCALEUNIT["unity",1], ## ID["EPSG",8805]], ## PARAMETER["False easting",500000, ## LENGTHUNIT["metre",1], ## ID["EPSG",8806]], ## PARAMETER["False northing",0, ## LENGTHUNIT["metre",1], ## ID["EPSG",8807]], ## ID["EPSG",16037]], ## CS[Cartesian,2], ## AXIS["(E)",east, ## ORDER[1], ## LENGTHUNIT["metre",1, ## ID["EPSG",9001]]], ## AXIS["(N)",north, ## ORDER[2], ## LENGTHUNIT["metre",1, ## ID["EPSG",9001]]]] ``` ] --- ## Замена системы координат Замена координатной системы требуется в тех случаях, когда слой не имеет пространственной привязки, или же она задана _некорректно_. В этом случае необходимо вызвать для слоя функцию `st_crs()` и перезаписать результат. ```r st_crs(countries) = NA st_crs(countries) ## Coordinate Reference System: NA st_crs(countries) = st_crs(4326) st_crs(countries) ## Coordinate Reference System: ## User input: EPSG:4326 ## wkt: ## GEOGCRS["WGS 84", ## ENSEMBLE["World Geodetic System 1984 ensemble", ## MEMBER["World Geodetic System 1984 (Transit)"], ## MEMBER["World Geodetic System 1984 (G730)"], ## MEMBER["World Geodetic System 1984 (G873)"], ## MEMBER["World Geodetic System 1984 (G1150)"], ## MEMBER["World Geodetic System 1984 (G1674)"], ## MEMBER["World Geodetic System 1984 (G1762)"], ## MEMBER["World Geodetic System 1984 (G2139)"], ## ELLIPSOID["WGS 84",6378137,298.257223563, ## LENGTHUNIT["metre",1]], ## ENSEMBLEACCURACY[2.0]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433]], ## CS[ellipsoidal,2], ## AXIS["geodetic latitude (Lat)",north, ## ORDER[1], ## ANGLEUNIT["degree",0.0174532925199433]], ## AXIS["geodetic longitude (Lon)",east, ## ORDER[2], ## ANGLEUNIT["degree",0.0174532925199433]], ## USAGE[ ## SCOPE["Horizontal component of 3D system."], ## AREA["World."], ## BBOX[-90,-180,90,180]], ## ID["EPSG",4326]] ``` --- ## Проецирование Для проецирования данных в другую систему координат следует использовать функцию `st_tranform(x, crs)`: .pull-left[ ```r # Проекция Миллера countries_mill = st_transform(countries, 'ESRI:54003') par(mar = c(2,16,2,16)) plot(st_geometry(countries_mill), col = 'lightgray', lwd = 0.5, graticule = TRUE, axes = TRUE) ``` ] .pull-right[ <!-- --> ] --- ## Проецирование Для проецирования данных в другую систему координат следует использовать функцию `st_tranform(x, crs)`: .pull-left[ ```r # Проекция Мольвейде countries_moll = countries |> st_transform('+proj=moll') plot(st_geometry(countries_moll), col = 'lightgray', lwd = 0.5, graticule = TRUE, axes = TRUE) ``` ] .pull-right[ <!-- --> ] --- ## Проецирование Для проецирования данных в другую систему координат следует использовать функцию `st_tranform(x, crs)`: .pull-left[ ```r # Коническая равнопромежуточная countries_eqdc = countries |> filter(continent == 'Europe' & sovereignt != 'Russia') |> st_transform('+proj=eqdc +lon_0=10 +lat_1=30 +lat_2=60 +datum=WGS84 +units=m') plot(st_geometry(countries_eqdc), col = 'lightgray', lwd = 0.5, graticule = TRUE, axes = TRUE) ``` ] .pull-right[ <!-- --> ] --- ## Система координат растровых данных Работа с системой координат __растровых данных__ также предполагает четыре возможных процедуры: чтение, создание, замена и проецирование. Для первых трех операций, как и прежде, используется `st_crs()`: ```r st_crs(dem) # читаем систему координат ## Coordinate Reference System: ## User input: WGS 84 ## wkt: ## GEOGCRS["WGS 84", ## ENSEMBLE["World Geodetic System 1984 ensemble", ## MEMBER["World Geodetic System 1984 (Transit)"], ## MEMBER["World Geodetic System 1984 (G730)"], ## MEMBER["World Geodetic System 1984 (G873)"], ## MEMBER["World Geodetic System 1984 (G1150)"], ## MEMBER["World Geodetic System 1984 (G1674)"], ## MEMBER["World Geodetic System 1984 (G1762)"], ## MEMBER["World Geodetic System 1984 (G2139)"], ## ELLIPSOID["WGS 84",6378137,298.257223563, ## LENGTHUNIT["metre",1]], ## ENSEMBLEACCURACY[2.0]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433]], ## CS[ellipsoidal,2], ## AXIS["geodetic latitude (Lat)",north, ## ORDER[1], ## ANGLEUNIT["degree",0.0174532925199433]], ## AXIS["geodetic longitude (Lon)",east, ## ORDER[2], ## ANGLEUNIT["degree",0.0174532925199433]], ## USAGE[ ## SCOPE["Horizontal component of 3D system."], ## AREA["World."], ## BBOX[-90,-180,90,180]], ## ID["EPSG",4326]] ``` --- ## Проецирование растровых данных Для проецирования растра в новую систему координат необходимо использовать функцию `st_warp()`. .pull-left[ ```r # Проекция Миллера img_mill = st_warp( img, crs = 'ESRI:54003' ) plot(st_rgb(img_mill), main = NULL, reset = FALSE) plot(st_geometry(countries_mill), border = rgb(1,1,1,0.5), lwd = 0.25, add = TRUE) ``` ] --- ## Проецирование растровых данных Для проецирования растра в новую систему координат необходимо использовать функцию `st_warp()`. .pull-left[ ```r # Проекция Мольвейде img_moll = st_warp( img, crs = st_crs('+proj=moll') ) plot(st_rgb(img_moll, probs = c(0.01, 0.99), stretch = "percent"), main = NULL, reset = FALSE) plot(st_geometry(countries_moll), border = rgb(1,1,1,0.5), lwd = 0.5, add = TRUE) ``` ] .pull-right[ <!-- --> ] --- ## Проецирование растровых данных Для проецирования растра в новую систему координат необходимо использовать функцию `st_warp()`. .pull-left[ ```r # Проекция коническая равнопромежуточная prj = '+proj=eqdc +lon_0=10 +lat_1=30 +lat_2=60 +datum=WGS84 +units=m' img_eqdc = st_warp( img, crs = st_crs(prj) ) img_eqdc_euro = img_eqdc[st_bbox(countries_eqdc)] plot(st_rgb(img_eqdc_euro), main = NULL, reset = FALSE) plot(st_geometry(countries_eqdc), border = rgb(1,1,1,0.5), lwd = 0.5, add = TRUE) ``` ] .pull-right[ <img src="09_SpatialData_files/figure-html/unnamed-chunk-101-1.png" height="400px" /> ] --- ## Рассмотренные темы 1. Модели пространственных данных 1. Векторные данные (simple features) 1. Растровые данные 1. Системы координат и проецирование Более подробно вопросы выбора проекций и построения сеток координат рассматриваются в следующей главе.