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)