Search code examples
rggplot2plotr-sfgeopackage

Fill Subdivisions of Chloropeth Map with multiple Colors


I want to create a plot of Germany, where each of its 16 Bundesländer (Subdivisons) is colored differently, depending on the fulfillment of a handful of criteria. At the bottom you find a representative sketch. The idea is, that you see right away which subdivision meets which criteria.

Each criterium stands for a different color, these are

library(RColorBrewer)
Colors <- brewer.pal(5, "Spectral") 
Colors
[1] "#D7191C" "#FDAE61" "#FFFFBF" "#ABDDA4" "#2B83BA"

Now, I have a list of 16, which contains the color(s) of each Bundesland, depending on which criteria they meet. For example:

Berlin
[1] "#D7191C" "#FDAE61" "#FFFFBF" "#ABDDA4" 
Bavaria
[1] "#D7191C" "#FDAE61"
Hamburg
[1] "#ABDDA4"

As you can see, they all differ in length.

Now, I was wondering, how I could pass these multiple color information to my chloropleth map? I have downloaded the Geopackage for Germany from GADM.

The following gives a plain plot of the map, where you can only see the borders of the Bundesländer (attached)

library(sf)
library(cartography)
plot(st_geometry(Germany), col = NA, border = "white", bg = "#aadaff")

Passing my list of colors as

"col = Colors"

Is, of course, not working, and does not fill the subdivisions with the desired color(s).

Desired Map


Solution

  • Update:

    OP states that they need multiple colors per Bundesländer. Assuming they have a named list or a dataframe with the colors, they can join it with their shapefile. Here, I am creating a set of random colors. I assign the colors to each state (assigning multiple colors to some of them).

    Then I would divide each state that has multiple colors to equal parts.

    Finally, I plot the divided up polygons (without a border) and would plot the original polygons with no color but only a white border (here I chose black to show them more clearly) on top of it to show the original subdivisions/state borders. I am not sure why the background color does not work.

    #### libraries ####
    
    library(sf)
    library(sp)
    library(dismo)
    library(deldir)
    library(dplyr) 
    
    #### downloading the data ####
    
    download.file("https://geodata.ucdavis.edu/gadm/gadm4.1/shp/gadm41_DEU_shp.zip",
                  "gadm41_DEU_shp.zip")
    unzip("gadm41_DEU_shp.zip", exdir= "DEU_adm")
    
    Germany <- st_read(dsn="DEU_adm\\gadm41_DEU_1.shp", quiet = TRUE)
    
    #### pre-processing the data and defining the split funciton ####
    
    ## creating a set of random colors, you should use your desired colors
    set.seed(123)
    cp <- sample(grDevices::colors()[grep('gr(a|e)y|white',
                 grDevices::colors(), invert = T)], 24)
    
    Germany$COLORS <- I(list(cp[1], cp[2:3], cp[3:5], cp[6:9], 
                             cp[10], cp[11], cp[12], cp[13],
                             cp[14], cp[15], cp[16], cp[17],
                             cp[18], cp[19], cp[20], cp[21:24]))
    
    wgs84 <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
    
    ## taken from https://gis.stackexchange.com/a/440698/93948
    ## split polygons 
    split_poly <- function(sf_poly, n_areas) {
      # Create random points
      points_rnd <- st_sample(sf_poly, size = 10000)
      # k-means clustering
      points <- do.call(rbind, st_geometry(points_rnd)) %>%
        as_tibble() %>% setNames(c("lon","lat"))
      k_means <- kmeans(points, centers = n_areas)
      # Create voronoi polygons
      voronoi_polys <- dismo::voronoi(k_means$centers, ext = sf_poly)
      # Clip to sf_poly
      crs(voronoi_polys) <- wgs84 ## hardcoding crs
      voronoi_sf <- st_as_sf(voronoi_polys)
      equal_areas <- st_intersection(voronoi_sf, sf_poly)
      #equal_areas$area <- st_area(equal_areas)
      return(equal_areas)
    }
    
    #### splitting the polygons and plotting ####
    
    ## split polygons for each Bundesländer based on number of colors provided
    pol_areas  <- lapply(seq_len(nrow(Germany)), function(i) 
            if (length(unlist(Germany[i,]$COLORS)) == 1) {Germany[i,]} 
            else {split_poly(Germany[i,], length(unlist(Germany[i,]$COLORS)))}) 
    
    ## combine the splited polygons to one sf object
    bind_rows(pol_areas) %>% 
      rowwise() %>% 
      mutate(Color = ifelse(is.na(id), COLORS, COLORS[id])) -> mapdata_deu 
    
    ## plot
    plot(mapdata_deu[,"geometry"],
         col = mapdata_deu$Color, 
         bg = "#aadaff", border = NA)
    plot(Germany[,"geometry"],
         col = NA, 
         bg = NA, border = "black", add = T)
    

    Created on 2024-04-01 with reprex v2.0.2


    I am not quite sure what are the colors that you want for each state (i.e. Bundesland), but the code below does the job. You can change the Colors variable as desired.

    library(sf)
    library(sp)
    library(raster)
    library(rgdal)
    
    Germany <- getData(country = "Germany", level = 1) 
    
    Colors <- c("#D7191C", "#FDAE61", "#FFFFBF", "#ABDDA4",
                "#D7191C", "#FDAE61", "#ABDDA4", "#D7191C", 
                "#FDAE61", "#FFFFBF", "#ABDDA4", "#2B83BA",
                "#D7191C", "#FDAE61", "#ABDDA4", "#D7191C")
    
    plot(Germany, 
         col = Colors[as.numeric(as.factor(Germany$NAME_1))], 
         border = "white", bg = "#aadaff")
    

    Created on 2024-03-24 with reprex v2.0.2