Joining attribute data with geofi data

Markus Kainu

2021-03-29

This vignettes provides few examples on how to join attribute data from common sources of attribute data. Here we are using data from three different sources of which two are from Statistics Finland PxWeb-api, namely municipality key figures and Paavo (Open data by postal code area). Third is Covid-19 tests and confirmed cases in Finland.

Installation

geofi can be installed from CRAN using

# install from CRAN
install.packages("geofi")

# Install development version from GitHub
remotes::install_github("ropengov/geofi")

Municipalities

Municipality data provided by get_municipalities()-function contains 77 indicators variables from each of 309 municipalities. Variables can be used either for aggregating data or as keys for joining attribute data.

PxWeb-data

In this first example we join municipality level data from Statistics Finland municipality key figures

library(geofi)
muni <- get_municipalities(year = 2019)

library(pxweb)
pxweb_query_list <-
  list("Alue 2020"=c("*"),
       "Tiedot"=c("*"),
       "Vuosi"=c("2019"))

px_raw <-
  pxweb_get(url = "https://pxnet2.stat.fi/PXWeb/api/v1/en/Kuntien_avainluvut/2020/kuntien_avainluvut_2020_aikasarja.px",
            query = pxweb_query_list)

library(dplyr)
library(janitor)
library(sf)
px_data <- as_tibble(
  as.data.frame(px_raw, 
                column.name.type = "text", 
                variable.value.type = "text")
  ) %>% setNames(make_clean_names(names(.)))
px_data
#> # A tibble: 12,800 x 4
#>    region_2020   information                           year  municipal_key_figu…
#>    <chr>         <chr>                                 <chr>               <dbl>
#>  1 WHOLE COUNTRY Degree of urbanisation, %             2019                 86.4
#>  2 WHOLE COUNTRY Population                            2019            5525292  
#>  3 WHOLE COUNTRY Population change from the previous … 2019                  0.1
#>  4 WHOLE COUNTRY Share of persons aged under 15 of th… 2019                 15.8
#>  5 WHOLE COUNTRY Share of persons aged 15 to 64 of th… 2019                 62  
#>  6 WHOLE COUNTRY Share of persons aged over 64 of the… 2019                 22.3
#>  7 WHOLE COUNTRY Share of Swedish-speakers of the pop… 2019                  5.2
#>  8 WHOLE COUNTRY Share of foreign citizens of the pop… 2019                  4.8
#>  9 WHOLE COUNTRY Excess of births, persons             2019              -8336  
#> 10 WHOLE COUNTRY Intermunicipal migration gain/loss, … 2019                  0  
#> # … with 12,790 more rows

Once we have the data in long format we can observe the region_2020-column.

count(px_data, region_2020)
#> # A tibble: 400 x 2
#>    region_2020              n
#>    <chr>                <int>
#>  1 Akaa                    32
#>  2 Alajärvi                32
#>  3 Alavieska               32
#>  4 Alavus                  32
#>  5 Asikkala                32
#>  6 Askola                  32
#>  7 Aura                    32
#>  8 Brändö                  32
#>  9 Central Finland         32
#> 10 Central Ostrobothnia    32
#> # … with 390 more rows

This is not obvious to all, but have the municipality names in Finnish among other regional breakdowns which allows us to combine the data with spatial data using municipality_name_fi-variable.

map_data <- right_join(muni, px_data, by = c("municipality_name_fi" = "region_2020"))

Now we can plot a map showing Share of Swedish-speakers of the population, % and Share of foreign citizens of the population, % on two panels sharing a scale.

library(ggplot2)
map_data %>% 
  filter(grepl('Share of Swedish-speakers of the population|Share of foreign citizens of the population', information)) %>% 
  ggplot(aes(fill = municipal_key_figures)) + 
  geom_sf() + 
  facet_wrap(~information) +
  theme(legend.position = "top")

Health Districts

in early 2021 we are still troubled by the COVID-19 and the health authorities are counting infections, deaths and vaccinated. Lets pull the daily data from API and compare the names of the health districts both in COVID-19 data and in municipality division from Statistics Finland.

if (FALSE){
  library(readr)
  cols(
    Area = col_character(),
    Time = col_date(format = ""),
    val = col_double()
  ) -> cov_cols
  
  thl_korona_api <- "https://sampo.thl.fi/pivot/prod/en/epirapo/covid19case/fact_epirapo_covid19case.csv?row=dateweek20200101-508804L&column=hcdmunicipality2020-445222L"
  status <- httr::status_code(httr::GET(thl_korona_api))

  xdf_raw <- read_csv2(thl_korona_api, col_types = cov_cols)
  xdf <- xdf_raw %>% 
    # filter(!grepl("Kaikki", Alue)) %>% 
    rename(date = Time, 
           shp = Area, 
           day_cases = val) %>% 
    group_by(shp) %>% 
    arrange(shp,date) %>% 
    filter(!is.na(day_cases)) %>% 
    mutate(total_cases = cumsum(day_cases)) %>% 
    ungroup() %>% 
    group_by(shp) %>% 
    filter(date == max(date, na.rm = TRUE)) %>% 
    ungroup()
} 
xdf <- structure(list(shp = c("Åland", "All areas", "Central Finland Hospital District", 
"Central Ostrobothnia Hospital District", "Helsinki and Uusimaa Hospital District", 
"Itä-Savo Hospital District", "Kainuu Hospital District", "Kanta-Häme Hospital District", 
"Kymenlaakso Hospital District", "Länsi-Pohja Hospital District", 
"Lappi Hospital District", "North Karelia Hospital District", 
"North Ostrobothnia Hospital District", "North Savo Hospital District", 
"Päijät-Häme Hospital District", "Pirkanmaa Hospital District", 
"Satakunta Hospital District", "South Karelia Hospital District", 
"South Ostrobothnia Hospital District", "South Savo Hospital District", 
"Southwest Finland Hospital District", "Vaasa Hospital District"
), date = structure(c(18674, 18674, 18674, 18674, 18674, 18674, 
18674, 18674, 18674, 18674, 18674, 18674, 18674, 18674, 18674, 
18674, 18674, 18674, 18674, 18674, 18674, 18674), class = "Date"), 
    day_cases = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0), total_cases = c(120, 51047, 1850, 177, 
    29519, 141, 227, 873, 854, 478, 465, 516, 2265, 850, 1243, 
    2662, 671, 348, 534, 573, 4803, 1878)), row.names = c(NA, 
-22L), class = c("tbl_df", "tbl", "data.frame"))

xdf %>% 
    count(shp)
#> # A tibble: 22 x 2
#>    shp                                        n
#>    <chr>                                  <int>
#>  1 All areas                                  1
#>  2 Central Finland Hospital District          1
#>  3 Central Ostrobothnia Hospital District     1
#>  4 Helsinki and Uusimaa Hospital District     1
#>  5 Itä-Savo Hospital District                 1
#>  6 Kainuu Hospital District                   1
#>  7 Kanta-Häme Hospital District               1
#>  8 Kymenlaakso Hospital District              1
#>  9 Lappi Hospital District                    1
#> 10 Länsi-Pohja Hospital District              1
#> # … with 12 more rows
  
muni <- get_municipalities(year = 2021) 
muni %>% 
  st_drop_geometry() %>% 
  count(sairaanhoitop_name_en)
#>                     sairaanhoitop_name_en  n
#> 1       Central Finland Hospital District 21
#> 2  Central Ostrobothnia Hospital District 10
#> 3  Helsinki and Uusimaa Hospital District 24
#> 4              Itä-Savo Hospital District  4
#> 5                Kainuu Hospital District  8
#> 6            Kanta-Häme Hospital District 11
#> 7           Kymenlaakso Hospital District  6
#> 8                 Lappi Hospital District 15
#> 9           Länsi-Pohja Hospital District  6
#> 10        North Karelia Hospital District 13
#> 11   North Ostrobothnia Hospital District 29
#> 12           North Savo Hospital District 18
#> 13            Pirkanmaa Hospital District 23
#> 14          Päijät-Häme Hospital District 12
#> 15            Satakunta Hospital District 16
#> 16        South Karelia Hospital District  9
#> 17   South Ostrobothnia Hospital District 18
#> 18           South Savo Hospital District  9
#> 19    Southwest Finland Hospital District 28
#> 20                Vaasa Hospital District 13
#> 21                                  Åland 16

The names look identical so we can join the two datasets and plot a map.

muni %>% 
  count(sairaanhoitop_name_en) %>% 
  left_join(xdf, by = c("sairaanhoitop_name_en" = "shp")) %>% 
  ggplot(aes(fill = total_cases)) +
  geom_sf() +
  geom_sf_text(aes(label = paste0(sairaanhoitop_name_en, "\n", total_cases)), 
               color = "white") +
  labs(title = "Number of total COVID-19 cases reported since January 2020", 
       fill = NULL)

Zipcode level

You can download data from Paavo (Open data by postal code area) using pxweb-package in a similar manner as in the first example.

library(pxweb)
# lets get all zipcodes and all variables
pxweb_query_list <- 
  list("Postinumeroalue"=c("*"),
                           "Tiedot"=c("*"))

# Download data 
px_raw <- 
 pxweb_get(url = "https://pxnet2.stat.fi/PXWeb/api/v1/en/Postinumeroalueittainen_avoin_tieto/2019/paavo_1_he_2019.px",
           query = pxweb_query_list)

px_data <- as_tibble(
  as.data.frame(px_raw, 
                column.name.type = "text", 
                variable.value.type = "text")
  ) %>% setNames(make_clean_names(names(.)))
px_data %>% 
  filter(postal_code_area != "Finland")
#> # A tibble: 72,624 x 3
#>    postal_code_area               data                paavo_open_data_by_postal…
#>    <chr>                          <chr>                                    <dbl>
#>  1 00100 Helsinki Keskusta - Etu… Inhabitants, total…                      18284
#>  2 00100 Helsinki Keskusta - Etu… Females, 2017 (HE)                        9613
#>  3 00100 Helsinki Keskusta - Etu… Males, 2017 (HE)                          8671
#>  4 00100 Helsinki Keskusta - Etu… Average age of inh…                         41
#>  5 00100 Helsinki Keskusta - Etu… 0-2 years, 2017 (H…                        434
#>  6 00100 Helsinki Keskusta - Etu… 3-6 years, 2017 (H…                        521
#>  7 00100 Helsinki Keskusta - Etu… 7-12 years, 2017 (…                        711
#>  8 00100 Helsinki Keskusta - Etu… 13-15 years, 2017 …                        274
#>  9 00100 Helsinki Keskusta - Etu… 16-17 years, 2017 …                        185
#> 10 00100 Helsinki Keskusta - Etu… 18-19 years, 2017 …                        264
#> # … with 72,614 more rows

Before we can join the data, we must extract the numerical postal code from postal_code_area-variable.

px_data$posti_alue <- sub(" .+$", "", px_data$postal_code_area)

# Lets join with spatial data and plot the area of each zipcode
zipcodes19 <- get_zipcodes(year = 2019) 
zipcodes_map <- left_join(zipcodes19, 
                          px_data %>% filter(data == "Average age of inhabitants, 2017 (HE)"))
ggplot(zipcodes_map) + 
  geom_sf(aes(fill = paavo_open_data_by_postal_code_area_2019), 
          color  = alpha("white", 1/3)) +
  labs(title = "Average age of inhabitants, 2017 (HE)", 
       fill = NULL)