Shapefiles in R using sf package

We will leave theory for some later posts and concentrate on a practical aspect for now. Our goal is to be able to create maps. To do that, we need to know the boundaries of areas i.e. countries, provinces or districts. That kind of data is usually stored in a special type of file, shapefile. This is clearly and oversimplification, but, as it is useful, lets stick to it for a moment.

There are a number of R packages to handle shapefiles. Most important are rgdal, sp; however we will focus on a relatively new addition - sf package that follows tidy data philosophy. Please look into the tutorial linked to get more info on how sf package works and learn a bit more about spatial data.

In our example we will use data on districts in Wroclaw. As mentioned before, data is in .shp file. Let’s read it into R and check what kind of object is it.

districts_wroclaw <- sf::read_sf('../../static/data/districts_wroclaw/districts_wroclaw.shp', 
                                 options = "ENCODING=UTF-8")
class(districts_wroclaw)

## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

head(districts_wroclaw)

## Simple feature collection with 6 features and 9 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 16.89252 ymin: 51.04515 xmax: 17.10723 ymax: 51.16599
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 10
##   OBJECTI NROSIED NAZWAOS  DATA   SHAPE_A SHAPE_L old_dst dstrct_  inhbtnt
##   <chr>     <int> <chr>    <chr>    <dbl>   <dbl> <chr>   <chr>      <dbl>
## 1 340          21 KrzykiP… 2016/… 5254965    9995 Krzyki  Krzyki-…  14.0  
## 2 341          24 GądówPo… 2016/… 3134569    7589 Fabryc… Gądów-P…  26.2  
## 3 348          42 Sołtyso… 2016/… 4547041    9004 Psie P… Sołtyso…   2.90 
## 4 349          18 Bieńkow… 2016/… 1433161    4760 Krzyki  Bieńkow…   0.500
## 5 351          32 Żerniki  2016/… 3908726    9133 Fabryc… Żerniki    3.40 
## 6 352          29 Muchobó… 2016/… 6879545   13228 Fabryc… Muchobó…   8.30 
## # ... with 1 more variable: geometry <POLYGON [°]>

Interesting! Apart from being special sf object, it is also a regular data.frame. Thanks to that, one can use packages like dplyr and tidyr directly with spatial data. This might not sounds that exciting, but it is a huge improvement compared to what was available in R year ago.

It is always nice, to be able to use old toolbox with new data. Luckily for us, plotting can be done using ggplot. Just make sure you use the new version from Github, and you can enjoy new geom_sf!

devtools::install_github('tidyverse/ggplot2')

Then everything becomes super easy.

ggplot() + 
  geom_sf(data = districts_wroclaw, aes(fill = inhbtnt), color='white') + 
  scale_fill_gradient(low = "white", high = "red")

Fill of districts represent number of inhabitants.

wroclaw_demo <- sf::read_sf('../../static/data/dem-rejurb-rejstat-shp/REJURB_20171231.shp')

wroclaw_demo_splited <- wroclaw_demo %>%
  mutate(below_44 = W0_2 +W3_6 +W7_12+ W13_15+ W16_18+ W19_24+ W25_34+ W35_44,
         above_44 = SUMA - below_44) %>%
  gather(group, value, below_44, above_44) %>%
  split(.$group)

install.packages(c('lwgeom'))

generate_samples <- function(data) 
  suppressMessages(st_sample(data, size = round(data$value / 100)))

points <- lapply(wroclaw_demo_splited, generate_samples)
points <- lapply(names(points), function(name){
  st_sf(data_frame(age = rep(name, length(points[[name]]))),
                      geometry = points[[name]])
})
points <- do.call(rbind, points)