Visualize a choropleth map of MSD prevalence in EU, using ggplot and the Eurostat API for data and administrative boundaries.

20220113_eurostat_hsw_pb5_subset.csv

vis_map_eurostat_msd.r

# ######################################################################### #
# Visualize a choropleth map of musculo-skeletal disorder prevalence in EU, #
#  using geom_sf() and Eurostat API for data and administrative boundaries. #
# ######################################################################### #

library(eurostat)
library(lubridate)
library(dplyr)
library(tidyr)
library(countrycode)
library(giscoR)
library(sf)
library(janitor)
library(ggplot2)
library(hrbrthemes)
library(readr)
library(assertr)
library(here)


# Query Eurostat for work-related health datasets
search_eurostat('work-related health')
# ...
# 6 Persons reporting a work-related health problem resulting in limitations of daily activities              hsw_pb4  datas… 12.10.2021 ...
# 7 Persons reporting a work-related health problem by sex, age and type of problem                           hsw_pb5  datas… 12.10.2021 ...
# ...

# query Eurostat for hsw_pb5: Persons reporting a work-related health problem by sex, age and type of problem
hsw_pb5 <- get_eurostat(
  id = 'hsw_pb5', 
  type = 'label', # get informative labels instead of (shorter) codes
  time_format = "date") %>%
  clean_names() # clean/standardise names

# subset data: only MSD and population totals
hsw_pb5_subset <- hsw_pb5 %>%
  filter(
    sex == 'Total', age == 'From 15 to 64 years', 
    diagnose == 'Musculo-skeletal disorders',
    time == '2020-01-01',
    unit == 'Percentage of persons reporting a health problem',
    !(geo %in% c('European Union - 27 countries (from 2020)', 
                 'Euro area - 19 countries  (from 2015)', 
                 'European Union - 28 countries (2013-2020)') )) %>%
  mutate(geo = if_else(geo == 'Germany (until 1990 former territory of the FRG)', 'Germany', geo)) %>%
  # when using labels, you dont have the ISO-countrycode -> add for merging with spatial boundaries
  mutate(iso3_code = countrycode(
    geo, origin = 'country.name', destination = 'iso3c')) %>%
  arrange(desc(values))


# export dataset for to reproduce data from API later if needed
hsw_pb5_subset %>%
  verify(dim(.) == c(30, 8) )%>%
  write_csv(here::here('20220113_eurostat_hsw_pb5_subset.csv'))


# get list of ISO 3-digit country codes for mainland EU-countries using countrycode-package
eu_mainland_countries_iso3c <- countrycode::codelist %>% 
  select(country.name.en, iso3c, region, continent) %>%
  filter(continent == 'Europe') %>%
  filter(iso3c != 'ALA') %>%  # Åland Islands not recognised by gisco_get_countries()
  filter(iso3c != 'RUS') %>% # # drop Russia
  pull(iso3c)

# Use list of ISO 3-digit country codes for mainland EU to query the Eurostat GISCO-service API 
#  for a sf dataframe with the administrative boundaries for these countries
eu_mainland <- gisco_get_countries(
  year = '2020', epsg = '4326', cache_dir = tempdir(),
  resolution = '20', 
  country = eu_mainland_countries_iso3c,
  spatialtype = 'RG') %>% # request region (multi)polygons
  clean_names() # clean/standardise names

# crop resulting multipolygons in sf dataframe to "mainland" (drop islands etc.)
eu_mainland <- st_crop(
  eu_full_noru, 
  c(xmin=-10, xmax=45, ymin=36, ymax=71))

# merge sf dataframe with Eurstat dataframe on ISO-countrycode
eu_mainland <- eu_mainland %>%
  left_join(hsw_pb5_subset, by = 'iso3_code')

# plot map using ggplot2::geom_sf()
p_eu_msa_prevalence <- ggplot(eu_mainland) +
  geom_sf(
    aes(fill = values), 
    color = 'white', lwd = 0.4) +  # thin, white country lines
  coord_sf(
    # limit map to 'mainland' EU
    xlim = c(2500000, 6500000), ylim =c(1500000, 5300000),
    label_axes = "",
    crs = 3035) +
  scale_fill_distiller(
    type = 'seq', palette = 'Reds', direction = 1,
    na.value="grey80",
    labels = scales::percent_format(scale = 1, accuracy = 1)) +
  theme_ipsum_rc() +
  theme(legend.position = c(.89, .85)) +
  labs(
    x = NULL, y = NULL, fill = NULL,
    title = "Musculo-skeletal disorders are by far the leading\nhealth problem for working Europeans",
    subtitle = "The percentage of working persons with a health problem,\nthat report a musculo-skeletal disorder (2020).",
    caption = 'Source: Eurostat/LFS (hsw_pb5, 2020), 15-64 year olds that are currently working or were working in\nthe past • Maarten Hermans (@hermansm) • © EuroGeographics for the administrative borders.')

ggsave(
  here::here('map_eu_msd_prevalence_en.png'), 
  p_eu_msa_prevalence, width = 7, height = 7)