cbs

cbs_2021/cbs_data_updaten.R

@@ -0,0 +1,233 @@
+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)
+