R, GIS и fuzzyjoin: восстанавливаем демографические данные для NUTS регионов Дании
В этом посте речь пойдет о том, как я восстанавливал демографические данные для регионов Дании, где после реформы территориального устройства 2007 года официальной гармонизации данных не проводилось. Это лишь небольшая часть гармонизации евростатовских данных, которую я выполнил в рамках своего phd проекта. Пост сперва опубликован в моем англоязычном блоге и в блоге Demotrends. Думаю, что он может быть интересен далеко не только демографам.
Что такое NUTS?
NUTS расшифровывается как Nomenclature of Territorial Units For Statistics. Это стандартизированная система административно-территориального деления, принятая странами Евросоюза. История вопроса уходит в 1970-е, когда родилась идея сделать регионы различных стран Европы сопоставимыми. В более или менее законченном и широко употребимом виде система появилась лишь на рубеже веков. Существуют три основных уровня NUTS (см. рис. 1), и наиболее распространенным в региональном анализе оказывается NUTS-2.
Рисунок 1. Иллюстрация принципа выделения регионов NUTS различного иерархического уровня
Главной задачей NUTS было создание сопоставимых критериев выделения регионов различного иерархического уровня. Тем не менее, в 2013 году численность населения регинов уровня NUTS-2 варьировалась от 28.5 тыс. чел. (Aland island, Finland) до почти 12 млн. чел (Ile-de-France, France — Париж и окрестности).
Не сопоставимые ряды данных
Довольно условное по своей сути, административно-территориальное деление не остается неизменным и постоянно трансформируется. Изменения территориальных систем создают значительные трудности для долгосрочного анализа процессов на региональном уровне, поскольку состыковка старых и новых регионов — трудоемкая и не всегда однозначно решаемая задача (если интересно, можно посмотреть про российские регионы в моей давней работе — приложения 1 и 2 в конце). Но несмотря на очевидные трудности для региональной статистики, границы изменяются постоянно в соответствии с требованиями государственных управленцев различного уровня. Евростат фиксирует все изменения, происходящие с регионами NUTS, и публикует детальные пояснения (рис. 2).
Рисунок 2. Изменения в системе NUTS между версиями 2006 и 2010
Несмотря на это, Евростат далеко не всегда производит перерасчет, чтобы состыковать современное территориальное деление с историческими стат данными. На практике это означает, что для наиболее актуальных версий NUTS данные содержат значительное количество пропусков в тех случаях, когда границы регионов претерпевали изменения. Поэтому исследователям приходится самостоятельно гармонизировать данные, чтобы исследовать достаточно длинные временные периоды. Разумеется, процесс восстановления данных изобилует допущениями, иногда довольно грубыми, которые позволяют оценить статистические индикаторы для регионов, никогда не существовавших в прошлом.
Кроме того, работа осложняется тем, что Евростат публикует данные лишь для актуальной версии NUTS (по крайней мере я не нашел, где скачать архивные датасеты — буду благодарен за наводку, если они где-то все же существуют). В рамках своего PhD проекта я провожу региональный анализ на уровне NUTS-2. Чтобы обеспечить наиболее длинный период наблюдения, в 2015 году, когда я производил основную работу по гармонизации данных, Я выбрал версию NUTS 2010 года, на которой базировался официальный региональный демографический прогноз EUROPOP2013. В целях воспроизводимости результатов я выложил на figshare копии исходных данных на уровне NUTS о возрастных структурах населения и смертях, скачанные в 2015 году.
Дания
Некоторым странам пришлось произвести достаточно масштабные реформы административно-территориального деления, чтобы соответствовать стандартам NUTS. Наиболее значительные изменения произошли одномоментно в Дании в 2007 году, когда 271 старых муниципалитетов преобразовали в 98 новых (см. научную статью о реформе). Та же реформа и ввела NUTS для Дании: 98 новых муниципалитетов объединялись в 11 регионов NUTS-3, которые в свою очередь объединялись в 5 регионов NUTS-2. Уровень NUTS-1 не выделяется, что типично для небольших стран вроде Дании.
Насколько мне известно, попытки восстановить данные для NUTS регионов Дании до 2007 года на официальном уровне не предпринималось. Типичная карта регионов Европы, публикуемая Евростатом, для периода до 2007 года выглядит так (рис. 3) — «no data available» для регионов Дании.
Рисунок 3. Ожидаемая продолжительность предстоящей жизни в NUTS-2 регионах Европы, 2006; скриншот евростатовского interactive data exploratory tool
Подобное отсутствие данных крайне удивительно для столь развитой во многих отношениях страны как Дания. Может быть довольно сложно состыковать старую и новую системы муниципального деления, но агрегировать муниципальные данные на значительно более высоких иерархических уровнях NUTS не так уж и сложно. И это именно то, что я сделал и желаю показать в данном посте. Я потратил довольно много усилий, чтобы выяснить, не делал ли кто-то еще этого до меня — к большому моему удивлению, ничего в открытом доступе я не нашел.
Задача, по сути, сводится к тому, чтобы соотнести старые муниципалитеты (271) с современными NUTS-3 регионами (11) и дальше просто агрегировать данные на уровне NUTS-3. Далее NUTS-3 элементарно агрегируется до NUTS-2. Подобная задача может занять долгие вечера, особенно если учесть всю радость от возни с очаровательными датскими названиями муниципалитетов. Но, к счастью, мы живем в эпоху GIS. Я использовал GIS, чтобы почти автоматически состыковать старые муниципалитеты с регионами NUTS-3. Далее весь процесс показан пошагово с кодом R.
Данные
Данные о возрастной структуре старых 271 муниципалитетов Дании взяты с официального сайта Statistics Denmark. Система позволяет скачать разом лишь 10К записей для незарегистрированных пользователей и 100К после регистрации. Учитывая детальность наших данных, процесс выгрузки данных представляет собой довольно нудную работу. Поэтому для задач своего исследования я скачал данные за 2001–2006 гг. При необходимости детальные муниципальные данные доступны аж с 1979 года. Данные, скачанные мной в 2015 году и немного отформатированные перед использованием лежат здесь.
Поиск файла с пространственными данными границ старых муниципалитетов Дании занял значительное время. Вернувшись к этому вопросу полтора года спустя, я не сумел найти исходный источник шейпфайла. Но я почти уверен, что скачал его отсюда. Сегодня там визит уведомление, что данные не доступны из-за плановых работ. Копия использованного шейпфайла здесь.
Наконец, для успешного воплощения задуманного нам нужны границы регионов NUTS-3. Тут все просто — данные доступны на сайте Евростата (Eurostat geodata repository). Использованный мной заархивированный шейпфайл называется «NUTS_2010_20M_SH.zip». Выборка 11 датских регионов лежит тут.
Географическая проекция использованная для обоих шейпфайлов — ESPG-3044, довольно стандартная для Дании.
Наконец, R код, готовящий сессию и качающий данные (комментарии по коду не перевожу).
# set locale and encoding parameters to read Danish
if(Sys.info()['sysname']=="Linux"){
Sys.setlocale("LC_CTYPE", "da_DK.utf8")
danish_ecnoding <- "WINDOWS-1252"
}else if(Sys.info()['sysname']=="Windows"){
Sys.setlocale("LC_CTYPE", "danish")
danish_ecnoding <- "Danish_Denmark.1252"
}
# load required packages (install first if needed)
library(tidyverse) # version: 1.0.0
library(ggthemes) # version: 3.3.0
library(rgdal) # version: 1.2-4
library(rgeos) # version: 0.3-21
library(RColorBrewer) # version: 1.1-2
mypal <- brewer.pal(11, "BrBG")
library(fuzzyjoin) # version: 0.1.2
library(viridis) # version: 0.3.4
# load Denmark pop structures for the old municipalities
df <- read_csv("https://ikashnitsky.github.io/doc/misc/nuts2-denmark/BEF1A.csv.gz")
# create a directory for geodata
ifelse(!dir.exists("geodata"), dir.create("geodata"), "Directory already exists")
# download, unzip and read Danish NUTS-3 geodata (31KB)
url_nuts <- "https://ikashnitsky.github.io/doc/misc/nuts2-denmark/denmark-nuts3-espg3044.tgz"
path_nuts <- "geodata/denmark-nuts3-espg3044.tgz"
ifelse(!file.exists(path_nuts), download.file(url_nuts, path_nuts, mode="wb"), 'file alredy exists')
# If there are problems downloading the data automatically, please download it manually from
# https://ikashnitsky.github.io/doc/misc/nuts2-denmark/denmark-nuts3-espg3044.tgz
untar(tarfile = path_nuts, exdir = "geodata")
sp_nuts3 <- readOGR(dsn = "geodata/.", layer = "denmark-nuts3-espg3044")
gd_nuts3 <- fortify(sp_nuts3, region = "NUTS_ID") # to the ggplot format
# download, unzip and read Danish old municipal geodata (6.0MB)
url_mun <- "https://ikashnitsky.github.io/doc/misc/nuts2-denmark/kommune2006win1252.tgz"
path_mun <- "geodata/kommune2006win1252.tgz"
ifelse(!file.exists(path_mun), download.file(url_mun, path_mun, mode="wb"), 'file alredy exists')
# If there are problems downloading the data automatically, please download it manually from
# https://ikashnitsky.github.io/doc/misc/nuts2-denmark/kommune2006utf8.tgz
untar(tarfile = path_mun, exdir = "geodata")
sp_mun <- readOGR(dsn = "geodata/.", layer = "kommune2006win1252", encoding = danish_ecnoding)
gd_mun <- fortify(sp_mun)
# coordinates of the municipalities
mun_coord <- bind_cols(as.data.frame(coordinates(sp_mun)), sp_mun@data[,1:3]) %>%
transmute(long = V1, lat = V2, enhedid, objectid, name = navn)Пространственная состыковка
Давайте сперва глянем на карту.
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()
Рисунок 4. Старые муниципалитеты и NUTS-3 регионы Дании
Видно, что границы муниципалитетов (светло-голубые) гораздо более детальны, нежели границы NUTS-3 регионов (светло-коричневые). Это не представляет трудностей до тех пор, пока центроиды муниципалитетов попадают в пределы соответствующих регионов NUTS-3. Это справедливо для всех центроидов за исключением самой восточной точки. Быстрый гуглеж говорит нам, что это Christiansø, крохотный островок со средневековой крепостью и богатейшей историей. Этот муниципалитет обладает особым статусом и не вкючен в систему NUTS. Для дальнейших наших манипуляций можем смело «прилепить» его к соседнему острову Bornholm.
Для выяснения, какие из центроидов попадают в какие регионы NUTS, я воспользовался функцией пространственного пересечения (over) из библиотеки sp. Тут я хочу принести благодарность Roger Bivand, человеку, без которого не было бы современных картографических методов для 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))Проверим, работает ли пространственная состыковка.
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))
Рисунок 5. Проверка пространственного пересечения между старыми муниципалитетами и регионами NUTS-3 Дании
Не плохо. Но есть несколько муниципалитетов, которые попали в категорию «NA», что означает, что для них не удалось найти соответствия на уровне NUTS. Сколько у нас таких случаев?
# 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
## 65 504188.4 6269329 105319
## 195 533446.8 6312770 47071 Всего 3. Я решил пофиксить их вручную.
# 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")Финальная визуальная проверка.
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))
Рисунок 6. Перепроверка пространственной состыковки после ручной корректировки для 3 муниципалитетов
Теперь все в норме.
Соединяем пространственные и статистические данные (fuzzy join)
Следующей задачей оказалась не тривиальная состыковка пространственных и статистических данных. Найденный мной шейпфайл муниципалитетов не содержал идентификационных кодов, используемых Statistics Denmark. Поэтому пришлось состыковывать муниципалитеты в пространственном и статистическом датасетах по названию — со всеми прелестями датских имен и западно-европейской кодировки. Названия муниципалитетов могут быть записаны чуть-чуть по-разному в двух датасетах. А муниципалитетов сотни. Казалось бы, достаточно скучная рутинная работа. Не тут-то было! На помощь приходит «Fuzzy String Matching» и его прекрасная имплементация в библиотеке fuzzyjoin, написанной замечательным разработчиком David Robinson.
Для начала я упростил данные имен муниципалитетов, переведя их в нижний регистр, разумеется, в обоих датасетах. Кроме того, первые попытки стыковки показали, что имеет смысл заменить букву «å» на альтернативное написание — «aa». В пространственном датасете я также убрал слово «Kommune» из имени каждого муниципалитета. Я скачал крохотный датасет с сайта Statistics Denmark, содержащий имена и коды старых муниципалитетов.
# 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/doc/misc/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"))Пробуем fuzzy join.
# first attempt
fuz_joined_1 <- regex_left_join(mun_geo, mun_stat, by = "name")На выходе получаем датафрейм с 278 строчками вместо 271. Это означает, что для некоторых муниципалитетов из пространственного датасета было найдено более одного соответствия в статистическом датасете. Найдем эти проблемные случаи.
# 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
##
## 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 Итак, для 7 муниципалитетов нашлось по 2 соответствия. Не идеальные соответствия мы выбросим при следующей итерации fuzzy join.
Другая проблема в том, что для некоторых муниципалитетов вообще не нашлось соответствия.
# show the non-matched cases
fuz_joined_1 %>% filter(is.na(code))## name.x objectid code name.y
## 1 faxe 105120
## 2 nykøbing falsters 46349
## 3 herstederne 46101 Поскольку таких случаев оказалось всего три, я исправил их вручную в пространственном датасете так, чтобы они соответствовали названиям в статистическом датасете. В двух случаях нестыковка произошла из-за неточностей в написании, а еще в одном случае произошла смена названия муниципалитета.
# 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"Теперь попробуем состыковать во второй раз (названия в пространственном датасете откорректированы).
# 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)Идеально. Теперь последний шаг — по полю «objectid» прикрепим данные о принадлежности регионам NUTS к старым статистическим кодам муниципалитетов.
# finally, attach the NUTS info to matched table
key <- left_join(fuz_joined_2, fixed, "objectid")Агрегирование муниципальных данных на уровнях NUTS
Предыдущие манипуляции дали нам небольшой, но очень ценный датафрейм, который состыковывает статистическое кода старых муниципалитетов с современными регионами NUTS. Нам осталось лишь агрегировать старые данные. И прикреплю «key» датасет к статистическим данным и агрегирую из на уровнях NUTS-3 и NUTS-2.
# 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()Давайте для пробы посчитаем долю населения в трудоспособном возрасте (15–64) в регионах NUTS-3 Дании в 2001 году и кинем эту информацию на карту.
# 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))
Рисунок 7. Доля населения в трудоспособном возрасте (15–64) в NUTS-3 регионах Дании в 2001
Результат выгляди реалистично — показатель выше в столичном регионе и тех регионах, где есть сравнительно большие города.
