Митя Серебренников, “Основы работы с пространственными данными в R”
План митапа
Типы пространственных данных;
Особенности геоданных;
R и пространственный анализ;
Практика - находим геоданные разных форматов и учимся с ними работать;
Бонус
Что можно почитать и посмотреть по теме.
Загрузим библиотеки
# Функция для загрузки\установки библиотек
package_loader <- function(pkgs) {
if (length(setdiff(pkgs, rownames(installed.packages()))) > 0) {
install.packages(setdiff(pkgs, rownames(installed.packages())))
}
lapply(pkgs, library, character.only = T)
}
# Загружаем библиотеки
package_loader(c('data.table', 'dplyr', 'readr', 'sf', 'geojsonsf', 'geojsonio', 'rnaturalearth', 'osmdata', 'mapview'))## [[1]]
## [1] "data.table" "stats" "graphics" "grDevices" "utils"
## [6] "datasets" "methods" "base"
##
## [[2]]
## [1] "dplyr" "data.table" "stats" "graphics" "grDevices"
## [6] "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "readr" "dplyr" "data.table" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "sf" "readr" "dplyr" "data.table" "stats"
## [6] "graphics" "grDevices" "utils" "datasets" "methods"
## [11] "base"
##
## [[5]]
## [1] "geojsonsf" "sf" "readr" "dplyr" "data.table"
## [6] "stats" "graphics" "grDevices" "utils" "datasets"
## [11] "methods" "base"
##
## [[6]]
## [1] "geojsonio" "geojsonsf" "sf" "readr" "dplyr"
## [6] "data.table" "stats" "graphics" "grDevices" "utils"
## [11] "datasets" "methods" "base"
##
## [[7]]
## [1] "rnaturalearth" "geojsonio" "geojsonsf" "sf"
## [5] "readr" "dplyr" "data.table" "stats"
## [9] "graphics" "grDevices" "utils" "datasets"
## [13] "methods" "base"
##
## [[8]]
## [1] "osmdata" "rnaturalearth" "geojsonio" "geojsonsf"
## [5] "sf" "readr" "dplyr" "data.table"
## [9] "stats" "graphics" "grDevices" "utils"
## [13] "datasets" "methods" "base"
##
## [[9]]
## [1] "mapview" "osmdata" "rnaturalearth" "geojsonio"
## [5] "geojsonsf" "sf" "readr" "dplyr"
## [9] "data.table" "stats" "graphics" "grDevices"
## [13] "utils" "datasets" "methods" "base"1. Типы пространственных данных
а. Векторные данные
Геометрии состоят из точек. Конвенционально x = широта, y = долгота. Кроме этого могут использоваться z = высота, M = мера изменчивости точки;
Виды геометрий:
+++ GEOMETRYCOLLECTION +++
Форматы записи данных:
Well-Known Text (WKT) - Запись точек вектором из геометрий
Well-Known Binary (WKB) - Запись координат в бинарных значениях. Используется в базах данных (т.к. увеличивает скорость работы с данными), но нечитаем и непонятен для человека.
Посмотрим в первом приближении:
# library(rnaturalearth)
world <- ne_countries(returnclass = "sf")
world[1:5, ]## Simple feature collection with 5 features and 94 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
## Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## featurecla scalerank labelrank sovereignt sov_a3
## 0 Admin-0 country 1 6 Fiji FJI
## 1 Admin-0 country 1 3 United Republic of Tanzania TZA
## 2 Admin-0 country 1 7 Western Sahara SAH
## 3 Admin-0 country 1 2 Canada CAN
## 4 Admin-0 country 1 2 United States of America US1
## adm0_dif level type admin adm0_a3 geou_dif
## 0 0 2 Sovereign country Fiji FJI 0
## 1 0 2 Sovereign country United Republic of Tanzania TZA 0
## 2 0 2 Indeterminate Western Sahara SAH 0
## 3 0 2 Sovereign country Canada CAN 0
## 4 1 2 Country United States of America USA 0
## geounit gu_a3 su_dif subunit su_a3 brk_diff
## 0 Fiji FJI 0 Fiji FJI 0
## 1 Tanzania TZA 0 Tanzania TZA 0
## 2 Western Sahara SAH 0 Western Sahara SAH 1
## 3 Canada CAN 0 Canada CAN 0
## 4 United States of America USA 0 United States USA 0
## name name_long brk_a3 brk_name brk_group
## 0 Fiji Fiji FJI Fiji <NA>
## 1 Tanzania Tanzania TZA Tanzania <NA>
## 2 W. Sahara Western Sahara B28 W. Sahara <NA>
## 3 Canada Canada CAN Canada <NA>
## 4 United States of America United States USA United States <NA>
## abbrev postal formal_en formal_fr name_ciawf
## 0 Fiji FJ Republic of Fiji <NA> Fiji
## 1 Tanz. TZ United Republic of Tanzania <NA> Tanzania
## 2 W. Sah. WS Sahrawi Arab Democratic Republic <NA> Western Sahara
## 3 Can. CA Canada <NA> Canada
## 4 U.S.A. US United States of America <NA> United States
## note_adm0 note_brk name_sort name_alt
## 0 <NA> <NA> Fiji <NA>
## 1 <NA> <NA> Tanzania <NA>
## 2 Self admin. Self admin.; Claimed by Morocco Western Sahara <NA>
## 3 <NA> <NA> Canada <NA>
## 4 <NA> <NA> United States of America <NA>
## mapcolor7 mapcolor8 mapcolor9 mapcolor13 pop_est pop_rank gdp_md_est
## 0 5 1 2 2 920938 11 8.374e+03
## 1 3 6 2 2 53950935 16 1.506e+05
## 2 4 7 4 4 603253 11 9.065e+02
## 3 6 6 2 2 35623680 15 1.674e+06
## 4 4 5 1 1 326625791 17 1.856e+07
## pop_year lastcensus gdp_year economy income_grp
## 0 2017 2007 2016 6. Developing region 4. Lower middle income
## 1 2017 2002 2016 7. Least developed region 5. Low income
## 2 2017 NA 2007 7. Least developed region 5. Low income
## 3 2017 2011 2016 1. Developed region: G7 1. High income: OECD
## 4 2017 2010 2016 1. Developed region: G7 1. High income: OECD
## wikipedia fips_10_ iso_a2 iso_a3 iso_a3_eh iso_n3 un_a3 wb_a2 wb_a3 woe_id
## 0 NA FJ FJ FJI FJI 242 242 FJ FJI 23424813
## 1 NA TZ TZ TZA TZA 834 834 TZ TZA 23424973
## 2 NA WI EH ESH ESH 732 732 <NA> <NA> 23424990
## 3 NA CA CA CAN CAN 124 124 CA CAN 23424775
## 4 0 US US USA USA 840 840 US USA 23424977
## woe_id_eh woe_note adm0_a3_is adm0_a3_us adm0_a3_un
## 0 23424813 Exact WOE match as country FJI FJI NA
## 1 23424973 Exact WOE match as country TZA TZA NA
## 2 23424990 Exact WOE match as country MAR SAH NA
## 3 23424775 Exact WOE match as country CAN CAN NA
## 4 23424977 Exact WOE match as country USA USA NA
## adm0_a3_wb continent region_un subregion
## 0 NA Oceania Oceania Melanesia
## 1 NA Africa Africa Eastern Africa
## 2 NA Africa Africa Northern Africa
## 3 NA North America Americas Northern America
## 4 NA North America Americas Northern America
## region_wb name_len long_len abbrev_len tiny homepart
## 0 East Asia & Pacific 4 4 4 NA 1
## 1 Sub-Saharan Africa 8 8 5 NA 1
## 2 Middle East & North Africa 9 14 7 NA 1
## 3 North America 6 6 4 NA 1
## 4 North America 24 13 6 NA 1
## min_zoom min_label max_label ne_id wikidataid name_ar name_bn
## 0 0.0 3.0 8.0 1159320625 Q712 <NA> <NA>
## 1 0.0 3.0 8.0 1159321337 Q924 <NA> <NA>
## 2 4.7 6.0 11.0 1159321223 Q6250 <NA> <NA>
## 3 0.0 1.7 5.7 1159320467 Q16 <NA> <NA>
## 4 0.0 1.7 5.7 1159321369 Q30 <NA> <NA>
## name_de name_en name_es
## 0 Fidschi Fiji Fiyi
## 1 Tansania Tanzania Tanzania
## 2 Westsahara Western Sahara Sahara Occidental
## 3 Kanada Canada Canadá
## 4 Vereinigte Staaten United States of America Estados Unidos
## name_fr name_el name_hi name_hu name_id
## 0 Fidji <NA> <NA> Fidzsi-szigetek Fiji
## 1 Tanzanie <NA> <NA> Tanzánia Tanzania
## 2 Sahara occidental <NA> <NA> Nyugat-Szahara Sahara Barat
## 3 Canada <NA> <NA> Kanada Kanada
## 4 États-Unis <NA> <NA> Amerikai Egyesült Államok Amerika Serikat
## name_it name_ja name_ko name_nl
## 0 Figi <NA> <NA> Fiji
## 1 Tanzania <NA> <NA> Tanzania
## 2 Sahara Occidentale <NA> <NA> Westelijke Sahara
## 3 Canada <NA> <NA> Canada
## 4 Stati Uniti d'America <NA> <NA> Verenigde Staten van Amerika
## name_pl name_pt name_ru name_sv
## 0 Fidzi Fiji <NA> Fiji
## 1 Tanzania Tanzânia <NA> Tanzania
## 2 Sahara Zachodnia Saara Ocidental <NA> Västsahara
## 3 Kanada Canadá <NA> Kanada
## 4 Stany Zjednoczone Estados Unidos <NA> USA
## name_tr name_vi name_zh geometry
## 0 Fiji Fiji <NA> MULTIPOLYGON (((180 -16.067...
## 1 Tanzanya Tanzania <NA> MULTIPOLYGON (((33.90371 -0...
## 2 Bati Sahra Tây Sahara <NA> MULTIPOLYGON (((-8.66559 27...
## 3 Kanada Canada <NA> MULTIPOLYGON (((-122.84 49,...
## 4 Amerika Birlesik Devletleri <NA> <NA> MULTIPOLYGON (((-122.84 49,...# Для иллюстрации
world <- world[1:5, c(4, ncol(world))]
world## Simple feature collection with 5 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
## Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## sovereignt geometry
## 0 Fiji MULTIPOLYGON (((180 -16.067...
## 1 United Republic of Tanzania MULTIPOLYGON (((33.90371 -0...
## 2 Western Sahara MULTIPOLYGON (((-8.66559 27...
## 3 Canada MULTIPOLYGON (((-122.84 49,...
## 4 United States of America MULTIPOLYGON (((-122.84 49,...as.data.frame(world)## sovereignt geometry
## 0 Fiji MULTIPOLYGON (((180 -16.067...
## 1 United Republic of Tanzania MULTIPOLYGON (((33.90371 -0...
## 2 Western Sahara MULTIPOLYGON (((-8.66559 27...
## 3 Canada MULTIPOLYGON (((-122.84 49,...
## 4 United States of America MULTIPOLYGON (((-122.84 49,...b. Растровые данные
Матрица значений пикселей пространственной области. В растрах хронят космические снимки и базы различных геологических данных. Также как и вектор - предназначены для определённого круга задач;
Основной формат: .tif
Основные библиотеки для обработки в R: raster и stars
Сравним:
На этом и последующих занятиях мы будем говорить только о векторных данных.
Главное, нужно помнить, что векторные пространственные данные - это (чаще всего) просто точки на плоскости! Создадим собственные пространственные данные:
p1 = st_point(c(7, 52)) # Рисуем точку 1
p2 = st_point(c(-30, 20)) # Рисуем точку 2
sp_obj = st_sfc(p1, p2, crs = 4326) # Преобразуем в пространственные данные
plot(sp_obj)
# Вы восхитительны!2. Особенности геоданных
Вернёмся к векторным геометриям:
world## Simple feature collection with 5 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
## Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## sovereignt geometry
## 0 Fiji MULTIPOLYGON (((180 -16.067...
## 1 United Republic of Tanzania MULTIPOLYGON (((33.90371 -0...
## 2 Western Sahara MULTIPOLYGON (((-8.66559 27...
## 3 Canada MULTIPOLYGON (((-122.84 49,...
## 4 United States of America MULTIPOLYGON (((-122.84 49,...Что такое CRS?
Coordinate Reference Systems (CRS) / Пространственная привязка и её проекции
Для иллюстрации обратимся к прекрасной работе Тасс.
Картинки чтобы окончательно закрепить идею проекции:
Как посмотреть CRS?
st_crs(world)## Coordinate Reference System:
## User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## wkt:
## BOUNDCRS[
## SOURCECRS[
## GEOGCRS["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]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Geocentric translations (geog2D domain)",
## ID["EPSG",9603]],
## PARAMETER["X-axis translation",0,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",0,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",0,
## ID["EPSG",8607]]]]
На самом деле всё не так сложно.
Пространственный объект можно перевести в другую CRS сразу и быстро, задав через st_transform() нужную проекцию одним из четырёх форматов (просто гуглите то, что нужно для ваших координат):
- Код EPSG (например, “4326” - универсальная проекция для всего мира или “31370” - кастомная проекция для Бельгии):
df <- st_as_sf(df) # <- убедитесь, что ваш объект типа "sf"
df <- st_transform(df,
crs = 4326)- Формула PROJ4 (например, “+proj=longlat +datum=WGS84 +no_defs”):
df <- st_transform(df,
crs = "+proj=longlat +datum=WGS84 +no_defs")- Просто переносом проекции с другого пространственного объекта в формате:
df1 <- st_transform(df1,
crs = st_crs(df2))Примеры
Afg = world[1,]
plot(Afg)
# Изменим crs
Afg <- st_transform(Afg, 2264)
plot(Afg)
# Приведём crs к стандартному формату
Afg <- st_transform(Afg, crs = "+proj=longlat +datum=WGS84 +no_defs")
plot(Afg)
…В геоданных есть ещё множество подводных камней, но этот - ключевой и его необходимо запомнить…
3. R и пространственный анализ
Почему R?
Библиотеки для работы с геоданными в R
(несть им числа, но основные)
- sf - state-of-the-art пространственного анализа в R.
Немного истории:
4. Практика
Что мы будем делать - делать карту мороженного в Израиле! :)
…
Я исхожу из того, что вы знаете dplyr или data.table на базовом уровне, но если у вас вызывает затруднения именно манипуляции с данными - обязательно задавайте вопросы!
Мы отработаем работу с данными на трёх основных форматах для векторных данных: * csv * geojson * shape-file
A. Датасет с мороженками в csv
Загрузим данные. К сожалению, в R есть известная проблема с кодировками и чтобы долго не превращать крокозябры в иврит - уберём колонки названий точек мороженного:
isr_icecream <- read.csv('https://raw.githubusercontent.com/AmitLevinson/Datasets/master/golda/golda_locations.csv')
isr_icecream <- isr_icecream |>
select(-c('city', 'street'))
class(isr_icecream)## [1] "data.frame"Превратим тип данных “data.frame” в “sf”:
isr_icecream <- st_as_sf(
isr_icecream,
dim = "XY",
remove = T,
na.fail = F,
coords = c("lon", "lat"),
crs = "+proj=longlat +datum=WGS84 +no_defs")
isr_icecream## Simple feature collection with 86 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 34.5701 ymin: 29.54952 xmax: 35.55034 ymax: 33.00677
## Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
## number geometry
## 1 03-7797577 POINT (34.86002 32.02051)
## 2 04-646-0705 POINT (34.91794 32.50436)
## 3 08-9428989 POINT (34.95968 29.54952)
## 4 08-9775157 POINT (34.63653 31.79119)
## 5 08-6869144 POINT (34.5701 31.68004)
## 6 08-9555764 POINT (34.82614 31.93624)
## 7 08-6445959 POINT (34.79946 31.23412)
## 8 09-9661614 POINT (34.86808 32.36156)
## 9 03-9441070 POINT (34.73899 32.01818)
## 10 03-7775169 POINT (34.85561 32.08402)plot(isr_icecream[1])
B. Границы Израиля в Shape-file
Идём сюда, ищем Израиль и загружаем архив. Что мы видим внутри?
.shp – файл с геометриями, которые мы подгружаем через st_read()
.dbf
.shx
.prj
Все они важны! По этой причине, когда вам нужно переслать shape-file отправляйте архив из всех четырёх файлов.
Загрузим его и отобразим:
isr_border <- st_read('ISR_adm/ISR_adm0.shp')## Reading layer `ISR_adm0' from data source
## `/builds/upravitelev/r_meetups/content/post/2022-05-14-eu_spb_meetup/ISR_adm/ISR_adm0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 70 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 34.26801 ymin: 29.49708 xmax: 35.90094 ymax: 33.36403
## Geodetic CRS: WGS 84isr_border## Simple feature collection with 1 feature and 70 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 34.26801 ymin: 29.49708 xmax: 35.90094 ymax: 33.36403
## Geodetic CRS: WGS 84
## ID_0 ISO NAME_0 OBJECTID_1 ISO3 NAME_ENGLI NAME_ISO NAME_FAO NAME_LOCAL
## 1 111 ISR Israel 113 ISR Israel ISRAEL Israel Yisra'el
## NAME_OBSOL NAME_VARIA NAME_NONLA NAME_FRENC NAME_SPANI NAME_RUSSI NAME_ARABI
## 1 <NA> <NA> <NA> Israël Israel Израиль إسرائيل
## NAME_CHINE WASPARTOF CONTAINS SOVEREIGN ISO2 WWW FIPS ISON VALIDFR VALIDTO
## 1 以色列 <NA> <NA> Israel IL <NA> IS 376 19480514 Present
## POP2000 SQKM POPSQKM UNREGION1 UNREGION2 DEVELOPING CIS Transition
## 1 6040427 20773.82 290.7711 Western Asia Asia 1 0 0
## OECD WBREGION WBINCOME WBDEBT WBOTHER CEEAC CEMAC
## 1 0 <NA> High income: nonOECD Debt not classified <NA> 0 0
## CEPLG COMESA EAC ECOWAS IGAD IOC MRU SACU UEMOA UMA PALOP PARTA CACM EurAsEC
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## Agadir SAARC ASEAN NAFTA GCC CSN CARICOM EU CAN ACP Landlocked AOSIS SIDS
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## Islands LDC geometry
## 1 0 0 POLYGON ((35.85013 33.2532,...plot(isr_border)
plot(isr_border[1])
plot(st_geometry(isr_border))
C. Точки городов Израиля в GeoJson
Json с геоданными. Достаточно тяжёлый с точки зрения размещения данных, но часто встречающийся формат.
Здесь стоит отвлечься и ознакомиться с “Википедией”” от мира пространственных данных - Open Street Map.
Если это “Википедия”, то с неё можно выкачать данные? Да. По тэгам…
Возьмём тэги городов.
Скачаем GeoJson с городами. Загрузим его. Внутри будет огромное количество колонок, которые сейчас нам не нужны. Оставим первую и последнюю:
isr_city <- st_read('city.geojson')## Reading layer `city' from data source
## `/builds/upravitelev/r_meetups/content/post/2022-05-14-eu_spb_meetup/city.geojson'
## using driver `GeoJSON'
## Simple feature collection with 51 features and 225 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 33.80459 ymin: 30.61197 xmax: 36.56875 ymax: 32.92817
## Geodetic CRS: WGS 84isr_city <- isr_city[,c(1,ncol(isr_city))]
isr_city## Simple feature collection with 51 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 33.80459 ymin: 30.61197 xmax: 36.56875 ymax: 32.92817
## Geodetic CRS: WGS 84
## First 10 features:
## id geometry
## 1 node/29090735 POINT (35.22576 31.77882)
## 2 node/196303399 POINT (34.65299 31.79773)
## 3 node/196307269 POINT (34.78041 32.01931)
## 4 node/196353200 POINT (34.81011 31.96357)
## 5 node/246258161 POINT (35.0089 32.26712)
## 6 node/246849053 POINT (35.17239 32.80458)
## 7 node/270873161 POINT (35.8493 32.55561)
## 8 node/275283306 POINT (35.03665 31.06866)
## 9 node/275547482 POINT (34.7712 31.6094)
## 10 node/275585094 POINT (35.21458 31.26122)plot(isr_city)
Альтернативный способ через пакет osmdata
# library(osmdata)
admin_osm = opq(bbox = 'Sankt-Peterburg') %>% # Попробуйте также 'Sankt-Peterburg' для примера
add_osm_feature(key = "admin_level", value = '5') %>%
osmdata_sf()
# Но есть проблемы с кодировкой
iconv(admin_osm$nodes$tags, from="UTF-8", to="UTF-8") ## character(0)iconv(admin_osm$nodes$tags, from="UTF-8", to="cp1251")## character(0)plot(admin_osm[["osm_multipolygons"]])
Объединяем слои и строим свою первую карту!
# Устанавливаем одинаковые проекции для всех объектов:
isr_border <- st_transform(isr_border, 4326)
isr_city = st_transform(isr_city, st_crs(isr_border))
isr_icecream = st_transform(isr_icecream, st_crs(isr_border))
# Делаем карту (st_geometry() - вывести чистую геометрию, без других колонок или значений в объекте)
{plot(isr_border %>% st_geometry())
plot(isr_city, col = 'red', add = T)
plot(isr_icecream, col = 'blue', add = T)}
Ура :)
Теперь сохраним данные (здесь есть свои подводные камни):
write_sf(isr_border,
"meetup_results_folder/isr_border.shp",
append = T,
layer_options = "ENCODING=UTF-8" # При сохранении данных с не-латинскими character обязательно прописывайте кодировку
)Загрузим в качестве проверки и убедимся, что всё работает
double_isr_border <- st_read("meetup_results_folder/isr_border.shp")## Reading layer `isr_border' from data source
## `/builds/upravitelev/r_meetups/content/post/2022-05-14-eu_spb_meetup/meetup_results_folder/isr_border.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 12 features and 70 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 34.26801 ymin: 29.49708 xmax: 35.90094 ymax: 33.36403
## Geodetic CRS: WGS 84plot(st_geometry(double_isr_border))
5. Бонус
Построим интерактивный график с помощью очень простой по своему синтаксису библиотеки mapview. Туториал по ней можно найти здесь. mapview предназначен для карт небольшого размера.
ЗАПУСКАТЬ ОСТОРОЖНО!!! (на Windows возможны вылеты RStudio - сохраните код перед запуском, чтобы понять, корректно ли установился пакет или нет, если нет - переустановите из исходников)
# library(mapview)
mapview(isr_border) +
mapview(isr_city %>% st_geometry(), col.regions = 'red', legend = FALSE) +
mapview(isr_icecream %>% st_geometry(), col.regions = 'blue', legend = FALSE)6. Где брать данные?
Open street map
naturalearthdata.com/downloads/
download.geofabrik.de/
…и много, много других ресурсов…
7. Что почитать?
Самсонов Т.Е. - Визуализация и анализ географических данных на языке R
Pebesma E., Bivand R. Spatial Data Science. 2020.
Lovelace R., Nowosad J., Muenchow J. Geocomputation with R. 2021
Dorman M. Using R for Spatial Data Analysis. 2021.
…
Можете посмотреть другие материалы из нашей подборки в закрепе Горячей линии