diff --git a/cbs_2021/cbs_data_updaten.R b/cbs_2021/cbs_data_updaten.R index 1b761f1..d91dfec 100644 --- a/cbs_2021/cbs_data_updaten.R +++ b/cbs_2021/cbs_data_updaten.R @@ -1,233 +1,235 @@ -library(dplyr) -library(sf) -# De %like% functie, zodat je kunt gebruiken: "Test" %like% "te" (= TRUE) -`%like%` <- function (x, pattern) grepl(pattern = tolower(pattern), x = tolower(x), fixed = FALSE) - - -# Bronnen ----------------------------------------------------------------- - -downloadmap <- "/Users/msberends/Downloads/" - -# download inwoners (gescheiden op geslacht) per postcode hier als 'CSV volgens tabelindeling': -# https://opendata.cbs.nl/#/CBS/nl/dataset/83503NED/table?dl=42FC6 -# verwijder dan de eerste rijen en de laatste rij ("Bron: CBS") -postcodes_bestand <- paste0(downloadmap, "Bevolking__geslacht__migratieachtergrond__viercijferige_postcode__1_januari_15102020_095747.csv") -# download inwoners per 5 jaar leeftijd en geslacht hier voor het huidige jaar: -# https://opendata.cbs.nl/#/CBS/nl/dataset/83502NED/table?dl=42FE0 -inwoners_bestand <- paste0(downloadmap, "Bevolking__leeftijd__postcode_15102020_102334.csv") -# gebiedsindelingen hier downloaden: -# https://www.cbs.nl/nl-nl/dossier/nederland-regionaal/geografische-data/cbs-gebiedsindelingen -gebiedsindelingen_bestand <- paste0(downloadmap, "cbsgebiedsindelingen_2021_v1.gpkg") -# download postcodes 4 onder 'Downloads' ('Numeriek deel van de postcode (PC4)') hier: -# https://www.cbs.nl/nl-nl/dossier/nederland-regionaal/geografische-data/gegevens-per-postcode -postcodes4_bestand <- paste0(downloadmap, "2019-CBS_PC4_2018_v1/CBS_PC4_2018_v1.shp") - - -# Helpfuncties ------------------------------------------------------------ - -kaart_fixen <- function(kaart) { - # CRS opnieuw instellen, bij nieuwe GDAL geeft dit anders problemen: - st_crs(kaart) <- clean_numeric(st_crs(kaart)$input) - if (is.na(st_crs(kaart))) { - # 'reserve' CRS - st_crs(kaart) <- 28992 - } - if (!"geometry" %in% colnames(kaart)) { - geom_naam <- c("shape", "geom") - geom_naam <- colnames(kaart)[colnames(kaart) %in% geom_naam][1] - geometry <- kaart[, colnames(kaart)[tolower(colnames(kaart)) == geom_naam], drop = TRUE] - kaart <- st_drop_geometry(kaart) - kaart <- st_set_geometry(kaart, geometry) - } - # geometrie als breedte- en lengtegraad, CRS transformeren van 28992 naar 4326 - kaart <- st_transform(kaart, crs = 4326) - - # alle ongeldige vormen geldig maken - if (any(!st_is_valid(kaart$geometry))) { - kaart$geometry <- st_make_valid(kaart$geometry) - } - - # CRS opnieuw instellen, bij nieuwe GDAL geeft dit anders problemen: - st_crs(kaart) <- clean_numeric(st_crs(kaart)$input) - if (is.na(st_crs(kaart))) { - # 'reserve' CRS - st_crs(kaart) <- 28992 - } - # geometrie als breedte- en lengtegraad, CRS transformeren van 28992 naar 4326 - kaart <- st_transform(kaart, crs = 4326) - # oppervlakte toevoegen - kaart$area_km2 <- as.double(st_area(st_geometry(kaart)) / 1000 ^ 2) - kaart -} - -lagen_beschikbaar <- sort(st_layers(gebiedsindelingen_bestand)$name) -downloaden_transformeren <- function(laag) { - zoeklaag <- sort(lagen_beschikbaar[lagen_beschikbaar %like% laag & - lagen_beschikbaar %like% "gegeneraliseerd" & - !lagen_beschikbaar %like% "niet"]) - if (length(zoeklaag) == 0) { - stop("Geen laag gevonden") - } - zoeklaag <- zoeklaag[length(zoeklaag)] - message("Laag '", zoeklaag, "' gevonden", appendLF = FALSE) - - kaart <- st_read(gebiedsindelingen_bestand, layer = zoeklaag, quiet = TRUE) - kaart <- kaart_fixen(kaart) - - message(", met ", nrow(kaart), " geometrieën") - colnames(kaart)[colnames(kaart) %like% "naam"] <- laag - # geen factor - kaart[, laag] <- as.character(kaart[, laag, drop = TRUE]) - - # alleen het type, oppervlakte en de geometrie behouden - kaart <- kaart[, c(laag, "area_km2", "geometry"), drop = FALSE] - kaart -} - - -# Inwoners per postcode/leeftijd/geslacht --------------------------------- - -inwoners_per_postcode_leeftijd <- read_csv2(inwoners_bestand) -inwoners_per_postcode_leeftijd <- inwoners_per_postcode_leeftijd %>% - filter(!leeftijd %like% "totaal", postcode %like% "[0-9]") %>% - mutate(postcode = as.double(postcode), - geslacht = case_when(geslacht %like% "totaal" ~ "inwoners", - geslacht %like% "man" ~ "inwoners_man", - geslacht %like% "vrouw" ~ "inwoners_vrouw"), - # bij CBS vinden ze 0-5 en 5-10 handig. Waarin zit iemand dan die 5 is?! We maken er 0-4 en 5-9 van: - leeftijd_min = as.numeric(gsub("([0-9]+).*", "\\1", leeftijd)), - leeftijd_max = as.numeric(gsub("([0-9]+)[^0-9]+([0-9]+)[^0-9]+", "\\2", leeftijd)), - leeftijd_nieuw = paste0(leeftijd_min, "-", leeftijd_max - 1), - leeftijd_nieuw = gsub("95-NA", "95+", leeftijd_nieuw), - leeftijd = factor(leeftijd_nieuw, levels = levels(age_groups(0, 5 * c(1:19))), ordered = TRUE)) %>% - select(postcode, geslacht, leeftijd, inwoners) %>% - pivot_wider(names_from = geslacht, values_from = inwoners) -# alle PC2 en PC3 toevoegen -inwoners_per_postcode_leeftijd <- inwoners_per_postcode_leeftijd %>% - bind_rows(inwoners_per_postcode_leeftijd %>% - group_by(postcode = clean_numeric(substr(postcode, 1, 2)), leeftijd) %>% - summarise_all(sum)) %>% - bind_rows(inwoners_per_postcode_leeftijd %>% - group_by(postcode = clean_numeric(substr(postcode, 1, 3)), leeftijd) %>% - summarise_all(sum)) %>% - arrange(postcode, leeftijd) - - -# Postcodes (wordt later alle referentiedata aan toegevoegd) -------------- - -# `postcodes` is hier de vorige versie die we als `postcodes` gebruikten, deze wordt vernieuwd -postcodes_plaats_gemeente <- postcodes %>% - filter(postcode > 999) %>% # alleen PC4 houden, wordt later weer aangevuld met PC2 en PC3 - select(postcode, plaats, gemeente) - -postcodes <- read_csv2(postcodes_bestand) -colnames(postcodes) <- gsub("(man|vrouw)n?en", "\\1", colnames(postcodes)) - - -# Postcode-4 kaart -------------------------------------------------------- - -# we gebruiken deze postcodekaart om te bepalen welke postcodes in welk gebied liggen met sf::st_intersects(). -# dus de geometrie van postcode 9251 valt in het snijvlak van de de geometrie van de gemeente Tytsjerksteradiel -# en dus is Tytsjerksteradiel de gemeente van postcode 9251 (en zo verder voor NUTS-3, GGD-regio, ...) -postcodes4 <- st_read(postcodes4_bestand) -postcodes4 <- kaart_fixen(postcodes4) -# alleen relevante kolommen houden -postcodes4 <- postcodes4 %>% - transmute(postcode = as.double(as.character(PC4)), - huishoudens = ifelse(AANTAL_HH < 0, NA_real_, AANTAL_HH), - huishouden_grootte = ifelse(GEM_HH_GR < 0, NA_real_, GEM_HH_GR), - area_km2 = as.double(st_area(geometry) / 1000 ^ 2), - geometry) - - -# Referentiewaarden aan `postcodes` toevoegen en kaarten opslaan ---------- - -relevante_lagen <- c("gemeente", - "provincie", - "nuts3", - "ggdregio", - "jeugdregio", - "veiligheidsregio") -for (i in 3:length(relevante_lagen)) { - message(">> zoeken naar ", relevante_lagen[i]) - kaart <- downloaden_transformeren(relevante_lagen[i]) - - if (!relevante_lagen[i] %in% c("plaats", "gemeente")) { - # referentiedata toevoegen aan 'postcodes' - p <- dplyr::progress_estimated(length(postcodes4$geometry)) - newvar <- character(length = nrow(postcodes4)) - for (pc in 1:nrow(postcodes4)) { - p$tick()$print() - suppressMessages( - verschillen <- as.double(st_area(st_difference(postcodes4 %>% slice(pc), - kaart)) / - st_area(postcodes4 %>% slice(pc))) - ) - if (any(verschillen < 1)) { - newvar[pc] <- as.character(kaart[, 1, drop = TRUE])[verschillen == min(verschillen)] - } else { - kaart_ind <- as.double(suppressMessages(st_intersects(postcodes4 %>% slice(pc), kaart)))[1] - newvar[pc] <- as.character(kaart[, 1, drop = TRUE])[kaart_ind] - } - } - newdf <- data.frame(postcode = as.double(postcodes4$postcode), - newvar = as.character(newvar), stringsAsFactors = FALSE) - postcodes <- postcodes %>% - left_join(newdf, by = "postcode") - colnames(postcodes)[colnames(postcodes) == "newvar"] <- relevante_lagen[i] - } - - object_naam <- case_when(relevante_lagen[i] == "gemeente" ~ "gemeenten", - relevante_lagen[i] == "nuts3" ~ "nuts3regios", - TRUE ~ paste0(relevante_lagen[i], "s")) - - # kaart opslaan in Global Environment - assign(x = object_naam, - value = kaart, - envir = globalenv()) - # kaart opslaan op schijf - saveRDS(object = kaart, - file = paste0(object_naam, ".rds"), - version = 2, - compress = "xz") -} - -# uit PC4-kaart van CBS ook nog wat kolommen halen, en die hoeven niet in dat kaartobject -postcodes <- postcodes %>% - left_join(postcodes_plaats_gemeente, by = "postcode") %>% - select(postcode, matches("inwoner"), "plaats", "gemeente", "provincie", everything()) %>% - left_join(postcodes4 %>% - as.data.frame(stringsAsFactors = FALSE) %>% - select(-area_km2, -geometry), - by = "postcode") - -# alles van PC2 en PC3 toevoegen -postcodes <- postcodes %>% - bind_rows(postcodes %>% - group_by(postcode = clean_numeric(substr(postcode, 1, 2))) %>% - summarise_all(function(x, ...) { - if (is.numeric(x) & any(x > 20)) { - # inwoners en aantal huishoudens - sum(x, na.rm = TRUE) - } else if (is.numeric(x)) { - # gemiddelde grootte van huishoudens - mean(x, na.rm = TRUE) - } else { - x[1] - } - })) %>% - bind_rows(postcodes %>% - group_by(postcode = clean_numeric(substr(postcode, 1, 3))) %>% - summarise_all(function(x, ...) { - if (is.numeric(x) & any(x > 20)) { - # inwoners en aantal huishoudens - sum(x, na.rm = TRUE) - } else if (is.numeric(x)) { - # gemiddelde grootte van huishoudens - mean(x, na.rm = TRUE) - } else { - x[1] - } - })) %>% - arrange(postcode) - +library(dplyr) +library(sf) +linbrary(cleaner) + +# De %like% functie, zodat je kunt gebruiken: "Test" %like% "te" (= TRUE) +`%like%` <- function (x, pattern) grepl(pattern = tolower(pattern), x = tolower(x), fixed = FALSE) + + +# Bronnen ----------------------------------------------------------------- + +downloadmap <- "/Users/msberends/Downloads/" + +# download inwoners (gescheiden op geslacht) per postcode hier als 'CSV volgens tabelindeling': +# https://opendata.cbs.nl/#/CBS/nl/dataset/83503NED/table?dl=42FC6 +# verwijder dan de eerste rijen en de laatste rij ("Bron: CBS") +postcodes_bestand <- paste0(downloadmap, "Bevolking__geslacht__migratieachtergrond__viercijferige_postcode__1_januari_15102020_095747.csv") +# download inwoners per 5 jaar leeftijd en geslacht hier voor het huidige jaar: +# https://opendata.cbs.nl/#/CBS/nl/dataset/83502NED/table?dl=42FE0 +inwoners_bestand <- paste0(downloadmap, "Bevolking__leeftijd__postcode_15102020_102334.csv") +# gebiedsindelingen hier downloaden: +# https://www.cbs.nl/nl-nl/dossier/nederland-regionaal/geografische-data/cbs-gebiedsindelingen +gebiedsindelingen_bestand <- paste0(downloadmap, "cbsgebiedsindelingen_2021_v1.gpkg") +# download postcodes 4 onder 'Downloads' ('Numeriek deel van de postcode (PC4)') hier: +# https://www.cbs.nl/nl-nl/dossier/nederland-regionaal/geografische-data/gegevens-per-postcode +postcodes4_bestand <- paste0(downloadmap, "2019-CBS_PC4_2018_v1/CBS_PC4_2018_v1.shp") + + +# Helpfuncties ------------------------------------------------------------ + +kaart_fixen <- function(kaart) { + # CRS opnieuw instellen, bij nieuwe GDAL geeft dit anders problemen: + st_crs(kaart) <- clean_numeric(st_crs(kaart)$input) + if (is.na(st_crs(kaart))) { + # 'reserve' CRS + st_crs(kaart) <- 28992 + } + if (!"geometry" %in% colnames(kaart)) { + geom_naam <- c("shape", "geom") + geom_naam <- colnames(kaart)[colnames(kaart) %in% geom_naam][1] + geometry <- kaart[, colnames(kaart)[tolower(colnames(kaart)) == geom_naam], drop = TRUE] + kaart <- st_drop_geometry(kaart) + kaart <- st_set_geometry(kaart, geometry) + } + # geometrie als breedte- en lengtegraad, CRS transformeren van 28992 naar 4326 + kaart <- st_transform(kaart, crs = 4326) + + # alle ongeldige vormen geldig maken + if (any(!st_is_valid(kaart$geometry))) { + kaart$geometry <- st_make_valid(kaart$geometry) + } + + # CRS opnieuw instellen, bij nieuwe GDAL geeft dit anders problemen: + st_crs(kaart) <- clean_numeric(st_crs(kaart)$input) + if (is.na(st_crs(kaart))) { + # 'reserve' CRS + st_crs(kaart) <- 28992 + } + # geometrie als breedte- en lengtegraad, CRS transformeren van 28992 naar 4326 + kaart <- st_transform(kaart, crs = 4326) + # oppervlakte toevoegen + kaart$area_km2 <- as.double(st_area(st_geometry(kaart)) / 1000 ^ 2) + kaart +} + +lagen_beschikbaar <- sort(st_layers(gebiedsindelingen_bestand)$name) +downloaden_transformeren <- function(laag) { + zoeklaag <- sort(lagen_beschikbaar[lagen_beschikbaar %like% laag & + lagen_beschikbaar %like% "gegeneraliseerd" & + !lagen_beschikbaar %like% "niet"]) + if (length(zoeklaag) == 0) { + stop("Geen laag gevonden") + } + zoeklaag <- zoeklaag[length(zoeklaag)] + message("Laag '", zoeklaag, "' gevonden", appendLF = FALSE) + + kaart <- st_read(gebiedsindelingen_bestand, layer = zoeklaag, quiet = TRUE) + kaart <- kaart_fixen(kaart) + + message(", met ", nrow(kaart), " geometrieën") + colnames(kaart)[colnames(kaart) %like% "naam"] <- laag + # geen factor + kaart[, laag] <- as.character(kaart[, laag, drop = TRUE]) + + # alleen het type, oppervlakte en de geometrie behouden + kaart <- kaart[, c(laag, "area_km2", "geometry"), drop = FALSE] + kaart +} + + +# Inwoners per postcode/leeftijd/geslacht --------------------------------- + +inwoners_per_postcode_leeftijd <- read_csv2(inwoners_bestand) +inwoners_per_postcode_leeftijd <- inwoners_per_postcode_leeftijd %>% + filter(!leeftijd %like% "totaal", postcode %like% "[0-9]") %>% + mutate(postcode = as.double(postcode), + geslacht = case_when(geslacht %like% "totaal" ~ "inwoners", + geslacht %like% "man" ~ "inwoners_man", + geslacht %like% "vrouw" ~ "inwoners_vrouw"), + # bij CBS vinden ze 0-5 en 5-10 handig. Waarin zit iemand dan die 5 is?! We maken er 0-4 en 5-9 van: + leeftijd_min = as.numeric(gsub("([0-9]+).*", "\\1", leeftijd)), + leeftijd_max = as.numeric(gsub("([0-9]+)[^0-9]+([0-9]+)[^0-9]+", "\\2", leeftijd)), + leeftijd_nieuw = paste0(leeftijd_min, "-", leeftijd_max - 1), + leeftijd_nieuw = gsub("95-NA", "95+", leeftijd_nieuw), + leeftijd = factor(leeftijd_nieuw, levels = levels(age_groups(0, 5 * c(1:19))), ordered = TRUE)) %>% + select(postcode, geslacht, leeftijd, inwoners) %>% + pivot_wider(names_from = geslacht, values_from = inwoners) +# alle PC2 en PC3 toevoegen +inwoners_per_postcode_leeftijd <- inwoners_per_postcode_leeftijd %>% + bind_rows(inwoners_per_postcode_leeftijd %>% + group_by(postcode = clean_numeric(substr(postcode, 1, 2)), leeftijd) %>% + summarise_all(sum)) %>% + bind_rows(inwoners_per_postcode_leeftijd %>% + group_by(postcode = clean_numeric(substr(postcode, 1, 3)), leeftijd) %>% + summarise_all(sum)) %>% + arrange(postcode, leeftijd) + + +# Postcodes (wordt later alle referentiedata aan toegevoegd) -------------- + +# `postcodes` is hier de vorige versie die we als `postcodes` gebruikten, deze wordt vernieuwd +postcodes_plaats_gemeente <- postcodes %>% + filter(postcode > 999) %>% # alleen PC4 houden, wordt later weer aangevuld met PC2 en PC3 + select(postcode, plaats, gemeente) + +postcodes <- read_csv2(postcodes_bestand) +colnames(postcodes) <- gsub("(man|vrouw)n?en", "\\1", colnames(postcodes)) + + +# Postcode-4 kaart -------------------------------------------------------- + +# we gebruiken deze postcodekaart om te bepalen welke postcodes in welk gebied liggen met sf::st_intersects(). +# dus de geometrie van postcode 9251 valt in het snijvlak van de de geometrie van de gemeente Tytsjerksteradiel +# en dus is Tytsjerksteradiel de gemeente van postcode 9251 (en zo verder voor NUTS-3, GGD-regio, ...) +postcodes4 <- st_read(postcodes4_bestand) +postcodes4 <- kaart_fixen(postcodes4) +# alleen relevante kolommen houden +postcodes4 <- postcodes4 %>% + transmute(postcode = as.double(as.character(PC4)), + huishoudens = ifelse(AANTAL_HH < 0, NA_real_, AANTAL_HH), + huishouden_grootte = ifelse(GEM_HH_GR < 0, NA_real_, GEM_HH_GR), + area_km2 = as.double(st_area(geometry) / 1000 ^ 2), + geometry) + + +# Referentiewaarden aan `postcodes` toevoegen en kaarten opslaan ---------- + +relevante_lagen <- c("gemeente", + "provincie", + "nuts3", + "ggdregio", + "jeugdregio", + "veiligheidsregio") +for (i in 3:length(relevante_lagen)) { + message(">> zoeken naar ", relevante_lagen[i]) + kaart <- downloaden_transformeren(relevante_lagen[i]) + + if (!relevante_lagen[i] %in% c("plaats", "gemeente")) { + # referentiedata toevoegen aan 'postcodes' + p <- dplyr::progress_estimated(length(postcodes4$geometry)) + newvar <- character(length = nrow(postcodes4)) + for (pc in 1:nrow(postcodes4)) { + p$tick()$print() + suppressMessages( + verschillen <- as.double(st_area(st_difference(postcodes4 %>% slice(pc), + kaart)) / + st_area(postcodes4 %>% slice(pc))) + ) + if (any(verschillen < 1)) { + newvar[pc] <- as.character(kaart[, 1, drop = TRUE])[verschillen == min(verschillen)] + } else { + kaart_ind <- as.double(suppressMessages(st_intersects(postcodes4 %>% slice(pc), kaart)))[1] + newvar[pc] <- as.character(kaart[, 1, drop = TRUE])[kaart_ind] + } + } + newdf <- data.frame(postcode = as.double(postcodes4$postcode), + newvar = as.character(newvar), stringsAsFactors = FALSE) + postcodes <- postcodes %>% + left_join(newdf, by = "postcode") + colnames(postcodes)[colnames(postcodes) == "newvar"] <- relevante_lagen[i] + } + + object_naam <- case_when(relevante_lagen[i] == "gemeente" ~ "gemeenten", + relevante_lagen[i] == "nuts3" ~ "nuts3regios", + TRUE ~ paste0(relevante_lagen[i], "s")) + + # kaart opslaan in Global Environment + assign(x = object_naam, + value = kaart, + envir = globalenv()) + # kaart opslaan op schijf + saveRDS(object = kaart, + file = paste0(object_naam, ".rds"), + version = 2, + compress = "xz") +} + +# uit PC4-kaart van CBS ook nog wat kolommen halen, en die hoeven niet in dat kaartobject +postcodes <- postcodes %>% + left_join(postcodes_plaats_gemeente, by = "postcode") %>% + select(postcode, matches("inwoner"), "plaats", "gemeente", "provincie", everything()) %>% + left_join(postcodes4 %>% + as.data.frame(stringsAsFactors = FALSE) %>% + select(-area_km2, -geometry), + by = "postcode") + +# alles van PC2 en PC3 toevoegen +postcodes <- postcodes %>% + bind_rows(postcodes %>% + group_by(postcode = clean_numeric(substr(postcode, 1, 2))) %>% + summarise_all(function(x, ...) { + if (is.numeric(x) & any(x > 20)) { + # inwoners en aantal huishoudens + sum(x, na.rm = TRUE) + } else if (is.numeric(x)) { + # gemiddelde grootte van huishoudens + mean(x, na.rm = TRUE) + } else { + x[1] + } + })) %>% + bind_rows(postcodes %>% + group_by(postcode = clean_numeric(substr(postcode, 1, 3))) %>% + summarise_all(function(x, ...) { + if (is.numeric(x) & any(x > 20)) { + # inwoners en aantal huishoudens + sum(x, na.rm = TRUE) + } else if (is.numeric(x)) { + # gemiddelde grootte van huishoudens + mean(x, na.rm = TRUE) + } else { + x[1] + } + })) %>% + arrange(postcode) +