class: center, middle, inverse, title-slide # GIS and mapping in R ## Introduction to the sf package ### Olivier Gimenez ### 2022-05-13 --- # Simple Features for R: the sf package <img src="img/sfcartoon.jpg" width="700px" style="display: block; margin: auto;" /> --- class: inverse, middle, center # Introduction --- # What's so nice about `sf`? * Easy to work with spatial data because the distinction between spatial data and other forms of data is minimized * Spatial objects are stored as **dataframes**, with the feature geometries stored in list-columns * This is similar to the way that **spatial databases** are structured * All functions begin with `st_` for easy autofill with RStudio tab * Functions are **pipe-friendly** * `dplyr` and `tidyr` verbs have been defined for the `sf` objects * `ggplot2` is able to plot `sf` objects directly <img src="img/animatedsf.gif" width="100px" style="display: block; margin: auto;" /> --- # Load packages ```r library(sf) # GIS package ``` ``` ## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1 ``` ```r library(tidyverse) # tidyverse packages, dplyr and ggplot2 among others ``` ``` ## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ── ``` ``` ## ✔ ggplot2 3.3.5 ✔ purrr 0.3.4.9000 ## ✔ tibble 3.1.6 ✔ dplyr 1.0.7 ## ✔ tidyr 1.1.4 ✔ stringr 1.4.0 ## ✔ readr 2.0.0 ✔ forcats 0.5.1 ``` ``` ## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ── ## ✖ dplyr::filter() masks stats::filter() ## ✖ dplyr::lag() masks stats::lag() ``` ```r theme_set(theme_minimal(base_size = 14)) # set ggplot theme ``` --- ## Vector layers in `sf` * The `sf` class is a hierarchical structure composed of 3 classes * In green, **sf** - Vector layer object, `data.frame` with `\(\geq 1\)` attribute columns and 1 geometry column * In red, **sfc** - Geometric part of vector layer - geometry column * In blue, **sfg** - Geometry of individual [simple feature](https://en.wikipedia.org/wiki/Simple_Features) <img src="img/sf_xfig.png" width="1497" style="display: block; margin: auto;" /> --- ## Simple feature geometry **`sfg`** <img src="index_files/figure-html/unnamed-chunk-5-1.png" width="550px" style="display: block; margin: auto;" /> --- class: inverse, middle, center # First steps --- # Case study <img src="img/oryxpaper.png" width="600px" style="display: block; margin: auto;" /> --- # Read in spatial data ```r studysites_raw <- st_read("shp/bearpyrenees.shp") ``` ``` ## Reading layer `bearpyrenees' from data source ## `/Users/oliviergimenez/Dropbox/OG/GITHUB/intro_spatialR/shp/bearpyrenees.shp' ## using driver `ESRI Shapefile' ## Simple feature collection with 138 features and 4 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 315722.7 ymin: 1704775 xmax: 644368 ymax: 1800721 ## Projected CRS: Lambert_Conformal_Conic ``` --- # Examine structure ```r glimpse(studysites_raw) ``` ``` ## Rows: 138 ## Columns: 5 ## $ Numero <dbl> 640101, 640102, 640203, 640204, 640202, 640201, 640301, 640… ## $ pres_08_12 <chr> "Absence", "Absence", "Occasionnelle", "Occasionnelle", "Ab… ## $ Perimeter <dbl> 53446.71, 60814.07, 38908.05, 32749.40, 36869.03, 37629.04,… ## $ Iti2012 <int> NA, NA, 1, 1, NA, NA, NA, 1, 1, 1, 1, NA, 0, 1, 1, 1, 1, 1,… ## $ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((315785.1 17..., MULTIPOLYGON (… ``` --- # Examine structure ```r studysites_raw ``` ``` ## Simple feature collection with 138 features and 4 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 315722.7 ymin: 1704775 xmax: 644368 ymax: 1800721 ## Projected CRS: Lambert_Conformal_Conic ## First 10 features: ## Numero pres_08_12 Perimeter Iti2012 geometry ## 1 640101 Absence 53446.71 NA MULTIPOLYGON (((315785.1 17... ## 2 640102 Absence 60814.07 NA MULTIPOLYGON (((326069.1 17... ## 3 640203 Occasionnelle 38908.05 1 MULTIPOLYGON (((328868.9 17... ## 4 640204 Occasionnelle 32749.40 1 MULTIPOLYGON (((338223.2 17... ## 5 640202 Absence 36869.03 NA MULTIPOLYGON (((348261.2 17... ## 6 640201 Absence 37629.04 NA MULTIPOLYGON (((347345.4 17... ## 7 640301 Absence 29586.76 NA MULTIPOLYGON (((350581.2 17... ## 8 640302 Absence 25304.44 1 MULTIPOLYGON (((353106.3 17... ## 9 640401 Occasionnelle 44181.28 1 MULTIPOLYGON (((353084 1783... ## 10 640402 Reguliere 43850.75 1 MULTIPOLYGON (((359081.4 17... ``` --- # Select relevant columns ```r studysites <- studysites_raw %>% * select('bearpresence' = pres_08_12, * 'idsite' = Numero) studysites ``` ``` ## Simple feature collection with 138 features and 2 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 315722.7 ymin: 1704775 xmax: 644368 ymax: 1800721 ## Projected CRS: Lambert_Conformal_Conic ## First 10 features: ## bearpresence idsite geometry ## 1 Absence 640101 MULTIPOLYGON (((315785.1 17... ## 2 Absence 640102 MULTIPOLYGON (((326069.1 17... ## 3 Occasionnelle 640203 MULTIPOLYGON (((328868.9 17... ## 4 Occasionnelle 640204 MULTIPOLYGON (((338223.2 17... ## 5 Absence 640202 MULTIPOLYGON (((348261.2 17... ## 6 Absence 640201 MULTIPOLYGON (((347345.4 17... ## 7 Absence 640301 MULTIPOLYGON (((350581.2 17... ## 8 Absence 640302 MULTIPOLYGON (((353106.3 17... ## 9 Occasionnelle 640401 MULTIPOLYGON (((353084 1783... ## 10 Reguliere 640402 MULTIPOLYGON (((359081.4 17... ``` --- # Our first map ```r studysites %>% * ggplot() + * geom_sf() + labs(title = 'Brown bear monitoring sites in the French Pyrenees mountains', subtitle = 'Source: French Biodiversity Agency') ``` <img src="index_files/figure-html/unnamed-chunk-11-1.png" style="display: block; margin: auto;" /> --- # Where did the species occur? ```r studysites %>% ggplot() + * geom_sf(aes(fill = bearpresence)) + labs(title = 'Brown bear presence in the French Pyrenees mountains', subtitle = 'Source: French Biodiversity Agency') ``` <img src="index_files/figure-html/unnamed-chunk-12-1.png" style="display: block; margin: auto;" /> --- # Where did the species occur? ```r studysites %>% ggplot() + geom_sf(aes(fill = bearpresence)) + labs(title = 'Brown bear presence in the French Pyrenees mountains', subtitle = 'Source: French Biodiversity Agency', * fill = "Presence") ``` <img src="index_files/figure-html/unnamed-chunk-13-1.png" style="display: block; margin: auto;" /> --- # Take control of your legends ```r studysites %>% ggplot() + geom_sf(aes(fill = bearpresence)) + labs(title = 'Brown bear presence in the French Pyrenees mountains', subtitle = 'Source: French Biodiversity Agency') + * scale_fill_manual(values = c('gray90','steelblue1','steelblue4'), * name = "Bear presence", * labels = c("Absent", "Occasional", "Regular")) ``` <img src="index_files/figure-html/unnamed-chunk-14-1.png" style="display: block; margin: auto;" /> --- class: inverse, middle, center # Spatial operations: transform, crop, intersect, join --- # Forest cover * Forest cover might be a driver of brown bear distribution * We use corine land cover (CLC) data (2012 version) to get forest cover * Data can be downloaded [from the official website](http://www.donnees.statistiques.developpement-durable.gouv.fr/donneesCLC/CLC/millesime/CLC12_FR_RGF_SHP.zip) * An explanation of what's in the data is available [here](https://www.statistiques.developpement-durable.gouv.fr/sites/default/files/2019-07/Nomenclature_CLC.pdf). --- # Read in data ```r clc2012 <- st_read("shp/CLC12_FR_RGF.shp") ``` ``` ## Reading layer `CLC12_FR_RGF' from data source ## `/Users/oliviergimenez/Dropbox/OG/GITHUB/intro_spatialR/shp/CLC12_FR_RGF.shp' ## using driver `ESRI Shapefile' ## Simple feature collection with 274573 features and 3 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: 73767.79 ymin: 6021170 xmax: 1267595 ymax: 7133490 ## Projected CRS: RGF93 / Lambert-93 ``` --- # Read in data ```r clc2012 ``` ``` ## Simple feature collection with 274573 features and 3 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: 73767.79 ymin: 6021170 xmax: 1267595 ymax: 7133490 ## Projected CRS: RGF93 / Lambert-93 ## First 10 features: ## ID CODE_12 AREA_HA geometry ## 1 FR-274573 523 4.879872e+06 POLYGON ((657640.6 7133490,... ## 2 FR-271961 423 8.555978e+02 POLYGON ((669050 7111012, 6... ## 3 FR-265902 331 3.092622e+02 POLYGON ((669132.3 7110834,... ## 4 FR-250624 322 7.593459e+01 POLYGON ((669181.5 7110656,... ## 5 FR-250623 322 1.826213e+02 POLYGON ((666715.4 7109580,... ## 6 FR-25381 112 3.282088e+02 POLYGON ((666715.4 7109580,... ## 7 FR-265901 331 2.503051e+01 POLYGON ((667326.7 7108770,... ## 8 FR-147959 242 6.157960e+01 POLYGON ((669097.7 7109171,... ## 9 FR-270463 333 2.949202e+01 POLYGON ((666773.8 7109345,... ## 10 FR-62513 211 9.808718e+05 POLYGON ((669737.1 7108874,... ``` --- # Extract forest codes ```r forest <- clc2012 %>% * filter(CODE_12 == '311' | CODE_12 == '312' | CODE_12 == '313') forest ``` ``` ## Simple feature collection with 69111 features and 3 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: 125951.9 ymin: 6021917 xmax: 1243919 ymax: 7100648 ## Projected CRS: RGF93 / Lambert-93 ## First 10 features: ## ID CODE_12 AREA_HA geometry ## 1 FR-211653 311 114.81939 POLYGON ((657867.2 7100629,... ## 2 FR-227070 312 29.01440 POLYGON ((632270.7 7100616,... ## 3 FR-211652 311 36.78086 POLYGON ((630186.3 7100313,... ## 4 FR-211651 311 90.71384 POLYGON ((625545 7098175, 6... ## 5 FR-211650 311 32.23925 POLYGON ((660260.8 7097323,... ## 6 FR-211649 311 43.00924 POLYGON ((607150.2 7090018,... ## 7 FR-211648 311 25.39949 POLYGON ((626532.2 7086540,... ## 8 FR-211647 311 45.10054 POLYGON ((669017.5 7085727,... ## 9 FR-211646 311 1005.69349 POLYGON ((618120.9 7082302,... ## 10 FR-211645 311 62.38459 POLYGON ((603423.4 7085365,... ``` --- # Use same coordinates system for map of the Pyrénées and forest layer ```r studysites <- studysites %>% * st_transform(crs = st_crs(forest)) studysites ``` ``` ## Simple feature collection with 138 features and 2 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 362039 ymin: 6139034 xmax: 690160.2 ymax: 6235505 ## Projected CRS: RGF93 / Lambert-93 ## First 10 features: ## bearpresence idsite geometry ## 1 Absence 640101 MULTIPOLYGON (((362110.5 62... ## 2 Absence 640102 MULTIPOLYGON (((372352.1 62... ## 3 Occasionnelle 640203 MULTIPOLYGON (((375125.4 62... ## 4 Occasionnelle 640204 MULTIPOLYGON (((384447.4 62... ## 5 Absence 640202 MULTIPOLYGON (((394521 6219... ## 6 Absence 640201 MULTIPOLYGON (((393646.7 62... ## 7 Absence 640301 MULTIPOLYGON (((396923.9 62... ## 8 Absence 640302 MULTIPOLYGON (((399383.2 62... ## 9 Occasionnelle 640401 MULTIPOLYGON (((399338.3 62... ## 10 Reguliere 640402 MULTIPOLYGON (((405268.8 62... ``` --- # Calculate area of each site ```r studysites %>% * mutate(area = st_area(.), .before = 1) ``` ``` ## Simple feature collection with 138 features and 3 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 362039 ymin: 6139034 xmax: 690160.2 ymax: 6235505 ## Projected CRS: RGF93 / Lambert-93 ## First 10 features: ## area bearpresence idsite geometry ## 1 79618924 [m^2] Absence 640101 MULTIPOLYGON (((362110.5 62... ## 2 122155210 [m^2] Absence 640102 MULTIPOLYGON (((372352.1 62... ## 3 64012619 [m^2] Occasionnelle 640203 MULTIPOLYGON (((375125.4 62... ## 4 36208552 [m^2] Occasionnelle 640204 MULTIPOLYGON (((384447.4 62... ## 5 67904470 [m^2] Absence 640202 MULTIPOLYGON (((394521 6219... ## 6 58455872 [m^2] Absence 640201 MULTIPOLYGON (((393646.7 62... ## 7 30310636 [m^2] Absence 640301 MULTIPOLYGON (((396923.9 62... ## 8 25866722 [m^2] Absence 640302 MULTIPOLYGON (((399383.2 62... ## 9 72321198 [m^2] Occasionnelle 640401 MULTIPOLYGON (((399338.3 62... ## 10 67797321 [m^2] Reguliere 640402 MULTIPOLYGON (((405268.8 62... ``` --- # Convert area in km<sup>2</sup> ```r studysites %>% mutate(.before = 1, area = st_area(.), * areakm2 = units::set_units(area, km^2)) ``` ``` ## Simple feature collection with 138 features and 4 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 362039 ymin: 6139034 xmax: 690160.2 ymax: 6235505 ## Projected CRS: RGF93 / Lambert-93 ## First 10 features: ## area areakm2 bearpresence idsite ## 1 79618924 [m^2] 79.61892 [km^2] Absence 640101 ## 2 122155210 [m^2] 122.15521 [km^2] Absence 640102 ## 3 64012619 [m^2] 64.01262 [km^2] Occasionnelle 640203 ## 4 36208552 [m^2] 36.20855 [km^2] Occasionnelle 640204 ## 5 67904470 [m^2] 67.90447 [km^2] Absence 640202 ## 6 58455872 [m^2] 58.45587 [km^2] Absence 640201 ## 7 30310636 [m^2] 30.31064 [km^2] Absence 640301 ## 8 25866722 [m^2] 25.86672 [km^2] Absence 640302 ## 9 72321198 [m^2] 72.32120 [km^2] Occasionnelle 640401 ## 10 67797321 [m^2] 67.79732 [km^2] Reguliere 640402 ## geometry ## 1 MULTIPOLYGON (((362110.5 62... ## 2 MULTIPOLYGON (((372352.1 62... ## 3 MULTIPOLYGON (((375125.4 62... ## 4 MULTIPOLYGON (((384447.4 62... ## 5 MULTIPOLYGON (((394521 6219... ## 6 MULTIPOLYGON (((393646.7 62... ## 7 MULTIPOLYGON (((396923.9 62... ## 8 MULTIPOLYGON (((399383.2 62... ## 9 MULTIPOLYGON (((399338.3 62... ## 10 MULTIPOLYGON (((405268.8 62... ``` --- # Define big sites (area > 300 km<sup>2</sup>) ```r studysites <- studysites %>% mutate(.before = 1, area = st_area(.), areakm2 = units::set_units(area, km^2), * bigsites = ifelse(as.numeric(areakm2) > 300, areakm2, NA)) studysites ``` ``` ## Simple feature collection with 138 features and 5 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 362039 ymin: 6139034 xmax: 690160.2 ymax: 6235505 ## Projected CRS: RGF93 / Lambert-93 ## First 10 features: ## area areakm2 bigsites bearpresence idsite ## 1 79618924 [m^2] 79.61892 [km^2] NA Absence 640101 ## 2 122155210 [m^2] 122.15521 [km^2] NA Absence 640102 ## 3 64012619 [m^2] 64.01262 [km^2] NA Occasionnelle 640203 ## 4 36208552 [m^2] 36.20855 [km^2] NA Occasionnelle 640204 ## 5 67904470 [m^2] 67.90447 [km^2] NA Absence 640202 ## 6 58455872 [m^2] 58.45587 [km^2] NA Absence 640201 ## 7 30310636 [m^2] 30.31064 [km^2] NA Absence 640301 ## 8 25866722 [m^2] 25.86672 [km^2] NA Absence 640302 ## 9 72321198 [m^2] 72.32120 [km^2] NA Occasionnelle 640401 ## 10 67797321 [m^2] 67.79732 [km^2] NA Reguliere 640402 ## geometry ## 1 MULTIPOLYGON (((362110.5 62... ## 2 MULTIPOLYGON (((372352.1 62... ## 3 MULTIPOLYGON (((375125.4 62... ## 4 MULTIPOLYGON (((384447.4 62... ## 5 MULTIPOLYGON (((394521 6219... ## 6 MULTIPOLYGON (((393646.7 62... ## 7 MULTIPOLYGON (((396923.9 62... ## 8 MULTIPOLYGON (((399383.2 62... ## 9 MULTIPOLYGON (((399338.3 62... ## 10 MULTIPOLYGON (((405268.8 62... ``` --- # Map again, with area on top of big sites ```r studysites %>% ggplot() + geom_sf() + * geom_sf_label(aes(label = round(bigsites))) + labs(title = 'Brown bear big monitoring sites in the French Pyrenees mountains', subtitle = 'Big sites have area > 300km2', caption = 'Data from: French Biodiversity Agency', x = "", y = "") ``` <img src="index_files/figure-html/unnamed-chunk-23-1.png" style="display: block; margin: auto;" /> --- # Crop forest to match study area boundaries ```r forest %>% * st_crop(st_bbox(studysites)) %>% as_tibble() ``` ``` ## # A tibble: 2,504 × 4 ## ID CODE_12 AREA_HA geometry ## <chr> <chr> <dbl> <GEOMETRY [m]> ## 1 FR-174449 311 2807. MULTIPOLYGON (((479636.9 6235505, 479719.6 6235430… ## 2 FR-174405 311 2409. POLYGON ((442252.5 6235505, 442429.7 6235438, 4425… ## 3 FR-174321 311 2843. MULTIPOLYGON (((396335.1 6235505, 396346.2 6235435… ## 4 FR-174252 311 476. POLYGON ((471577.9 6235505, 471467.3 6235289, 4713… ## 5 FR-174222 311 588. POLYGON ((490987.1 6235505, 491006.8 6235337, 4909… ## 6 FR-174216 311 354. MULTIPOLYGON (((484715.2 6235505, 484694.9 6235499… ## 7 FR-174202 311 579. MULTIPOLYGON (((500091.4 6235505, 500058.3 6235481… ## 8 FR-174187 311 357. POLYGON ((509149.7 6235505, 509148.6 6235486, 5093… ## 9 FR-174186 311 6202. MULTIPOLYGON (((379133.9 6235505, 379132.4 6235483… ## 10 FR-174174 311 292. POLYGON ((451256.3 6235505, 451260.6 6235447, 4512… ## # … with 2,494 more rows ``` --- # Then intersect the two layers ```r forest %>% st_crop(st_bbox(studysites)) %>% * st_intersection(studysites) %>% as_tibble() ``` ``` ## # A tibble: 2,177 × 9 ## ID CODE_12 AREA_HA area areakm2 bigsites bearpresence idsite ## <chr> <chr> <dbl> [m^2] [km^2] <dbl> <chr> <dbl> ## 1 FR-173974 311 243. 79618924. 79.6 NA Absence 640101 ## 2 FR-173874 311 18041. 79618924. 79.6 NA Absence 640101 ## 3 FR-173842 311 289. 79618924. 79.6 NA Absence 640101 ## 4 FR-173792 311 235. 79618924. 79.6 NA Absence 640101 ## 5 FR-173749 311 153. 79618924. 79.6 NA Absence 640101 ## 6 FR-173715 311 83.7 79618924. 79.6 NA Absence 640101 ## 7 FR-173687 311 26.9 79618924. 79.6 NA Absence 640101 ## 8 FR-173686 311 2697. 79618924. 79.6 NA Absence 640101 ## 9 FR-173667 311 49.2 79618924. 79.6 NA Absence 640101 ## 10 FR-173664 311 52.6 79618924. 79.6 NA Absence 640101 ## # … with 2,167 more rows, and 1 more variable: geometry <GEOMETRY [m]> ``` --- # Get forest area for each intersected `sfg` ```r forest %>% st_crop(st_bbox(studysites)) %>% st_intersection(studysites) %>% * mutate(area = st_area(.)) %>% as_tibble() ``` ``` ## # A tibble: 2,177 × 9 ## ID CODE_12 AREA_HA area areakm2 bigsites bearpresence idsite ## <chr> <chr> <dbl> [m^2] [km^2] <dbl> <chr> <dbl> ## 1 FR-173974 311 243. 640945. 79.6 NA Absence 640101 ## 2 FR-173874 311 18041. 5354968. 79.6 NA Absence 640101 ## 3 FR-173842 311 289. 2584902. 79.6 NA Absence 640101 ## 4 FR-173792 311 235. 567599. 79.6 NA Absence 640101 ## 5 FR-173749 311 153. 1530535. 79.6 NA Absence 640101 ## 6 FR-173715 311 83.7 837290. 79.6 NA Absence 640101 ## 7 FR-173687 311 26.9 268517. 79.6 NA Absence 640101 ## 8 FR-173686 311 2697. 21933124. 79.6 NA Absence 640101 ## 9 FR-173667 311 49.2 492325. 79.6 NA Absence 640101 ## 10 FR-173664 311 52.6 526027. 79.6 NA Absence 640101 ## # … with 2,167 more rows, and 1 more variable: geometry <GEOMETRY [m]> ``` --- # Sum forest over all study sites ```r forestpyrenees <- forest %>% st_crop(st_bbox(studysites)) %>% st_intersection(studysites) %>% mutate(area = st_area(.)) %>% * group_by(idsite) %>% # groups a data frame by variables * summarise(areaforest = sum(area)) %>% # perform group-wise summaries as_tibble() %>% select(-geometry) forestpyrenees ``` ``` ## # A tibble: 138 × 2 ## idsite areaforest ## <dbl> [m^2] ## 1 11403 157762282. ## 2 11404 32526319. ## 3 90101 137218858. ## 4 90102 86944708. ## 5 90103 82872880. ## 6 90104 42241267. ## 7 90201 19224150. ## 8 90202 26309123. ## 9 90203 31962484. ## 10 90301 8958589. ## # … with 128 more rows ``` --- # Join `sf` and `tibble` objects ## More info [here](https://r-spatial.github.io/sf/reference/tidyverse.html) ```r studysites %>% * inner_join(forestpyrenees, by = 'idsite') ``` ``` ## Simple feature collection with 138 features and 6 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 362039 ymin: 6139034 xmax: 690160.2 ymax: 6235505 ## Projected CRS: RGF93 / Lambert-93 ## First 10 features: ## area areakm2 bigsites bearpresence idsite ## 1 79618924 [m^2] 79.61892 [km^2] NA Absence 640101 ## 2 122155210 [m^2] 122.15521 [km^2] NA Absence 640102 ## 3 64012619 [m^2] 64.01262 [km^2] NA Occasionnelle 640203 ## 4 36208552 [m^2] 36.20855 [km^2] NA Occasionnelle 640204 ## 5 67904470 [m^2] 67.90447 [km^2] NA Absence 640202 ## 6 58455872 [m^2] 58.45587 [km^2] NA Absence 640201 ## 7 30310636 [m^2] 30.31064 [km^2] NA Absence 640301 ## 8 25866722 [m^2] 25.86672 [km^2] NA Absence 640302 ## 9 72321198 [m^2] 72.32120 [km^2] NA Occasionnelle 640401 ## 10 67797321 [m^2] 67.79732 [km^2] NA Reguliere 640402 ## areaforest geometry ## 1 42105194 [m^2] MULTIPOLYGON (((362110.5 62... ## 2 60136750 [m^2] MULTIPOLYGON (((372352.1 62... ## 3 26652318 [m^2] MULTIPOLYGON (((375125.4 62... ## 4 21003797 [m^2] MULTIPOLYGON (((384447.4 62... ## 5 42325894 [m^2] MULTIPOLYGON (((394521 6219... ## 6 22478565 [m^2] MULTIPOLYGON (((393646.7 62... ## 7 13955020 [m^2] MULTIPOLYGON (((396923.9 62... ## 8 19648844 [m^2] MULTIPOLYGON (((399383.2 62... ## 9 31347244 [m^2] MULTIPOLYGON (((399338.3 62... ## 10 21175525 [m^2] MULTIPOLYGON (((405268.8 62... ``` --- # Calculate forest cover ```r covariates <- studysites %>% inner_join(forestpyrenees, by = 'idsite') %>% mutate(.before = 1, * forestcover = areaforest / area) ``` --- # Map forest cover ```r covariates %>% ggplot() + aes(fill = as.numeric(forestcover)) + geom_sf(lwd = 0.1) + scale_fill_viridis_c( labels = scales::percent_format(), #<< format percentage name = 'Cover', alpha = 0.7) + #<< control transparency labs(title = 'Map of forest cover in the Pyrenees mountains', subtitle = 'Source: Corine Land Cover 2012') ``` --- # Map forest cover <img src="index_files/figure-html/unnamed-chunk-31-1.png" style="display: block; margin: auto;" /> --- # Interactive map with `mapview` ## More info [here](https://r-spatial.github.io/mapview/) ```r library(mapview) *mapview(covariates, zcol = "bearpresence") ```
--- # Interactive map with `mapview` ```r covariates <- covariates %>% mutate(forestcover = as.numeric(forestcover)) *mapview(covariates, zcol = "forestcover", map.types = "OpenTopoMap") ```
```r # map.types = "CartoDB.Positron", "CartoDB.DarkMatter", "OpenStreetMap", # "Esri.WorldImagery" ``` --- # Human density * Human density might be a driver of brown bear distribution * We use data on population size of cities in France * Population data is available from **IGN Admin-Express** database [here](https://geoservices.ign.fr/documentation/diffusion/telechargement-donnees-libres.html#admin-express) --- # Human density ```r load("shp/france_population_data_2016.RData") glimpse(df.france) ``` ``` ## Rows: 35,798 ## Columns: 18 ## $ ID_GEOFLA <fct> COMMUNE00000000000000001, COMMUNE00000000000000002, COMMUNE… ## $ CODE_COM <fct> 216, 033, 009, 225, 890, 018, 113, 319, 097, 070, 046, 524,… ## $ INSEE_COM <fct> 32216, 47033, 32009, 38225, 62890, 08018, 32113, 10319, 060… ## $ NOM_COM <fct> LOURTIES-MONBRUN, BOUDY-DE-BEAUREGARD, ARMOUS-ET-CAU, AUTRA… ## $ STATUT <fct> Commune simple, Commune simple, Commune simple, Commune sim… ## $ X_CHF_LIEU <int> 500820, 516424, 472979, 898640, 640049, 824246, 461332, 746… ## $ Y_CHF_LIEU <int> 6264958, 6384852, 6278963, 6450689, 7028672, 6908952, 63007… ## $ X_CENTROID <int> 500515, 515575, 473004, 898625, 640115, 824391, 460721, 747… ## $ Y_CENTROID <int> 6265413, 6385938, 6278937, 6451597, 7029900, 6908954, 63022… ## $ Z_MOYEN <int> 252, 112, 221, 1234, 79, 125, 134, 167, 752, 438, 1276, 362… ## $ SUPERFICIE <dbl> 966, 1019, 932, 3371, 1023, 438, 919, 1904, 2217, 2667, 306… ## $ POPULATION <dbl> 139, 414, 95, 2973, 178, 80, 97, 362, 296, 901, 10, 166, 94… ## $ CODE_ARR <fct> 3, 3, 3, 1, 4, 4, 2, 3, 2, 2, 2, 3, 2, 1, 1, 3, 2, 1, 3, 3,… ## $ CODE_DEPT <fct> 32, 47, 32, 38, 62, 08, 32, 10, 06, 42, 31, 71, 53, 16, 38,… ## $ NOM_DEPT <fct> GERS, LOT-ET-GARONNE, GERS, ISERE, PAS-DE-CALAIS, ARDENNES,… ## $ CODE_REG <fct> 76, 75, 76, 84, 32, 44, 76, 44, 93, 84, 76, 27, 52, 75, 84,… ## $ NOM_REG <fct> LANGUEDOC-ROUSSILLON-MIDI-PYRENEES, AQUITAINE-LIMOUSIN-POIT… ## $ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((499484.6 62..., MULTIPOLYGON (… ``` --- # Transform into lower case ```r colnames(df.france) <- colnames(df.france) %>% * str_to_lower() colnames(df.france) ``` ``` ## [1] "id_geofla" "code_com" "insee_com" "nom_com" "statut" ## [6] "x_chf_lieu" "y_chf_lieu" "x_centroid" "y_centroid" "z_moyen" ## [11] "superficie" "population" "code_arr" "code_dept" "nom_dept" ## [16] "code_reg" "nom_reg" "geometry" ``` --- # Have a look to the distribution ```r df.france$population %>% * summary() ``` ``` ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0 197 442 1779 1100 458298 ``` --- # Calculate density ```r df.france <- df.france %>% * mutate(density = population/superficie*100) as_tibble(df.france) ``` ``` ## # A tibble: 35,798 × 19 ## id_geofla code_com insee_com nom_com statut x_chf_lieu y_chf_lieu x_centroid ## <fct> <fct> <fct> <fct> <fct> <int> <int> <int> ## 1 COMMUNE00… 216 32216 LOURTI… Commu… 500820 6264958 500515 ## 2 COMMUNE00… 033 47033 BOUDY-… Commu… 516424 6384852 515575 ## 3 COMMUNE00… 009 32009 ARMOUS… Commu… 472979 6278963 473004 ## 4 COMMUNE00… 225 38225 AUTRAN… Commu… 898640 6450689 898625 ## 5 COMMUNE00… 890 62890 WILLEM… Commu… 640049 7028672 640115 ## 6 COMMUNE00… 018 08018 ARDEUI… Commu… 824246 6908952 824391 ## 7 COMMUNE00… 113 32113 CRAVEN… Commu… 461332 6300782 460721 ## 8 COMMUNE00… 319 10319 RIGNY-… Commu… 746925 6790005 747181 ## 9 COMMUNE00… 097 06097 PIERRE… Commu… 1028827 6315717 1027327 ## 10 COMMUNE00… 070 42070 CORDEL… Commu… 782215 6538794 782159 ## # … with 35,788 more rows, and 11 more variables: y_centroid <int>, ## # z_moyen <int>, superficie <dbl>, population <dbl>, code_arr <fct>, ## # code_dept <fct>, nom_dept <fct>, code_reg <fct>, nom_reg <fct>, ## # geometry <MULTIPOLYGON [m]>, density <dbl> ``` --- # Show density ```r df.france %>% pull(density) %>% head() ``` ``` ## [1] 14.38923 40.62807 10.19313 88.19341 17.39980 18.26484 ``` --- # Sum population size over sites ```r df.pyrenees <- df.france %>% st_transform(crs = st_crs(forest)) %>% st_crop(st_bbox(covariates)) %>% st_intersection(covariates) %>% st_set_geometry(NULL) %>% * group_by(idsite) %>% * summarise(humpop = sum(population)) df.pyrenees ``` ``` ## # A tibble: 138 × 2 ## idsite humpop ## <dbl> <dbl> ## 1 11403 10711 ## 2 11404 7576 ## 3 90101 33154 ## 4 90102 17391 ## 5 90103 19126 ## 6 90104 28000 ## 7 90201 8203 ## 8 90202 2554 ## 9 90203 3505 ## 10 90301 285 ## # … with 128 more rows ``` --- # Join, then calculate density ```r covariates <- covariates %>% * inner_join(df.pyrenees, by = 'idsite') %>% mutate(.before = 1, * humdens = humpop / (area/1000000)) as_tibble(covariates) ``` ``` ## # A tibble: 138 × 10 ## humdens forestcover area areakm2 bigsites bearpresence idsite areaforest ## [1/m^2] <dbl> [m^2] [km^2] <dbl> <chr> <dbl> [m^2] ## 1 15.7 0.529 7.96e7 79.6 NA Absence 640101 42105194. ## 2 16.4 0.492 1.22e8 122. NA Absence 640102 60136750. ## 3 10.0 0.416 6.40e7 64.0 NA Occasionnel… 640203 26652318. ## 4 51.8 0.580 3.62e7 36.2 NA Occasionnel… 640204 21003797. ## 5 30.4 0.623 6.79e7 67.9 NA Absence 640202 42325894. ## 6 68.1 0.385 5.85e7 58.5 NA Absence 640201 22478565. ## 7 58.8 0.460 3.03e7 30.3 NA Absence 640301 13955020. ## 8 70.6 0.760 2.59e7 25.9 NA Absence 640302 19648844. ## 9 39.8 0.433 7.23e7 72.3 NA Occasionnel… 640401 31347244. ## 10 12.7 0.312 6.78e7 67.8 NA Reguliere 640402 21175525. ## # … with 128 more rows, and 2 more variables: humpop <dbl>, ## # geometry <MULTIPOLYGON [m]> ``` --- # Map human population density ```r covariates %>% ggplot() + aes(fill = as.numeric(humdens)) + geom_sf(lwd = 0.1) + scale_fill_viridis_c( name = bquote('Density\n(people per km'^2*')'), alpha = 0.7) + labs(title = 'Human population density', subtitle = 'Source: IGN GEOFLA 2016') ``` --- # Map human population density <img src="index_files/figure-html/unnamed-chunk-42-1.png" style="display: block; margin: auto;" /> --- class: inverse, middle, center # Spatial operations: distance --- # Distance to highways * Distance to highways might be a driver of brown bear distribution * We use data from [Route500 database](https://geoservices.ign.fr/documentation/diffusion/telechargement-donnees-libres.html#route-500) --- # Read in data ```r roads <- st_read("shp/TRONCON_ROUTE.shp") ``` ``` ## Reading layer `TRONCON_ROUTE' from data source ## `/Users/oliviergimenez/Dropbox/OG/GITHUB/intro_spatialR/shp/TRONCON_ROUTE.shp' ## using driver `ESRI Shapefile' ## Simple feature collection with 339250 features and 12 fields ## Geometry type: LINESTRING ## Dimension: XY ## Bounding box: xmin: 100076.5 ymin: 6021027 xmax: 1329738 ymax: 7120638 ## Projected CRS: RGF93 / Lambert-93 ``` --- # Read in data ```r as_tibble(roads) ``` ``` ## # A tibble: 339,250 × 13 ## ID_RTE500 VOCATION NB_CHAUSSE NB_VOIES ETAT ACCES RES_VERT SENS RES_EUROPE ## <int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> ## 1 1 "Liaison… "1 chauss… "1 voie… "Rev… Libre N'appar… Doub… <NA> ## 2 2 "Liaison… "1 chauss… "1 voie… "Rev… Libre N'appar… Doub… <NA> ## 3 3 "Liaison… "1 chauss… "1 voie… "Rev… Libre N'appar… Doub… <NA> ## 4 4 "Liaison… "1 chauss… "1 voie… "Rev… Libre N'appar… Doub… <NA> ## 5 5 "Liaison… "1 chauss… "1 voie… "Rev… Libre N'appar… Doub… <NA> ## 6 6 "Liaison… "1 chauss… "1 voie… "Rev… Libre N'appar… Doub… <NA> ## 7 7 "Liaison… "1 chauss… "2 voie… "Rev… Libre N'appar… Sens… <NA> ## 8 8 "Liaison… "1 chauss… "1 voie… "Rev… Libre N'appar… Doub… <NA> ## 9 9 "Liaison… "1 chauss… "1 voie… "Rev… Libre N'appar… Doub… <NA> ## 10 10 "Liaison… "1 chauss… "1 voie… "Rev… Libre N'appar… Doub… <NA> ## # … with 339,240 more rows, and 4 more variables: NUM_ROUTE <chr>, ## # CLASS_ADM <chr>, LONGUEUR <dbl>, geometry <LINESTRING [m]> ``` --- # Focus on highways ```r highways <- roads %>% * filter(CLASS_ADM == "Autoroute") as_tibble(highways) ``` ``` ## # A tibble: 4,506 × 13 ## ID_RTE500 VOCATION NB_CHAUSSE NB_VOIES ETAT ACCES RES_VERT SENS RES_EUROPE ## <int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> ## 1 32 Type aut… "2 chauss… "Sans o… "Rev… "A p… Apparti… Doub… <NA> ## 2 33 Type aut… "1 chauss… "2 voie… "Rev… "A p… Apparti… Doub… <NA> ## 3 1041 Type aut… "1 chauss… "1 voie… "Rev… "Lib… Apparti… Sens… <NA> ## 4 2460 Type aut… "2 chauss… "Sans o… "Rev… "Lib… N'appar… Doub… E60 ## 5 3074 Type aut… "1 chauss… "1 voie… "Rev… "Lib… N'appar… Sens… <NA> ## 6 3198 Type aut… "1 chauss… "1 voie… "Rev… "A p… Apparti… Sens… E712 ## 7 3200 Type aut… "2 chauss… "Sans o… "Rev… "Lib… Apparti… Doub… <NA> ## 8 3201 Type aut… "2 chauss… "Sans o… "Rev… "Lib… Apparti… Doub… <NA> ## 9 3202 Type aut… "2 chauss… "Sans o… "Rev… "Lib… Apparti… Doub… <NA> ## 10 3363 Type aut… "2 chauss… "Sans o… "Rev… "A p… Apparti… Doub… E713 ## # … with 4,496 more rows, and 4 more variables: NUM_ROUTE <chr>, ## # CLASS_ADM <chr>, LONGUEUR <dbl>, geometry <LINESTRING [m]> ``` --- # Reproject and crop to match France extent ```r highways <- highways %>% * st_transform(crs = st_crs(forest)) %>% * st_crop(st_bbox(df.france)) as_tibble(highways) ``` ``` ## # A tibble: 4,494 × 13 ## ID_RTE500 VOCATION NB_CHAUSSE NB_VOIES ETAT ACCES RES_VERT SENS RES_EUROPE ## <int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> ## 1 32 Type aut… "2 chauss… "Sans o… "Rev… "A p… Apparti… Doub… <NA> ## 2 33 Type aut… "1 chauss… "2 voie… "Rev… "A p… Apparti… Doub… <NA> ## 3 1041 Type aut… "1 chauss… "1 voie… "Rev… "Lib… Apparti… Sens… <NA> ## 4 2460 Type aut… "2 chauss… "Sans o… "Rev… "Lib… N'appar… Doub… E60 ## 5 3074 Type aut… "1 chauss… "1 voie… "Rev… "Lib… N'appar… Sens… <NA> ## 6 3198 Type aut… "1 chauss… "1 voie… "Rev… "A p… Apparti… Sens… E712 ## 7 3200 Type aut… "2 chauss… "Sans o… "Rev… "Lib… Apparti… Doub… <NA> ## 8 3201 Type aut… "2 chauss… "Sans o… "Rev… "Lib… Apparti… Doub… <NA> ## 9 3202 Type aut… "2 chauss… "Sans o… "Rev… "Lib… Apparti… Doub… <NA> ## 10 3363 Type aut… "2 chauss… "Sans o… "Rev… "A p… Apparti… Doub… E713 ## # … with 4,484 more rows, and 4 more variables: NUM_ROUTE <chr>, ## # CLASS_ADM <chr>, LONGUEUR <dbl>, geometry <LINESTRING [m]> ``` --- # Get centroids of each monitoring sites ```r centroids <- covariates %>% * st_centroid() as_tibble(centroids) ``` ``` ## # A tibble: 138 × 10 ## humdens forestcover area areakm2 bigsites bearpresence idsite areaforest ## [1/m^2] <dbl> [m^2] [km^2] <dbl> <chr> <dbl> [m^2] ## 1 15.7 0.529 7.96e7 79.6 NA Absence 640101 42105194. ## 2 16.4 0.492 1.22e8 122. NA Absence 640102 60136750. ## 3 10.0 0.416 6.40e7 64.0 NA Occasionnel… 640203 26652318. ## 4 51.8 0.580 3.62e7 36.2 NA Occasionnel… 640204 21003797. ## 5 30.4 0.623 6.79e7 67.9 NA Absence 640202 42325894. ## 6 68.1 0.385 5.85e7 58.5 NA Absence 640201 22478565. ## 7 58.8 0.460 3.03e7 30.3 NA Absence 640301 13955020. ## 8 70.6 0.760 2.59e7 25.9 NA Absence 640302 19648844. ## 9 39.8 0.433 7.23e7 72.3 NA Occasionnel… 640401 31347244. ## 10 12.7 0.312 6.78e7 67.8 NA Reguliere 640402 21175525. ## # … with 128 more rows, and 2 more variables: humpop <dbl>, ## # geometry <POINT [m]> ``` --- # Then distance from centroids to highways ```r dtohighways <- highways %>% * st_distance(centroids, by_element = F) head(dtohighways) ``` ``` ## Units: [m] ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] ## [1,] 490490.3 488265.3 493003.4 492821.9 487237.4 481060.4 481121.5 486089.2 ## [2,] 515331.5 507719.9 508080.6 503082.1 499732.8 494803.3 490230.3 494329.4 ## [3,] 718550.3 716685.9 721602.9 721642.1 716013.6 709840.7 710080.5 715040.7 ## [4,] 898734.7 890836.7 890814.1 885459.0 882398.2 877681.9 872780.4 876701.3 ## [5,] 919059.7 915277.9 918978.3 917377.9 912208.2 906119.3 904794.9 909730.3 ## [6,] 679572.2 671341.9 670914.4 665250.8 662457.8 657966.6 652762.1 656514.7 ## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] ## [1,] 492001.5 498394.1 486791.1 480913.1 477857.5 486211.6 494968.5 502324.0 ## [2,] 496296.6 499690.1 490796.0 486687.6 480646.9 485431.1 491383.7 498446.8 ## [3,] 721040.1 727480.1 715870.2 709978.7 707015.7 715405.3 724179.1 731511.8 ## [4,] 878203.4 881160.3 872809.8 868989.3 862787.5 867045.2 872482.5 879296.8 ## [5,] 914698.5 920441.5 909152.7 903547.5 899397.6 907014.2 915278.5 922884.7 ## [6,] 657646.3 660295.5 652347.9 648767.0 642453.8 646315.8 651414.2 658068.7 ## [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] ## [1,] 503492.4 496943.5 485118.6 475420.8 476127.3 483837.8 489587.5 497561.7 ## [2,] 496659.6 489140.4 480224.3 474727.6 468192.9 475111.9 478057.0 484235.5 ## [3,] 732726.8 726207.8 714378.6 704653.7 705421.4 713128.9 718882.7 726857.0 ## [4,] 877162.5 869748.3 861490.3 856646.4 849448.8 856046.2 858513.1 864218.4 ## [5,] 923129.9 915942.4 904508.9 895707.7 894016.9 901863.1 906940.6 914692.3 ## [6,] 655757.1 648417.1 640554.2 636163.0 628543.0 634931.2 637135.9 642599.0 ## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] ## [1,] 502169.7 481567.1 494352.4 474992.7 481480.9 491781.3 499425.2 499168.0 ## [2,] 490669.7 469656.0 476853.7 461019.4 463382.6 471086.6 476109.0 472448.7 ## [3,] 731455.3 710861.2 723633.3 704268.0 710728.3 721024.2 728653.8 728335.1 ## [4,] 870699.8 850345.0 856472.2 841725.9 843418.0 850425.0 854812.8 850718.4 ## [5,] 920160.9 898366.7 909839.9 890702.8 896063.3 905972.7 913051.2 911505.9 ## [6,] 649084.4 629112.4 634714.8 620529.0 621881.8 628569.4 632707.6 628477.0 ## [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] ## [1,] 473359.6 479866.0 497112.3 507856.6 507231.0 494691.9 488694.2 483485.9 ## [2,] 454703.7 453431.2 467322.2 476102.5 470137.5 459471.4 458986.0 448385.6 ## [3,] 702571.7 708940.2 726197.0 736947.8 736155.7 723591.0 717731.1 712297.4 ## [4,] 834985.8 832571.6 845269.9 853278.5 846545.5 836778.2 837326.5 826258.2 ## [5,] 887280.1 891245.8 908168.2 918722.5 915966.5 903490.5 899318.1 891690.5 ## [6,] 613592.4 610688.8 622940.8 630709.8 623801.3 614284.5 615143.6 603956.2 ## [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] ## [1,] 472425.6 474027.4 488915.6 498881.3 506604.3 502202.7 506817.8 491798.5 ## [2,] 444035.2 433861.1 445765.1 455722.7 459705.6 450874.4 450573.0 439475.0 ## [3,] 701390.0 702510.7 717416.6 727485.9 735103.1 730408.3 734823.0 719836.0 ## [4,] 823267.2 811533.3 822217.3 831629.6 834577.6 825257.8 823794.0 814319.2 ## [5,] 882612.4 879602.2 894162.6 904699.2 911295.2 904767.3 907540.2 893345.5 ## [6,] 601448.0 589212.0 599540.7 608806.6 611544.2 602158.7 600523.1 591323.7 ## [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] ## [1,] 475689.4 486361.4 496396.0 505461.2 508708.4 502688.9 500488.5 490572.9 ## [2,] 424506.8 429328.2 438305.6 446520.2 444575.9 440943.0 435418.2 425741.0 ## [3,] 703638.9 714160.8 724241.6 733348.4 736385.9 730425.2 728033.4 718018.8 ## [4,] 800498.9 803698.2 811844.1 819315.0 816130.5 813379.1 807353.6 798415.8 ## [5,] 876778.2 885620.4 895773.9 904969.2 906144.8 900841.1 897066.2 886693.2 ## [6,] 577771.0 580642.5 588639.9 595996.4 592678.6 590026.9 583954.6 575119.6 ## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] ## [1,] 491709.5 492957.5 487799.3 485401.4 483950.1 481832.5 476828.7 471875.6 ## [2,] 430317.6 422020.5 416952.2 420583.3 415038.0 416914.8 414854.3 409226.4 ## [3,] 719353.2 720058.9 714835.3 712785.3 711060.6 709165.2 704269.5 699212.6 ## [4,] 803557.0 793332.5 788670.2 793623.5 787416.5 790187.9 788996.4 783592.1 ## [5,] 889402.5 886479.0 881057.7 881223.8 877854.7 877396.4 873383.9 867825.1 ## [6,] 580326.3 569891.7 565279.4 570385.1 564104.3 566989.5 565935.9 560582.9 ## [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] ## [1,] 471133.2 478433.0 467224.9 478267.6 483336.1 488707.0 495988.4 496239.6 ## [2,] 399364.7 406573.1 385972.3 398082.5 404211.6 411855.8 418506.7 416067.2 ## [3,] 697846.1 705258.1 693080.8 704433.1 709681.5 715323.2 722670.3 722716.4 ## [4,] 772177.6 778808.1 757321.8 768701.7 774622.2 782294.5 788225.0 785182.4 ## [5,] 862932.5 870655.9 854311.5 866587.7 872478.1 879247.4 886663.8 885660.5 ## [6,] 548953.8 555495.4 533950.7 545236.0 551128.2 558791.0 564657.4 561578.9 ## [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80] ## [1,] 497351.8 498488.9 494125.2 490926.0 488453.9 484933.1 491710.4 499237.4 ## [2,] 415251.3 413444.7 408765.0 410831.2 403627.3 397313.1 401986.7 409372.4 ## [3,] 723690.0 724597.3 720129.6 717321.9 714402.8 710555.6 717264.9 724922.2 ## [4,] 783847.1 781276.5 776908.1 780406.1 772389.3 765776.1 769352.1 776009.4 ## [5,] 885918.3 885703.2 880930.2 880069.9 875177.3 870075.6 876229.6 884132.2 ## [6,] 560215.4 557610.7 553260.5 556838.6 548779.6 542157.1 545677.5 552306.1 ## [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] ## [1,] 505515.8 505438.5 470916.3 483643.9 489460.8 496111.3 497431.5 507472.6 ## [2,] 412933.8 410735.6 379468.7 391089.7 392533.5 398195.4 394546.3 405002.5 ## [3,] 731056.5 730767.8 695849.1 708748.2 714227.7 720915.4 721706.2 732013.8 ## [4,] 778323.6 775599.8 748415.3 758568.2 758416.7 763168.7 758117.5 767632.3 ## [5,] 889416.7 888267.0 853254.0 866253.9 870215.7 876761.0 875580.2 886418.1 ## [6,] 554595.6 551867.9 524849.4 534908.1 534704.6 539440.1 534385.1 543912.1 ## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] ## [1,] 506054.1 498338.2 514893.4 514576.8 506721.1 504218.2 499282.7 494772.2 ## [2,] 408174.2 402349.1 408192.4 404061.2 400625.9 395190.2 389825.2 383362.2 ## [3,] 731067.6 723394.5 739103.6 738316.2 730825.8 727907.2 722800.1 717925.1 ## [4,] 772163.4 767580.9 768863.5 763662.3 762345.9 756370.3 751438.7 744954.9 ## [5,] 887298.6 880104.5 892013.9 889628.9 883711.4 879496.8 874044.2 868209.0 ## [6,] 548430.8 543853.9 545195.0 540029.9 538641.2 532675.0 527729.5 521242.2 ## [,97] [,98] [,99] [,100] [,101] [,102] [,103] [,104] ## [1,] 488371.7 483855.2 477691.3 473978.4 496337.7 502441.2 514891.1 506367.4 ## [2,] 380227.6 383135.3 375615.0 366123.3 378882.7 385990.2 398476.6 384521.1 ## [3,] 711784.3 708064.8 701585.2 697041.1 718709.0 725112.1 737873.3 728376.7 ## [4,] 743419.4 748678.0 741485.8 731017.7 738598.1 745257.4 756274.6 741689.9 ## [5,] 863182.6 862304.7 855037.0 848019.0 866538.8 873557.1 886721.0 874673.9 ## [6,] 519686.9 524963.7 517783.8 507299.9 514931.7 521615.8 532720.0 518144.1 ## [,105] [,106] [,107] [,108] [,109] [,110] [,111] [,112] ## [1,] 519938.4 527220.8 522444.3 524730.6 514998.4 520210.1 527786.4 501791.0 ## [2,] 396238.3 397922.8 392269.6 392259.2 385965.6 387611.9 393333.5 376421.7 ## [3,] 742023.7 748650.8 743618.1 745601.8 736155.1 740942.5 748418.5 723151.7 ## [4,] 751115.8 749975.1 744645.4 743526.2 739716.4 739446.5 743476.2 733090.0 ## [5,] 887920.4 892321.2 886786.9 887826.3 879625.6 883000.9 889846.9 867821.5 ## [6,] 527740.9 526833.3 521456.6 520435.4 516397.1 516288.6 520495.6 509555.7 ## [,113] [,114] [,115] [,116] [,117] [,118] [,119] [,120] ## [1,] 508003.7 503868.6 512530.8 514341.4 522197.3 530088.4 522084.0 520447.0 ## [2,] 375353.1 366072.6 370142.3 377286.7 383421.2 390751.5 377493.1 371014.0 ## [3,] 728391.5 723270.9 731387.1 734178.2 741961.6 749970.3 740824.1 738243.0 ## [4,] 728845.6 718422.2 719594.6 728378.6 732770.5 738765.7 724686.5 716589.1 ## [5,] 870117.3 862663.5 868972.2 874148.3 881356.7 889290.0 877634.2 872848.8 ## [6,] 505514.4 495148.4 496584.7 505248.5 509827.1 515991.7 501936.7 493981.3 ## [,121] [,122] [,123] [,124] [,125] [,126] [,127] [,128] ## [1,] 507456.2 509163.1 522078.4 530625.7 535061.9 537928.8 547655.2 557358.3 ## [2,] 358845.6 353092.4 365308.8 373897.1 383797.5 392367.0 392536.2 395338.5 ## [3,] 725007.3 725300.9 738486.1 747299.3 752931.8 756931.1 765125.9 773698.3 ## [4,] 706803.2 697915.2 707659.5 714839.0 726331.8 736807.8 731423.8 729457.2 ## [5,] 859728.1 856708.0 869786.9 878749.8 886964.6 893612.4 897572.2 903008.2 ## [6,] 483884.3 475287.0 485393.2 492809.0 504163.3 514452.4 509808.9 508586.0 ## [,129] [,130] [,131] [,132] [,133] [,134] [,135] [,136] ## [1,] 545520.1 549730.6 552381.4 540689.9 538668.6 535793.2 539028.1 488084.1 ## [2,] 382045.9 401708.2 411238.4 401014.0 404229.0 404722.4 411295.1 464599.6 ## [3,] 761232.2 768594.2 772517.5 760787.7 759575.9 757160.1 760976.0 717268.7 ## [4,] 717458.2 743375.1 755404.9 747401.3 752903.2 755022.9 762386.6 843750.4 ## [5,] 889718.5 904360.1 911474.6 900177.5 901258.3 900263.7 905621.5 901042.0 ## [6,] 496270.1 521435.2 533189.9 524858.7 530118.4 532071.1 539380.9 621839.7 ## [,137] [,138] ## [1,] 511083.6 498086.9 ## [2,] 341971.2 326549.7 ## [3,] 724481.3 710457.4 ## [4,] 681047.7 667126.5 ## [5,] 849928.1 834630.9 ## [6,] 459013.4 444788.6 ``` --- # Convert distance to highways into numeric values and keep only minimal distance to highways ```r units(dtohighways) <- units::as_units("km") dtohighwaysnum <- matrix(as.numeric(dtohighways), nrow = nrow(dtohighways), ncol = ncol(dtohighways)) dtohighwaysnum <- apply(dtohighwaysnum,2,min) head(dtohighwaysnum) ``` ``` ## [1] 48.67387 49.63361 53.25704 50.65845 45.78141 39.89380 ``` --- # Map the distance to highways ```r covariates <- covariates %>% add_column(dtohighwaysnum) covariates %>% ggplot() + geom_sf(lwd = 0.1, aes(fill = dtohighwaysnum)) + scale_fill_viridis_c(name = 'distance\n(km)',alpha = 0.7) + geom_sf(data = highways, aes(color = 'red'), show.legend = "line") + scale_color_manual(values = "red", labels = "", name = "highway") + coord_sf(xlim = st_bbox(covariates)[c(1,3)], ylim = st_bbox(covariates)[c(2,4)]) + # what if you turn this off? labs(title = 'Distance to highways in the Pyrénées', subtitle = 'Source: Route500') ``` --- # Map the distance to highways <img src="index_files/figure-html/unnamed-chunk-51-1.png" style="display: block; margin: auto;" /> --- class: inverse, middle, center # Moooooore --- # Generate geometries from (X,Y) columns * We use [EPSG code = 4326](https://epsg.io/4326) as CRS for WGS84 coordinates * Let's make a tibble with 10 points ```r random_pts <- tibble( id_point = seq_len(10), lat = c(42.96, 42.71, 42.72, 42.95, 42.96, 42.72, 42.68, 42.82, 42.79, 42.85), lon = c(0.31, 0.66, 0.75, 1.60, 0.58, 1.87, 1.07, 1.05, 0.36, 0.06) ) *sf_pts <- st_as_sf( * random_pts, coords = c("lon", "lat"), crs=4326, remove=FALSE) ``` --- # Generate geometries from (X,Y) columns ```r sf_pts ``` ``` ## Simple feature collection with 10 features and 3 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: 0.06 ymin: 42.68 xmax: 1.87 ymax: 42.96 ## Geodetic CRS: WGS 84 ## # A tibble: 10 × 4 ## id_point lat lon geometry ## * <int> <dbl> <dbl> <POINT [°]> ## 1 1 43.0 0.31 (0.31 42.96) ## 2 2 42.7 0.66 (0.66 42.71) ## 3 3 42.7 0.75 (0.75 42.72) ## 4 4 43.0 1.6 (1.6 42.95) ## 5 5 43.0 0.58 (0.58 42.96) ## 6 6 42.7 1.87 (1.87 42.72) ## 7 7 42.7 1.07 (1.07 42.68) ## 8 8 42.8 1.05 (1.05 42.82) ## 9 9 42.8 0.36 (0.36 42.79) ## 10 10 42.8 0.06 (0.06 42.85) ``` --- # Display our points ```r ggplot() + geom_sf(data = studysites) + geom_sf(data = sf_pts) + labs(title = 'Study area and some random points') ``` <img src="index_files/figure-html/unnamed-chunk-54-1.png" style="display: block; margin: auto;" /> --- # Find points that intersect polygons * `st_intersects` tells us which polygons intersect which points * All objects should have same CRS * If there is no intersection, we get an empty result ```r sf_pts <- sf_pts %>% st_transform(st_crs(studysites)) *result <- st_intersects(sf_pts, studysites) result ``` ``` ## Sparse geometry binary predicate list of length 10, where the predicate ## was `intersects' ## 1: 42 ## 2: 53 ## 3: (empty) ## 4: 97 ## 5: 50 ## 6: 109 ## 7: (empty) ## 8: 74 ## 9: 46 ## 10: 38 ``` --- # Find points that intersect polygons * `st_intersects` tells us which polygons intersect which points * All objects should have same CRS * If there is no intersection, we get an empty result ```r sf_pts <- sf_pts %>% st_transform(st_crs(studysites)) result <- st_intersects(sf_pts, studysites) *subset_pts <- sf_pts %>% * add_column(nb_intersected = sapply(result, length)) %>% * filter(nb_intersected > 0) ``` --- # Highlight points intersecting polygons ```r ggplot() + geom_sf(data = studysites) + geom_sf(data = sf_pts) + geom_sf(data = subset_pts, color = "red", size = 1.3) + labs(title = 'Study area and some random points') ``` <img src="index_files/figure-html/unnamed-chunk-57-1.png" style="display: block; margin: auto;" /> --- # How to make a grid * Let's create an hexagonal grid with 10 km cells covering the study area ```r *hex_grid <- st_make_grid(studysites, cellsize = 10000, square = FALSE) ggplot() + geom_sf(data = studysites, color = NA, fill = "blue", alpha = 0.4) + geom_sf(data = hex_grid, fill = NA, lwd = 0.2) + labs(title = 'A hexagonal grid') ``` <img src="index_files/figure-html/unnamed-chunk-58-1.png" style="display: block; margin: auto;" /> --- class: inverse, middle, center # Wrap up --- ## Geometric calculations **Geometric operations** on vector layers can conceptually be divided into **three groups** according to their output: * **Numeric** values: Functions that summarize geometrical properties of: * A **single layer** (e.g. area, length) * A **pair of layers** (e.g. distance) * **Logical** values: Functions that evaluate whether a certain condition holds true, regarding: * A **single layer** (e.g. geometry is valid) * A **pair of layers** (e.g. feature A intersects feature B) * **Spatial** layers: Functions that create a new layer based on: * A **single layer** (e.g. centroids) * A **pair of layers** (e.g. intersection area) --- ## Numeric * Several functions to calculate **numeric geometric properties** of vector layers: * `st_length` * `st_area` * `st_distance` * `st_bbox` * ... --- ## Logical * Given two layers, `x` and `y`, the following **logical geometric functions** check whether each feature in `x` maintains the specified **relation** with each feature in `y`: * `st_intersects` * `st_disjoint` * `st_touches` * `st_crosses` * `st_within` * `st_contains` * `st_overlaps` * `st_covers` * `st_equals` * ... --- ## Spatial * Common **geometry-generating** functions applicable to **individual** geometries: * `st_centroid` * `st_buffer` * `st_union` * `st_sample` * `st_convex_hull` * `st_voronoi` * ... --- ## All `sf` methods ``` ## [1] [ [[<- $<- ## [4] aggregate anti_join arrange ## [7] as.data.frame cbind coerce ## [10] dbDataType dbWriteTable distinct ## [13] dplyr_reconstruct filter full_join ## [16] gather group_by group_split ## [19] identify initialize inner_join ## [22] left_join mapView merge ## [25] mutate nest pivot_longer ## [28] plot print rbind ## [31] rename right_join rowwise ## [34] sample_frac sample_n select ## [37] semi_join separate_rows separate ## [40] show slice slotsFromS3 ## [43] spread st_agr st_agr<- ## [46] st_area st_as_s2 st_as_sf ## [49] st_bbox st_boundary st_buffer ## [52] st_cast st_centroid st_collection_extract ## [55] st_convex_hull st_coordinates st_crop ## [58] st_crs st_crs<- st_difference ## [61] st_filter st_geometry st_geometry<- ## [64] st_inscribed_circle st_interpolate_aw st_intersection ## [67] st_intersects st_is_valid st_is ## [70] st_join st_line_merge st_m_range ## [73] st_make_valid st_nearest_points st_node ## [76] st_normalize st_point_on_surface st_polygonize ## [79] st_precision st_reverse st_sample ## [82] st_segmentize st_set_precision st_shift_longitude ## [85] st_simplify st_snap st_sym_difference ## [88] st_transform st_triangulate st_union ## [91] st_voronoi st_wrap_dateline st_write ## [94] st_z_range st_zm summarise ## [97] transform transmute ungroup ## [100] unite unnest ## see '?methods' for accessing help and source code ``` --- class: inverse, middle, center # To go further --- # To dive even deeper into sf * Detailed sf package [vignettes](https://r-spatial.github.io/sf/articles/) * Blog posts: [here](https://www.r-spatial.org/r/2016/02/15/simple-features-for-r.html), [here](https://www.r-spatial.org/r/2016/07/18/sf2.html), [here](https://www.r-spatial.org/r/2016/11/02/sfcran.html), [here](https://www.r-spatial.org/r/2017/01/12/newssf.html) and [there](https://statnmap.com/fr/2018-07-14-initiation-a-la-cartographie-avec-sf-et-compagnie/) (in French) * [wiki page](https://github.com/r-spatial/sf/wiki/Migrating) describing sp-sf migration * Awesome online book [Geocomputation with R](https://geocompr.robinlovelace.net/) by Lovelace, Nowosad and Muenchow <img src="img/lovelacebookcover.png" width="150px" style="display: block; margin: auto;" /> --- # The [RStudio Cheat Sheets](https://www.rstudio.com/resources/cheatsheets/) <img src="img/sf_Page_1.png" width="600px" style="display: block; margin: auto;" /> --- # The [RStudio Cheat Sheets](https://www.rstudio.com/resources/cheatsheets/) <img src="img/sf_Page_2.png" width="600px" style="display: block; margin: auto;" /> --- class: title-slide-final, middle background-size: 55px background-position: 9% 15% # Thanks! ### I created these slides with [xaringan](https://github.com/yihui/xaringan) and [RMarkdown](https://rmarkdown.rstudio.com/) using the [rutgers css](https://github.com/jvcasillas/ru_xaringan) that I slightly modified. ### Credits: @cyber_nard contributed a few slides, and I used material from @StrimasMackey, @jafflerbach, @StatnMap, @SharpSightLabs and @edzerpebesma