owned this note
owned this note
Published
Linked with GitHub
# INBO CODING CLUB
26 February, 2019
Welcome
## Share your code snippet
If you want to share your code snippet, copy paste your snippet within a section of three backticks (```):
As an **example**:
```r
library(tidyverse)
...
```
(*you can copy paste this example and add your code further down, but do not fill in your code in this section*)
Your snippets:library
```r
```
## Problems loading `sf` and `tidyverse` in same session
Everybody today had problems loading correctly `sf` and `tidyverse`. We all got this:
```r
> library(sf)
> library(tidyverse)
Error in get(genname, envir = envir) : object 'group_map' not found
In addition: Warning message:
package ‘tidyverse’ was built under R version 3.5.2
Error in get(genname, envir = envir) : object 'group_split' not found
-- Attaching packages --------------------------------------- tidyverse 1.2.1 --
v ggplot2 3.1.0 v purrr 0.3.0
v tibble 2.0.99.9000 v dplyr 0.7.8
v tidyr 0.8.2 v stringr 1.3.1
v readr 1.3.1 v forcats 0.3.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
Warning messages:
1: package ‘readr’ was built under R version 3.5.2
2: package ‘purrr’ was built under R version 3.5.2
```
Well, actually we are not the only ones: check [here](https://github.com/r-spatial/sf/issues/969#issuecomment-467367771)! And check the answer provided by the head developer of `sf`! Just update `dplyr` to version 0.8.0!
Restart R (Ctrl+Shift+10) and then:
```
devtools::install_github("tidyverse/dplyr")
```
if you load tidyverse now, you can notice the new version: 0.8.0.9002.
```r
> library(tidyverse)
-- Attaching packages --------------------------------------- tidyverse 1.2.1 --
v ggplot2 3.1.0 v purrr 0.3.0
v tibble 2.0.1.9001 v dplyr 0.8.0.9002
v tidyr 0.8.2 v stringr 1.3.1
v readr 1.3.1 v forcats 0.3.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
Warning messages:
1: package ‘tidyverse’ was built under R version 3.5.2
2: package ‘readr’ was built under R version 3.5.2
3: package ‘purrr’ was built under R version 3.5.2
```
Now, if you load `sf` no errors anymore occur:
```r
> library(sf)
Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
Warning message:
package ‘sf’ was built under R version 3.5.2
```
## Alternative way to specify file paths
```r
filename <- "mydata.csv"
paste0("../data/", filename)
file.path("../data", filename)
file.path("..", "data", filename)
```
```r
library(here)
filename <- here("data", "mydata.csv")
```
## Trick to run unzip only when needed
```r
# file.exists checks if a file exists
# ! is the negation
if (!file.exists("../data/be_10km.shp")) {
# this is only run when "be_10km.shp" doesn't exists in "../data"
unzip("../data/20190226_utm10_bel.zip", exdir = "../data")
}
```
## Challenge 1
### Peter
1. Load EEA reference grid:
```r
grid <- st_read("../data/20190226_utm10_bel", "be_10km.shp")
```
2. Create an sf data.frame from crab data (= adding a geometry column to original data.frame):
```r
spatial_crab_df <- st_as_sf(
crab_df,
coords = c("decimalLongitude", "decimalLatitude"),
crs = 4326
)
```
Note that `x` = longitude (horizontal on map) and `y` = latitude (vertical on map). Don't forget to define the CRS (`4326` for WGS84).
3. Compare CRS:
```r
st_crs(grid)
st_crs(spatial_crab_df)
```
### Floris
```r
# 1. read shapefile UTM
utm10bel <- st_read("../data/20190226_utm10_bel/be_10km.shp")
# 2. set CRS of crab_df
str(crab_df, give.attr = FALSE)
class(crab_df)
spatial_crab_df <- st_as_sf(crab_df,
coords = c("decimalLongitude", "decimalLatitude"))
spatial_crab_df
st_crs(spatial_crab_df) <- 4326
# or at once:
# spatial_crab_df <- st_as_sf(crab_df,
# coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)
spatial_crab_df
class(spatial_crab_df)
spatial_crab_df$geometry
# 3. compare crs
st_crs(spatial_crab_df)
st_crs(utm10bel)
```
## Challenge 2
### floris
```r
# 1. How many occurrences of Chinese mitten crab in each cell SINCE 2015?
crabcounts <-
utm10bel %>%
st_join(spatial_crab_df %>%
# hint: use the CRS of the POLYGONS!
st_transform(crs = 3035),
left = FALSE) %>%
filter(year >= 2015) %>%
group_by(CELLCODE) %>%
summarise(number_occ = n()) %>%
st_drop_geometry
crabcounts
# 2. Find the 5 cells with the highest number of occurrences.
crabcounts %>%
arrange(desc(number_occ)) %>%
slice(1:5)
## Alternative
crabcounts %>%
top_n(5, number_occ)
```
Spatial joins:
st_join(sf-object1, sf-object2, join = )
- "join" argument is used for the topological operator
- default topological operator used by st_join() is st_intersects()
- when more than one row of the 2nd object topologically match the 1st, then rows of the 1st are duplicated!
Resource:
- https://geocompr.robinlovelace.net/
- https://geocompr.robinlovelace.net/spatial-operations.html#spatial-vec
### Damiano
Thanks Floris to show how to use the powerful and general function `st_join()`! Here below how to do it by using `st_intersection()` which also relies on `st_intersects()`:
```r
# get CRS of the grid
crs_belgium <- st_crs(belgium_grid)
# convert CRS of points to CRS of grid
spatial_crab_df <-
spatial_crab_df %>%
filter(year >= 2015) %>%
st_transform(crs_belgium)
# points on grid
assign_grid_to_pts <-
st_intersection(
x = spatial_crab_df,
y = belgium_grid)
# or via pipe
assign_grid_to_pts <-
spatial_crab_df %>%
st_intersection(eugrid_shape)
# counts pts in grid
n_pts_grid <-
assign_grid_to_pts %>%
group_by(CELLCODE) %>%
count() %>%
ungroup()
# get the 5 cells with the highest number of pts
n_pts_grid_top5 <-
n_pts_grid %>%
top_n(n = 5, wt = n)
n_pts_grid_top5
```
## Challenge 3
### Ivy
```r
# 1. Read areas from first layer ("ps_hbtrl")
st_layers("data/20190226_ps_hbtrl")
HabRL <- st_read("data/20190226_ps_hbtrl/ps_hbtrl.shp", crs = 31370)
# 2. Create a buffer of 1km around the protected areas
st_crs(HabRL)$units # in "m", dus 1km = 1000m
HabRL_buffer <- st_buffer(HabRL, dist = 1000)
# 3. How far is the nearest Chinese mitten crab for each area?
spatial_crab_df3 <- st_transform(spatial_crab_df, crs = 31370)
Afstand <- spatial_crab_df3 %>%
st_distance(HabRL_buffer) %>%
# distance matrix met punten in de rijen en areas in de kolommen
as.data.frame()
```
### Damiano
```r
# Get layers in folder 20190226_ps_hbtrl
st_layers("../data/20190226_ps_hbtrl")
# Read Habitatrichtlijngebieden from folder
habitats <- st_read("../data/20190226_ps_hbtrl", layer = "ps_hbtrl")
# you can read the right layer directly by passing the right file: ps_hbtrl.shp
habitats <- st_read("../data/20190226_ps_hbtrl/ps_hbtrl.shp")
# CRS not directly found (EPSG: NA) although the rest is present:
st_crs(habitats)
# assign CRS from metadata: Belge 1972 / Belgian Lambert 72 (EPSG:31370) found in the link
habitats <-
habitats %>%
st_transform(crs = 31370)
st_crs(habitats)
# check class of geomtry: point? polygon? multipolygon?
class(habitats$geometry)
# Create buffer of 1km (check the measurement unit you are using by checking the last
# part of message returned by st_crs(habitats)), or by running st_crs(habitats)$units. it will return "m", isn't?
st_crs(habitats)$units
habitats_buffer <-
habitats %>%
st_buffer(dist = 1000)
# Change CRS of crab data
spatial_crab_df_2015 <-
spatial_crab_df_2015 %>%
st_transform(crs = 31370)
# Get an idea what I am doing
ggplot() +
geom_sf(data = habitats_buffer) +
geom_sf(data = spatial_crab_df_2015)
# Measure distances from all pts to all polygons
distances <- st_distance(spatial_crab_df_2015, habitats_buffer, by_element = FALSE)
distances_df <- as_data_frame(distances)
# add the least distance of the crabs to each buffered area as column least_dist
habitats_buffer <- habitats_buffer %>%
mutate(least_dist = map_dbl(distances_df, min, na.rm = TRUE))
# get a look
View(habitats_buffer)
```
