Возможности R для работы с пространственными данными


Мастер-класс на IV Конференции сообщества природоохранных ГИС в России
Национальный парк «Валдайский», 05 октября 2019 г.


Никита Платонов (ИПЭЭ РАН)

Обновлено: 2019-09-28 23:55

Abstract

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

Планирование (до 27 сентября)

Предлагаемый подход к формированию темы мастер-класса

Планируется, что основные материалы для проведения мастер-класса появятся здесь к 27 сентября. Предполагается, что некоторые дополнения могут быть внесены и в более поздний срок, но с учетом того, чтобы эти изменения не потребовали существенного заимствования времени от занятия. Основной режим проведения мастер-класса - это копирование/вставка фрагмента кода, поэтому этот материал может быть использован дистанционными участниками.

Факультативно задумывается и осваивается вариант занятия с использованием Jupyter R Notebook. Установку необходимых программных компонентов нужно выполнить самостоятельно.

Пожелания по затрагиванию какой-либо темы могут быть приняты ✉ не позднее 27 сентября.

Подготовка (до 05 октября)

Операционные системы пользователя

Не для всех операционных систем есть скомпилированные модули (ядро R, библиотеки). В таком случае модули компилируются, и на это уходит какое-то время. Поэтому если ОС не OS Windows, то этот этап нужно пройти заранее.

Установка базового R

Установить или обновить R, например, отсюда. Актуальная версия 3.6.1. Не ниже версии 3.0.

Установка необходимых библиотек

## Loading required namespace: rgdal
## Loading required namespace: sf
## Loading required namespace: raster
## Loading required namespace: mapview
## Loading required namespace: mapedit
## Loading required namespace: jpeg
##     rgdal        sf    raster   ggplot2   leaflet   mapview   mapedit     knitr 
##      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE 
## rmarkdown      jpeg       png      ursa 
##      TRUE      TRUE      TRUE      TRUE

Если отображается TRUE для всех библиотек, то подготовка к занятию осуществлена успешно.

## Everything is ready? 
##                 TRUE

Если где-то выскочило FALSE (например, для библиотеки “foo”), то можно попробовать его установить заново функцией install.packages("foo").

Установка дополнительного программного обеспечения

Pandoc

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

Ссылка на страницу для скачивания. Для пользователей Windows достаточно перейти к скачиванию актуального релиза, и выбрать либо установщик (*.msi), либо архив (*.zip). Запомнить путь, куда произведена установка и где находится файл pandoc.exe и добавить этот путь в переменную окружения %PATH%, например: WindowsKey+Q, ввести “Переменные среды/Environment Variables”, попасть в окошко “Свойства Системы/System Properties”, нажать на кнопку “Переменные среды/Environment Variables” и отредактировать пользовательскую или системную переменную PATH, добавив путь к pandoc.exe.

Чтобы проверить правильно ли установлен Pandoc, в новой R-сессии:

## [1] TRUE

Jupyter R Notebook

Jupyter Notebook для работы с кодом в браузере.

  1. Загрузить и установить менеджер Miniconda (проверено на Windows 64bit Python 3.7).

  2. Запустить conda (На Windows попробовать WinKey+Q и начать вводить “Anaconda”)

  3. Последовательно выполнить команды (в среде Conda, не в среде R):

Команда Описание
conda install jupyter установка Jupyter Notebook
conda install -c r r-recommended r-irkernel Установка R, базовых библиотек и библиотеки для использования R в Jupyter Notebook
conda install -c conda-forge jupytext Установка преобразователя кода R в формат Jupyter Notebook

Предварительно скачайте main.R или только код codeOnly.R, если будете использовать на занятии Jupyter Notebook. Если будет проблема с кодировкой, возьмите codeOnly1251.R Предупреждение: описательный текст (markdown) местами не отформатирован, в частности, не будут отображаться таблицы; также не будут отображаться html-виджеты.

Как использовать материал

Предполагается, что во время занятий будет использоваться интерактивный режим: либо простая среда R (R, запущенный без команд), либо графическая среда R GUI, либо RStudio.

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

К примеру, ниже приведен пример вывода числа π: сверху – команда, снизу – выход.

## [1] 3.141593

Начало работы

Текущий путь, где мы находимся. Здесь появятся файлы, созданные в процессе занятия.

## [1] "D:/RAS/2019/SCGIS"

Его можно поменять через меню RGui/RStudio или с помощью setwd().

Зафиксируем генератор псевдослучайных чисел на определенную последовательность

##  [1]  2  7 10  6  4  1  3  9  5  8

Команда для проверки кириллицы:

## Здесь кириллица? 
##              Да!

Если вывод не читаемый, то возможные пути решения:

  1. Задать кириллицу для символьной локали:
  1. Если скрипты подключаются через source(), то можно настроить запуск R через конфигурационные файлы. Например, скачать .Rprofile, разместить в рабочей директории, перегрузить R. Файл имеет следующую структуру:
local({
   options(encoding="UTF-8")
   Sys.setlocale("LC_CTYPE","Russian")
 })
  1. Если скрипт bar.R с кириллицей в кодировке UTF-8 запускается из командной строки, то использовать следующие параметры запуска.
 R --encoding UTF-8 -f bar.R

Занятие (05 октября)

R как ГИС?

Структуры данных

Числа и имена значений

## [1] 1 2 3 4 5 6 7
##  int [1:7] 1 2 3 4 5 6 7
##  Saturday    Sunday    Monday   Tuesday Wednesday  Thursday    Friday 
##         1         2         3         4         5         6         7
##  Named int [1:7] 1 2 3 4 5 6 7
##  - attr(*, "names")= chr [1:7] "Saturday" "Sunday" "Monday" "Tuesday" ...

Целые числа, числа с плавающей точкой

##  Named int [1:7] 1 2 3 4 5 6 7
##  - attr(*, "names")= chr [1:7] "Saturday" "Sunday" "Monday" "Tuesday" ...
##  Named num [1:7] 1 2 3 4 5 6 7
##  - attr(*, "names")= chr [1:7] "Saturday" "Sunday" "Monday" "Tuesday" ...
## [1] "integer"
## [1] "double"

Логические значения, строки

## [1] FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE
## [1] "logical"
##  logi [1:7] FALSE FALSE TRUE TRUE TRUE FALSE ...
## [1] "Saturday"  "Sunday"    "Monday"    "Tuesday"   "Wednesday" "Thursday" 
## [7] "Friday"
## [1] "character"
##  chr [1:7] "Saturday" "Sunday" "Monday" "Tuesday" "Wednesday" "Thursday" ...

Списки одинаковой длины

## $num
## [1] 1 3 4 5 6 7
## 
## $char
## [1] "Saturday"  "Monday"    "Tuesday"   "Wednesday" "Thursday"  "Friday"
## List of 2
##  $ num : int [1:6] 1 3 4 5 6 7
##  $ char: chr [1:6] "Saturday" "Monday" "Tuesday" "Wednesday" ...
## [1] "list"
## [1] 2
##  num char 
##    6    6

Таблицы

##   num      char
## 1   1  Saturday
## 2   3    Monday
## 3   4   Tuesday
## 4   5 Wednesday
## 5   6  Thursday
## 6   7    Friday
## 'data.frame':    6 obs. of  2 variables:
##  $ num : int  1 3 4 5 6 7
##  $ char: chr  "Saturday" "Monday" "Tuesday" "Wednesday" ...
## [1] "data.frame"
## [1] 6 2

Списки различной длины

## $x
## [1] 7 3 2
## 
## $y
## [1] 4 7 6 5 1
## List of 2
##  $ x: int [1:3] 7 3 2
##  $ y: int [1:5] 4 7 6 5 1
## [1] "list"
## [1] 2
## x y 
## 3 5

Матрицы

##      [,1] [,2] [,3] [,4]
## [1,]    6   23   11   20
## [2,]    2   10   13   21
## [3,]    8    1    4    3

Размерность массива

## [1] 3 4

Структура данных массива

##  int [1:3, 1:4] 6 2 8 23 10 1 11 13 4 20 ...

Массивы

## , , 1
## 
##      [,1] [,2] [,3] [,4]
## [1,]   15   21    3   16
## [2,]   14   11   23   10
## [3,]   18   24    6   22
## 
## , , 2
## 
##      [,1] [,2] [,3] [,4]
## [1,]   19    9    4   20
## [2,]    1   13   17    5
## [3,]    2    8    7   12
##  int [1:3, 1:4, 1:2] 15 14 18 21 11 24 3 23 6 16 ...
## [1] "array"

Факторы

## [1] Monday    Friday    Wednesday Sunday    Tuesday   Thursday  Saturday 
## Levels: Saturday Sunday Monday Tuesday Wednesday Thursday Friday
##  Factor w/ 7 levels "Saturday","Sunday",..: 3 7 5 2 4 6 1
## [1] "factor"
## [1] Tuesday   Saturday  Thursday  Sunday    Monday    Wednesday Friday   
## 7 Levels: Saturday < Sunday < Monday < Tuesday < Wednesday < ... < Friday
##  Ord.factor w/ 7 levels "Saturday"<"Sunday"<..: 4 1 6 2 3 5 7
## [1] "ordered" "factor"

Визуализация данных

Базовые средства R для визуализации (библиотека graphics) пространственных данных:

  • points() – отображение точек (геометрия POINT)

  • lines(), segments(), contour() – отображение линий (геометрия LINESTRING)

  • polygon(), polypath() – отображение полигонов (геометрия POLYGON)

  • image(), rasterImage() – растровые изображения

  • text(), legend() – аннотациии

  • axis(), mtext() – рамочное оформление

data(volcano)

[]{style="display: none;"}

Это вулкан Maunga Whau (Mt Eden).

[]{style="display: none;"}

Манипуляции с файлами пространственных данных

Векторные данные

Библиотека Формат Импорт Экспорт
sp GDAL vector drivers rgdal::readOGR() rgdal::writeOGR()
sf GDAL vector drivers st_read() st_write()

Растровые данные

Библиотека Формат Импорт Экспорт
sp GDAL raster drivers rgdal::readGDAL() rgdal::writeGDAL()
raster GDAL raster drivers raster(), brick(), stack() writeRaster()
ncdf4 NetCDF nc_open() nc_create()
ursa ENVI read_envi() write_envi()

Особенности форматов данных

Растровые данные

  • Для хранения многослойные растров используется BSQ/BIL/BIP чередование слоев/строк/пикелей. Самый неэффективный - это BIP. При пространственно-временном анализе можно выбрать BIL, для большинства случаев - BSQ.
  • Целочисленный GeoTIFF быстро пишется и читается при использовании функций из библиотеки rgdal

Векторные данные

  • Хорошую скорость чтения и записи демонстрирует формат “SQLite” при использовании библиотеки sf.

  • “GeoJSON” не очень быстрый.

  • При записи использовать опции, предусмотренные для выбранного формата данных

  • При записи “ESRI Shapefile” обращать внимания на *.prj, так как у QGIS и ESRI-продуктов разные восприятия файлов проекций.

sf или sp?

  • В пользу sf больше аргументов:

    • удобнее, активно развивается, поддерживается

    • в sf объекты класса S3, в sp объекты класса S4 (строже, но медленнее)

    • Эффективность с геометрией POINT выше у sp из-за представления атрибутивной таблицы и геометрии в одной таблице.

    • Ряд библиотек завязаны на формат данных sp, например adehabitatHR.

    • есть sf::as_Spatial() для преобразования в объекты sp.

Импорт пространственных данных

Векторные данные – rgdal/sp

## Source: "C:\Software\Rtools\library\rgdal\vectors\scot_BNG.shp", layer: "scot_BNG"
## Driver: ESRI Shapefile; number of rows: 56 
## Feature type: wkbPolygon with 2 dimensions
## Extent: (7094.552 529495) - (468285.5 1218342)
## CRS: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs  
## LDID: 0 
## Number of fields: 13 
##      name type length typeName
## 1   SP_ID    4      5   String
## 2    NAME    4     13   String
## 3    ID_x    2     19     Real
## 4   COUNT    2     19     Real
## 5     SMR    2     19     Real
## 6    LONG    2     19     Real
## 7     LAT    2     19     Real
## 8      PY    2     19     Real
## 9    EXP_    2     19     Real
## 10    AFF    2     19     Real
## 11 X_COOR    2     19     Real
## 12 Y_COOR    2     19     Real
## 13   ID_y    2     19     Real
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Software\Rtools\library\rgdal\vectors\scot_BNG.shp", layer: "scot_BNG"
## with 56 features
## It has 13 fields
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  2 obs. of  13 variables:
##   .. ..$ SP_ID : chr [1:2] "0" "1"
##   .. ..$ NAME  : chr [1:2] "Sutherland" "Nairn"
##   .. ..$ ID_x  : num [1:2] 12 13
##   .. ..$ COUNT : num [1:2] 5 3
##   .. ..$ SMR   : num [1:2] 279 278
##   .. ..$ LONG  : num [1:2] 58.1 57.5
##   .. ..$ LAT   : num [1:2] 4.64 3.98
##   .. ..$ PY    : num [1:2] 37521 29374
##   .. ..$ EXP_  : num [1:2] 1.8 1.1
##   .. ..$ AFF   : num [1:2] 16 10
##   .. ..$ X_COOR: num [1:2] 247757 289592
##   .. ..$ Y_COOR: num [1:2] 924954 842305
##   .. ..$ ID_y  : num [1:2] 12 13
##   ..@ polygons   :List of 2
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. .. .. ..@ Polygons :List of 1
##   .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
##   .. .. .. .. .. .. ..@ labpt  : num [1:2] 247754 924885
##   .. .. .. .. .. .. ..@ area   : num 3.92e+09
##   .. .. .. .. .. .. ..@ hole   : logi FALSE
##   .. .. .. .. .. .. ..@ ringDir: int 1
##   .. .. .. .. .. .. ..@ coords : num [1:73, 1:2] 254218 254359 252991 244974 259133 ...
##   .. .. .. ..@ plotOrder: int 1
##   .. .. .. ..@ labpt    : num [1:2] 247754 924885
##   .. .. .. ..@ ID       : chr "0"
##   .. .. .. ..@ area     : num 3.92e+09
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. .. .. ..@ Polygons :List of 1
##   .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
##   .. .. .. .. .. .. ..@ labpt  : num [1:2] 289590 842235
##   .. .. .. .. .. .. ..@ area   : num 4.7e+08
##   .. .. .. .. .. .. ..@ hole   : logi FALSE
##   .. .. .. .. .. .. ..@ ringDir: int 1
##   .. .. .. .. .. .. ..@ coords : num [1:9, 1:2] 292542 301399 301520 289304 288882 ...
##   .. .. .. ..@ plotOrder: int 1
##   .. .. .. ..@ labpt    : num [1:2] 289590 842235
##   .. .. .. ..@ ID       : chr "1"
##   .. .. .. ..@ area     : num 4.7e+08
##   ..@ plotOrder  : int [1:2] 1 2
##   ..@ bbox       : num [1:2, 1:2] 204670 825310 306671 974566
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps"| __truncated__
## [1] TRUE
## [1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

Векторные данные – sf

## Reading layer `scot_BNG' from data source `C:\Software\Rtools\library\rgdal\vectors\scot_BNG.shp' using driver `ESRI Shapefile'
## Simple feature collection with 56 features and 13 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 7094.552 ymin: 529495 xmax: 468285.5 ymax: 1218342
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs
## Classes 'sf' and 'data.frame':   56 obs. of  14 variables:
##  $ SP_ID   : chr  "0" "1" "2" "3" ...
##  $ NAME    : chr  "Sutherland" "Nairn" "Inverness" "Banff-Buchan" ...
##  $ ID_x    : num  12 13 19 2 17 16 21 50 15 25 ...
##  $ COUNT   : num  5 3 9 39 2 9 16 6 17 19 ...
##  $ SMR     : num  279 278 163 450 187 ...
##  $ LONG    : num  58.1 57.5 57.2 57.6 57.1 ...
##  $ LAT     : num  4.64 3.98 4.73 2.36 4.09 3 2.98 3.2 3.1 3.3 ...
##  $ PY      : num  37521 29374 162867 231337 27075 ...
##  $ EXP_    : num  1.8 1.1 5.5 8.7 1.1 4.6 10.5 19.6 7.8 15.5 ...
##  $ AFF     : num  16 10 7 16 10 16 7 1 7 1 ...
##  $ X_COOR  : num  247757 289592 249685 385776 280492 ...
##  $ Y_COOR  : num  924954 842305 825855 852378 801036 ...
##  $ ID_y    : num  12 13 19 2 17 16 21 50 15 25 ...
##  $ geometry:sfc_MULTIPOLYGON of length 56; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:73, 1:2] 254218 254359 252991 244974 259133 ...
##   ..- attr(*, "class")= chr  "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  "SP_ID" "NAME" "ID_x" "COUNT" ...

Растровые данные – rgdal/sp

Получение данных напрямую

## C:/Software/Rtools/library/rgdal/pictures/cea.tif has GDAL driver GTiff 
## and has 515 rows and 514 columns
## Formal class 'SpatialGridDataFrame' [package "sp"] with 4 slots
##   ..@ data       :'data.frame':  264710 obs. of  1 variable:
##   .. ..$ band1: int [1:264710] 0 0 0 0 0 0 0 0 0 0 ...
##   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
##   .. .. ..@ cellcentre.offset: Named num [1:2] -28463 4225003
##   .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
##   .. .. ..@ cellsize         : num [1:2] 60 60
##   .. .. ..@ cells.dim        : int [1:2] 514 515
##   ..@ bbox       : num [1:2, 1:2] -28493 4224973 2358 4255885
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=cea +lon_0=-117.333333333333 +lat_ts=33.75 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs +ellps=clrk66 +na"| __truncated__
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0    58.0    99.0   103.1   140.0   255.0

Получение данных более гибким способом

## Warning in rgdal::GDALinfo(tifname): statistics not supported by this driver
##          rows       columns         bands          ll.x          ll.y 
##     515.00000     514.00000       1.00000  -28493.16678 4224973.14326 
##         res.x         res.y     oblique.x     oblique.y 
##      60.02214      60.02214       0.00000       0.00000
##  int [1:514, 1:515] 0 0 0 0 0 0 0 0 0 0 ...
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0    58.0    99.0   103.1   140.0   255.0

Растровые данные – raster

## class      : RasterBrick 
## dimensions : 515, 514, 264710, 1  (nrow, ncol, ncell, nlayers)
## resolution : 60.02214, 60.02214  (x, y)
## extent     : -28493.17, 2358.212, 4224973, 4255885  (xmin, xmax, ymin, ymax)
## crs        : +proj=cea +lon_0=-117.333333333333 +lat_ts=33.75 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat 
## source     : C:/Software/Rtools/library/rgdal/pictures/cea.tif 
## names      : cea 
## min values :   0 
## max values : 255
##      d3      v3 
##    8256 1059160
##  int [1:264710, 1] 0 0 0 0 0 0 0 0 0 0 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr "cea"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0    58.0    99.0   103.1   140.0   255.0

Характеристики пространственных данных

Атрибутивная таблица

sp:

##   SP_ID         NAME ID_x COUNT   SMR  LONG  LAT     PY EXP_ AFF   X_COOR
## 0     0   Sutherland   12     5 279.3 58.06 4.64  37521  1.8  16 247756.6
## 1     1        Nairn   13     3 277.8 57.47 3.98  29374  1.1  10 289592.2
## 2     2    Inverness   19     9 162.7 57.24 4.73 162867  5.5   7 249685.5
## 3     3 Banff-Buchan    2    39 450.3 57.56 2.36 231337  8.7  16 385776.1
## 4     4     Bedenoch   17     2 186.9 57.06 4.09  27075  1.1  10 280491.6
## 5     5   Kincardine   16     9 197.8 57.00 3.00 111665  4.6  16 351350.5
##     Y_COOR ID_y
## 0 924953.6   12
## 1 842305.2   13
## 2 825854.5   19
## 3 852378.2    2
## 4 801035.6   17
## 5 792161.7   16

sf:

##   SP_ID         NAME ID_x COUNT   SMR  LONG  LAT     PY EXP_ AFF   X_COOR
## 1     0   Sutherland   12     5 279.3 58.06 4.64  37521  1.8  16 247756.6
## 2     1        Nairn   13     3 277.8 57.47 3.98  29374  1.1  10 289592.2
## 3     2    Inverness   19     9 162.7 57.24 4.73 162867  5.5   7 249685.5
## 4     3 Banff-Buchan    2    39 450.3 57.56 2.36 231337  8.7  16 385776.1
## 5     4     Bedenoch   17     2 186.9 57.06 4.09  27075  1.1  10 280491.6
## 6     5   Kincardine   16     9 197.8 57.00 3.00 111665  4.6  16 351350.5
##     Y_COOR ID_y
## 1 924953.6   12
## 2 842305.2   13
## 3 825854.5   19
## 4 852378.2    2
## 5 801035.6   17
## 6 792161.7   16

Геометрия

sp:

## Formal class 'SpatialPolygons' [package "sp"] with 4 slots
##   ..@ polygons   :List of 2
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. .. .. ..@ Polygons :List of 1
##   .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
##   .. .. .. .. .. .. ..@ labpt  : num [1:2] 247754 924885
##   .. .. .. .. .. .. ..@ area   : num 3.92e+09
##   .. .. .. .. .. .. ..@ hole   : logi FALSE
##   .. .. .. .. .. .. ..@ ringDir: int 1
##   .. .. .. .. .. .. ..@ coords : num [1:73, 1:2] 254218 254359 252991 244974 259133 ...
##   .. .. .. ..@ plotOrder: int 1
##   .. .. .. ..@ labpt    : num [1:2] 247754 924885
##   .. .. .. ..@ ID       : chr "0"
##   .. .. .. ..@ area     : num 3.92e+09
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. .. .. ..@ Polygons :List of 1
##   .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
##   .. .. .. .. .. .. ..@ labpt  : num [1:2] 289590 842235
##   .. .. .. .. .. .. ..@ area   : num 4.7e+08
##   .. .. .. .. .. .. ..@ hole   : logi FALSE
##   .. .. .. .. .. .. ..@ ringDir: int 1
##   .. .. .. .. .. .. ..@ coords : num [1:9, 1:2] 292542 301399 301520 289304 288882 ...
##   .. .. .. ..@ plotOrder: int 1
##   .. .. .. ..@ labpt    : num [1:2] 289590 842235
##   .. .. .. ..@ ID       : chr "1"
##   .. .. .. ..@ area     : num 4.7e+08
##   ..@ plotOrder  : int [1:2] 1 2
##   ..@ bbox       : num [1:2, 1:2] 204670 825310 306671 974566
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps"| __truncated__

sf:

## Geometry set for 2 features 
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 204669.7 ymin: 825310 xmax: 306671.2 ymax: 974566.2
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs
## MULTIPOLYGON (((254218.3 967027, 254358.6 96696...
## MULTIPOLYGON (((292542.2 859596.3, 301398.8 836...
## sfc_MULTIPOLYGON of length 56; first list element: List of 1
##  $ :List of 1
##   ..$ : num [1:73, 1:2] 254218 254359 252991 244974 259133 ...
##  - attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"

Проекция

В R данные проекции представляются в виде PROJ.4, поэтому при импорте/экспорте данных осуществляется двойное WKT -> PROJ4 -> WKT преобразование

sp:

## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894"

sf:

## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs"

Пространственный охват

sp:

##          min       max
## x   7094.552  468285.5
## y 529495.039 1218342.5

sf:

##        xmin        ymin        xmax        ymax 
##    7094.552  529495.039  468285.495 1218342.493

Создание пространственных данных

Выйдем из здания вокзала ст. Валдай и создадим точечный Spatial-объект:

[]{style="display: none;"}

Смотрим на структуру данных:

## Formal class 'SpatialPoints' [package "sp"] with 3 slots
##   ..@ coords     : num [1, 1:2] 33.2 58
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr "1"
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   ..@ bbox       : num [1:2, 1:2] 33.2 58 33.2 58
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Изменим проекцию для расчета расстояний по координатам:

Извлечем координаты:

Решаем перемещаться минутными сегментами в течение часа:

На основе координат создаем таблицу перемещения. По умолчанию весь период стоим на месте и смотрим на север:

Скорость перемещений (в м мин-1) в течение минуты случайна:

##  num [1:60] 33.97 13.93 23.8 9.04 21.29 ...

Задаем, что если следующий сегмент длинный, то отклонение направления от предыдущего сегмента будет меньше:

##  num [1:60] 0.978 2.263 0.686 -2.839 2.318 ...

Заполняем последующие шаги на основе предыдущих:

##   step     look         x         y
## 1    0 1.570796 -103834.6 -1982.122
## 2    1 2.548995 -103862.8 -1963.148
## 3    2 4.811903 -103861.4 -1977.006
## 4    3 5.497473 -103844.6 -1993.837
## 5    4 2.657984 -103852.6 -1989.634
## 6    5 4.976429 -103847.0 -2010.182

Создадим линейный sf-объект

Матрица (с двумя столбцами) координат задается как LINESTRING:

## List of 60
##  $ : 'XY' num [1:2, 1:2] -103835 -103863 -1982 -1963
##  $ : 'XY' num [1:2, 1:2] -103863 -103861 -1963 -1977
##  $ : 'XY' num [1:2, 1:2] -103861 -103845 -1977 -1994
##  $ : 'XY' num [1:2, 1:2] -103845 -103853 -1994 -1990
##  $ : 'XY' num [1:2, 1:2] -103853 -103847 -1990 -2010
##  $ : 'XY' num [1:2, 1:2] -103847 -103808 -2010 -2067
##  $ : 'XY' num [1:2, 1:2] -103808 -103748 -2067 -2074
##  $ : 'XY' num [1:2, 1:2] -103748 -103733 -2074 -2098
##  $ : 'XY' num [1:2, 1:2] -103733 -103769 -2098 -2129
##  $ : 'XY' num [1:2, 1:2] -103769 -103803 -2129 -2149
##  $ : 'XY' num [1:2, 1:2] -103803 -103797 -2149 -2154
##  $ : 'XY' num [1:2, 1:2] -103797 -103737 -2154 -2160
##  $ : 'XY' num [1:2, 1:2] -103737 -103700 -2160 -2164
##  $ : 'XY' num [1:2, 1:2] -103700 -103696 -2164 -2116
##  $ : 'XY' num [1:2, 1:2] -103696 -103701 -2116 -2120
##  $ : 'XY' num [1:2, 1:2] -103701 -103719 -2120 -2191
##  $ : 'XY' num [1:2, 1:2] -103719 -103732 -2191 -2179
##  $ : 'XY' num [1:2, 1:2] -103732 -103785 -2179 -2125
##  $ : 'XY' num [1:2, 1:2] -103785 -103792 -2125 -2145
##  $ : 'XY' num [1:2, 1:2] -103792 -103827 -2145 -2213
##  $ : 'XY' num [1:2, 1:2] -103827 -103809 -2213 -2224
##  $ : 'XY' num [1:2, 1:2] -103809 -103760 -2224 -2248
##  $ : 'XY' num [1:2, 1:2] -103760 -103699 -2248 -2229
##  $ : 'XY' num [1:2, 1:2] -103699 -103682 -2229 -2169
##  $ : 'XY' num [1:2, 1:2] -103682 -103642 -2169 -2133
##  $ : 'XY' num [1:2, 1:2] -103642 -103657 -2133 -2110
##  $ : 'XY' num [1:2, 1:2] -103657 -103639 -2110 -2094
##  $ : 'XY' num [1:2, 1:2] -103639 -103598 -2094 -2086
##  $ : 'XY' num [1:2, 1:2] -103598 -103565 -2086 -2023
##  $ : 'XY' num [1:2, 1:2] -103565 -103529 -2023 -2016
##  $ : 'XY' num [1:2, 1:2] -103529 -103536 -2016 -2009
##  $ : 'XY' num [1:2, 1:2] -103536 -103520 -2009 -1977
##  $ : 'XY' num [1:2, 1:2] -103520 -103502 -1977 -1928
##  $ : 'XY' num [1:2, 1:2] -103502 -103464 -1928 -1906
##  $ : 'XY' num [1:2, 1:2] -103464 -103463 -1906 -1847
##  $ : 'XY' num [1:2, 1:2] -103463 -103434 -1847 -1859
##  $ : 'XY' num [1:2, 1:2] -103434 -103403 -1859 -1916
##  $ : 'XY' num [1:2, 1:2] -103403 -103345 -1916 -1949
##  $ : 'XY' num [1:2, 1:2] -103345 -103291 -1949 -1929
##  $ : 'XY' num [1:2, 1:2] -103291 -103232 -1929 -1964
##  $ : 'XY' num [1:2, 1:2] -103232 -103194 -1964 -2010
##  $ : 'XY' num [1:2, 1:2] -103194 -103186 -2010 -2033
##  $ : 'XY' num [1:2, 1:2] -103186 -103192 -2033 -2034
##  $ : 'XY' num [1:2, 1:2] -103192 -103201 -2034 -2017
##  $ : 'XY' num [1:2, 1:2] -103201 -103244 -2017 -1982
##  $ : 'XY' num [1:2, 1:2] -103244 -103230 -1982 -1931
##  $ : 'XY' num [1:2, 1:2] -103230 -103225 -1931 -1923
##  $ : 'XY' num [1:2, 1:2] -103225 -103210 -1923 -1927
##  $ : 'XY' num [1:2, 1:2] -103210 -103182 -1927 -1882
##  $ : 'XY' num [1:2, 1:2] -103182 -103176 -1882 -1901
##  $ : 'XY' num [1:2, 1:2] -103176 -103162 -1901 -1927
##  $ : 'XY' num [1:2, 1:2] -103162 -103140 -1927 -1973
##  $ : 'XY' num [1:2, 1:2] -103140 -103168 -1973 -2032
##  $ : 'XY' num [1:2, 1:2] -103168 -103206 -2032 -2054
##  $ : 'XY' num [1:2, 1:2] -103206 -103281 -2054 -2080
##  $ : 'XY' num [1:2, 1:2] -103281 -103300 -2080 -2078
##  $ : 'XY' num [1:2, 1:2] -103300 -103369 -2078 -2074
##  $ : 'XY' num [1:2, 1:2] -103369 -103437 -2074 -2040
##  $ : 'XY' num [1:2, 1:2] -103437 -103507 -2040 -2064
##  $ : 'XY' num [1:2, 1:2] -103507 -103527 -2064 -2068

Набор сегментов объединяем в геометрию:

## sfc_LINESTRING of length 60; first list element:  'XY' num [1:2, 1:2] -103835 -103863 -1982 -1963

С геометрией связываем атрибутивную таблицу.

## Classes 'sf' and 'data.frame':   60 obs. of  3 variables:
##  $ step    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ segment : num  33.97 13.93 23.8 9.04 21.29 ...
##  $ geometry:sfc_LINESTRING of length 60; first list element:  'XY' num [1:2, 1:2] -103835 -103863 -1982 -1963
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
##   ..- attr(*, "names")= chr  "step" "segment"

Статическая визуализация

[]{style="display: none;"}

Графическое отображение переопределенной под класс объекта функцией plot() средствами библиотеки sp(факультативно) и sf развито слабо…

[]{style="display: none;"}

Возвращаясь к ранее загруженным данным:

[]{style="display: none;"}

[]{style="display: none;"}

… Поэтому используются возможности других библиотек.

[]{style="display: none;"}

Интерактивная визуализация

Экспорт пространственных данных

Проверим, появился ли файл:

## [1] "afterTrain.geojson"
[]{style="display: none;"} 

## [1] "scotland.sqlite"
## Note: unable to make bold caption for "категория"
[]{style="display: none;"} 

## [1] "track.dbf" "track.prj" "track.shp" "track.shx"
[]{style="display: none;"}

Рисование

Этот раздел предлагается пройти самостоятельно

Воспроизводимые вычисления и публикация результата

Ниже приведены фрагмента кода, которые не были выполнены при составлении программы занятия. Их предлагается выполнить самостоятельно. Также необходимо установленный Pandoc.

Содержимое занятия сгенерировано из файла main.R. Загрузим его, переименовав:

На основе этого файла создадим markdown-документ:

Затем из markdown-документа сформируем html-документ:

Полученный html-документ можно открыть в браузере.

И можно скопировать на флешку на память:))

Успехов!

Дополнительная информация

Как узнать об R поглубже

  • R-Bloggers – современные тенденции R

  • R’s Spatial and SpatioTemporal Task Views

  • R-Spatial – веб-сайт и блог для интересующихся использованием R для анализа пространственных и пространственно-временных данных

  • Stackoverflow с тегом #R.

  • Stackexchange о ГИС, в т.ч. с использованием R.

  • Package’s vignettes – обобщенное знакомство с библиотекой. Обычно содержат воспроизводимый код.

Ведущий

Данные после занятия остались на диске. Если не нужны, выполнить следующее: