R, GIS, and fuzzyjoin to reconstruct demographic data for NUTS regions of Denmark

NUTS

NUTS stands for the Nomenclature of Territorial Units For Statistics. The history of NUTS dates back to the beginning of 1970s, when European countries developed unified standards for systems of administrative geography. It was not until the beginning of this century when such a system finally became widely used. There are three main hierarchical levels of NUTS, and the most commonly used for regional analysis is NUTS-2.

Figure 1. Illustration of the principle of NUTS hierarchical system

One of the objectives of NUTS was to provide more or less comparable administrative divisions for all countries of Europe. Nevertheless, in 2013, population figures for single NUTS-2 regions ranged from 28.5 thousands in Aland island (Finland) to almost 12 million in Ile-de-France (Paris and surroundings, France).

The broken time series

Quite arbitrary in its essence, territorial division tends to evolve. Changes in administrative boundaries can cause problems for regional analysis as they break the time series and therefore make it harder to analyze trends. Despite this inconvenience, the boundaries of regions actually change quite often based on the needs and interests of local or national governmenta. Eurostat tracks all modifications providing detailed explanations of all the changes that happen between versions of NUTS (figure 2).

Figure 2. Changes in NUTS between versions 2006 and 2010

Despite this, Eurostat does not recalculate historic demographic data to match the most recent NUTS version. This means that, for the most recent version of NUTS, there is missing data for all years before the latest administrative change. So researchers have to reconstruct historical data manually to obtain a long time series. Of course, crude assumptions often have to be accepted in order to approximate the population figures for the current regions that did not exist in the past.

To make thing even more complex, Eurostat provides the data only for the latest version of NUTS (at least, I did not work out how to download previous versions). In my PhD project I carry out regional analysis for the NUTS-2 regions of European Union. To have the longest possible time series, when I did the data preparation in 2015, I chose the 2010 version of NUTS, on which the regional demographic projection EUROPOP2013 is based. For reproducibility, I uploaded the precise versions of the Eurostat data at NUTS-2 level on population age structures and deaths, as downloaded in 2015, to figshare.

Spatial matching

Let’s first have a look at the map.

ggplot()+
        geom_polygon(data = gd_nuts3, aes(long, lat, group = group), 
                     color = brbg[3], fill = "grey90", size = 1)+
        geom_point(data = mun_coord, aes(long, lat), 
                   color = brbg[10], size = 1)+
        theme_map()

Figure 4. Reference map of the old municipalities and NUTS-3 regions of Denmark

We can easily see that the boundaries of the municipalities (light blue) are much more precise than that of the NUTS-3 regions (orange/brown). This is not a problem as long as all the centroids of the municipalities fall within the boundaries of the NUTS-3 regions, which seems to be true for all municipalities except for the easternmost one. A quick check reveals that this is Christiansø, a tiny fortified island, whose history goes back to the Middle Ages. It has a special status and is not included into the NUTS system. For further manipulations, Christiansø can safely merge it with the close-by Bornholm.

To identify which municipalities fall into which NUTS regions, I use the spatial overlap function (over) from sp package. Here I should thank Roger Bivand, a person who made it possible to do any spatial analysis in R.

# municipality coordinates to Spatial
mun_centr <- SpatialPoints(coordinates(sp_mun), proj4string = CRS(proj4string(sp_nuts3)))

# spatial intersection with sp::over
inter <- bind_cols(mun_coord, over(mun_centr, sp_nuts3[,"NUTS_ID"])) %>% 
        transmute(long, lat, objectid,
                  nuts3 = as.character(NUTS_ID),
                  nuts2 = substr(nuts3, 1, 4))

Let’s again check visually if the spatial matching worked okay.

ggplot()+
        geom_polygon(data = gd_mun, aes(long, lat, group = group), 
                     color = brbg[9], fill = "grey90", size = .1)+
        geom_polygon(data = gd_nuts3, aes(long, lat, group = group), 
                     color = brbg[3], fill = NA, size = 1)+
        geom_point(data = inter, aes(long, lat, color = nuts3), size = 1)+
        geom_point(data = inter[is.na(inter$nuts3),], 
                   aes(long, lat), color = "red", size = 7, pch = 1)+
        theme_map(base_size = 15)+
        theme(legend.position = c(1, 1),
              legend.justification = c(1, 1))

Figure 5. Checking the spatial intersection between the old municipalities and NUTS-3 regions of Denmark

Not bad. But there is an “NA” category that represents all the cases where the spatial match failed. How many such cases do we have?

# how many failed cases do we have
sum(is.na(inter$nuts3))

## [1] 3

# where the intersection failed
inter[is.na(inter$nuts3),]

##         long     lat objectid nuts3 nuts2
## 23  892474.0 6147918    46399  <NA>  <NA>
## 65  504188.4 6269329   105319  <NA>  <NA>
## 195 533446.8 6312770    47071  <NA>  <NA>

As there are only 3 cases, I decided to fix them manually.

# fix the three cases manually
fixed <- inter
fixed[fixed$objectid=="46399", 4:5] <- c("DK014", "DK01")
fixed[fixed$objectid=="105319", 4:5] <- c("DK041", "DK04")
fixed[fixed$objectid=="47071", 4:5] <- c("DK050", "DK05")

The final visual check.

ggplot()+
        geom_polygon(data = gd_mun, aes(long, lat, group = group), 
                     color = brbg[9], fill = "grey90", size = .1)+
        geom_polygon(data = gd_nuts3, aes(long, lat, group = group), 
                     color = brbg[3], fill = NA, size = 1)+
        geom_point(data = fixed, aes(long, lat, color = nuts3), size = 1)+
        theme_map(base_size = 15)+
        theme(legend.position = c(1, 1),
              legend.justification = c(1, 1))

Figure 6. Re-checking the spatial intersection between the old municipalities and NUTS-3 regions of Denmark

Now everything seems okay.

Joining spatial and statistical data (fuzzy join)

The next task is to join the spatial data and statistical data together. The spatial layer for municipalities does not contain the codes that are used by Statistics Denmark, so I have to match municipalities in the two datasets by their names. This is quite a difficult task. Names can be written slightly differently, there are some special characters in Danish alphabet, and some municipalities may have experienced a change of name. To solve the task most efficiently, I used the ‘Fuzzy String Matching’ approach which is implemented in the fuzzyjoin package by David Robinson.

First, I simplify the names in both datasets turning them into lowercase, replacing the character “å” with “aa”, and removing the “Kommune” word in the spatial dataset names. Please note that I downloaded (separately) a small selection from Statistics Denmark to have a lightweight dataframe with municipal codes and names.

# simplify municipalities names
mun_geo <- mun_coord %>%         
        transmute(name = sub(x = name, " Kommune", replacement = ""), objectid) %>% 
        mutate(name = gsub(x = tolower(name), "å", "aa"))

mun_stat <- read.csv2("https://ikashnitsky.github.io/share/1703-nuts2-denmark/stat-codes-names.csv", 
                      fileEncoding = danish_ecnoding) %>% 
        select(name) %>% 
        separate(name, into = c("code", "name"), sep = " ", extra = "merge") %>% 
        mutate(name = gsub("\\s*\\([^\\)]+\\)", "", x = name)) %>% 
        mutate(name = gsub(x = tolower(name), "å", "aa"))

Let’s try fuzzy join.

# first attempt
fuz_joined_1 <- regex_left_join(mun_geo, mun_stat, by = "name")

The resulting dataframe has 278 rows instead of 271. That means that for some municipalities in the spatial dataset there was more than one match. Let’s identify them.

# identify more that 1 match (7 cases) and select which to drop
fuz_joined_1 %>% group_by(objectid) %>% mutate(n = n()) %>% filter(n > 1)

## Source: local data frame [14 x 5]
## Groups: objectid [7]
## 
##         name.x objectid  code      name.y     n
##          <chr>    <dbl> <chr>       <chr> <int>
## 1       haslev   105112   313      haslev     2
## 2       haslev   105112   403       hasle     2
## 3  brønderslev    47003   739       rønde     2
## 4  brønderslev    47003   805 brønderslev     2
## 5    hirtshals    47037   817        hals     2
## 6    hirtshals    47037   819   hirtshals     2
## 7      rønnede    46378   385     rønnede     2
## 8      rønnede    46378   407       rønne     2
## 9     hvidebæk    46268   317    hvidebæk     2
## 10    hvidebæk    46268   681     videbæk     2
## 11    ryslinge    46463   477    ryslinge     2
## 12    ryslinge    46463   737          ry     2
## 13     aarslev    46494   497     aarslev     2
## 14     aarslev    46494   861        aars     2

So, for 7 municipalities, two matches were found. I will drop the imperfect match variants in the next iteration of fuzzy join.

The other issue is the municipalities for which no match was found in that statistical data.

# show the non-matched cases
fuz_joined_1 %>% filter(is.na(code))

##              name.x objectid code name.y
## 1              faxe   105120 <NA>   <NA>
## 2 nykøbing falsters    46349 <NA>   <NA>
## 3       herstederne    46101 <NA>   <NA>

As there are only three such cases, I corrected them manually in the spatial data to match the statistical data. There are two cases of a difference in the way the name of municipality are written and one case of name change.

# correct the 3 non-matching geo names
mun_geo_cor <- mun_geo

mun_geo_cor[mun_geo_cor$name=="faxe", "name"] <- "fakse"
mun_geo_cor[mun_geo_cor$name=="nykøbing falsters", "name"] <- "nykøbing f."
mun_geo_cor[mun_geo_cor$name=="herstederne", "name"] <- "albertslund"

Now the second attempt to match the datasets (spatial dataset is corrected).

# second attempt
fuz_joined_2 <- regex_left_join(mun_geo_cor, mun_stat, by = "name")

# drop non-perfect match
fuz_joined_2 <- fuz_joined_2 %>% 
        group_by(objectid) %>% 
        mutate(n = n()) %>% 
        ungroup() %>% 
        filter(n < 2 | name.x==name.y)

fuz_joined_2 <- fuz_joined_2 %>% transmute(name = name.x, objectid, code)

The output looks perfect. Now, the last step – using the matched “objectid” field, I will finally attach the NUTS data to statistical codes.

# finally, attach the NUTS info to matched table
key <- left_join(fuz_joined_2, fixed, "objectid")

Aggregate old municipal data to NUTS levels

The previous manipulations yielded a dataframe that links statistical codes of the old municipalities with the corresponding NUTS regions. The last thing that has to be done is aggregation. I will attach the “key” dataset to a statistical dataset and aggregate the data at NUTS-3 and NUTS-2 levels.

# finally, we only need to aggregate the old stat data
df_agr <- left_join(key, df, "code") %>% 
        filter(!is.na(name)) %>% 
        gather("year", "value", y2001:y2006)

df_nuts3 <- df_agr %>% 
        group_by(year, sex, age, nuts3) %>% 
        summarise(value = sum(value)) %>% 
        ungroup()

df_nuts2 <- df_agr %>% 
        group_by(year, sex, age, nuts2) %>% 
        summarise(value = sum(value)) %>% 
        ungroup()

Let’s now calculate the shares of working age population in Danish NUTS-3 regions in 2001 and map the information.

# total population in 2001 by NUTS-3 regions
tot_01 <- df_nuts3 %>% 
        filter(year=="y2001") %>% 
        group_by(nuts3) %>% 
        summarise(tot = sum(value, na.rm = TRUE)) %>% 
        ungroup()

# working-age population in 2001 by NUTS-3 regions
working_01 <- df_nuts3 %>% 
        filter(year=="y2001", age %in% paste0("a0", 15:64)) %>% 
        group_by(nuts3) %>% 
        summarise(work = sum(value, na.rm = TRUE)) %>% 
        ungroup()

# calculate the shares of working age population
sw_01 <- left_join(working_01, tot_01, "nuts3") %>% 
        mutate(sw = work / tot * 100)

# map the shares of working age population in 2001 by NUTS-3 regions
ggplot()+
        geom_polygon(data = gd_nuts3 %>% left_join(sw_01, c("id" = "nuts3")),
                     aes(long, lat, group = group, fill = sw), 
                     color = "grey50", size = 1) +
        scale_fill_viridis()+
        theme_map(base_size = 15)+
        theme(legend.position = c(1, 1),
              legend.justification = c(1, 1))

Figure 7. The share of working age (15-64) population by NUTS-3 regions of Denmark in 2001

The result (thankfully!) looks realistic, with higher shares of the working-age population in the capital region, and in other regions that have relatively big cities.