R, GIS и fuzzyjoin: восстанавливаем статистические данные для регионов NUTS

    В этом посте речь пойдет о том, как я восстанавливал демографические данные для регионов Дании, где после реформы территориального устройства 2007 года официальной гармонизации данных не проводилось. Это лишь небольшая часть гармонизации евростатовских данных, которую я выполнил в рамках своего phd проекта. Пост сперва опубликован в моем англоязычном блоге и в блоге Demotrends. Думаю, что он может быть интересен далеко не только демографам.


    Что такое NUTS?


    NUTS расшифровывается как Nomenclature of Territorial Units For Statistics. Это стандартизированная система административно-территориального деления, принятая странами Евросоюза. История вопроса уходит в 1970-е, когда родилась идея сделать регионы различных стран Европы сопоставимыми. В более или менее законченном и широко употребимом виде система появилась лишь на рубеже веков. Существуют три основных уровня NUTS (см. рис. 1), и наиболее распространенным в региональном анализе оказывается NUTS-2.


    fig1
    Рисунок 1. Иллюстрация принципа выделения регионов NUTS различного иерархического уровня


    Главной задачей NUTS было создание сопоставимых критериев выделения регионов различного иерархического уровня. Тем не менее, в 2013 году численность населения регинов уровня NUTS-2 варьировалась от 28.5 тыс. чел. (Aland island, Finland) до почти 12 млн. чел (Ile-de-France, France — Париж и окрестности).


    Не сопоставимые ряды данных


    Довольно условное по своей сути, административно-территориальное деление не остается неизменным и постоянно трансформируется. Изменения территориальных систем создают значительные трудности для долгосрочного анализа процессов на региональном уровне, поскольку состыковка старых и новых регионов — трудоемкая и не всегда однозначно решаемая задача (если интересно, можно посмотреть про российские регионы в моей давней работе — приложения 1 и 2 в конце). Но несмотря на очевидные трудности для региональной статистики, границы изменяются постоянно в соответствии с требованиями государственных управленцев различного уровня. Евростат фиксирует все изменения, происходящие с регионами NUTS, и публикует детальные пояснения (рис. 2).


    fig2
    Рисунок 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” для регионов Дании.


    fig3
    Рисунок 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()

    fig4
    Рисунок 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))

    fig5
    Рисунок 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  <NA>  <NA>
    ## 65  504188.4 6269329   105319  <NA>  <NA>
    ## 195 533446.8 6312770    47071  <NA>  <NA>

    Всего 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))

    fig6
    Рисунок 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
    ##          <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

    Итак, для 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 <NA>   <NA>
    ## 2 nykøbing falsters    46349 <NA>   <NA>
    ## 3       herstederne    46101 <NA>   <NA>

    Поскольку таких случаев оказалось всего три, я исправил их вручную в пространственном датасете так, чтобы они соответствовали названиям в статистическом датасете. В двух случаях нестыковка произошла из-за неточностей в написании, а еще в одном случае произошла смена названия муниципалитета.


    # 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))

    fig7
    Рисунок 7. Доля населения в трудоспособном возрасте (15-64) в NUTS-3 регионах Дании в 2001


    Результат выгляди реалистично — показатель выше в столичном регионе и тех регионах, где есть сравнительно большие города.



    Поделиться публикацией
    Реклама помогает поддерживать и развивать наши сервисы

    Подробнее
    Реклама
    Комментарии 0

    Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.