Color polygons in a map so that adjacent polygons have different colors

Question

I made the following map:

library(sf)  
library(leaflet)
library(leafgl)
library(colourvalues)
library(leaflet.extras)


nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) %>% 
  st_transform(st_crs(4326)) %>% 
  st_cast('POLYGON')

leaflet(data = nc) %>% addPolygons( stroke = FALSE) %>% addTiles(group = "OSM") %>%  addProviderTiles(provider = providers$OpenStreetMap) %>% addPolygons(data = nc, weight=1, popup = ~NAME,
                label = ~NAME, group = "name", col = 'blue') %>% 
    addSearchFeatures(targetGroups  = 'name', options = searchFeaturesOptions(zoom=10, openPopup=TRUE))

I wanted to color the polygons different colors so that they are a bit easier to see - I did this by randomly assigning colors to the polygons:

nc$color <- sample(c("red", "blue", "green", "yellow", "purple"), nrow(nc), replace = TRUE)

leaflet(data = nc) %>% 
    addTiles(group = "OSM") %>% 
    addProviderTiles(provider = providers$OpenStreetMap) %>% 
    addPolygons(data = nc, weight=1, popup = ~NAME,
                label = ~NAME, group = "name", fillColor = ~color, fillOpacity = 0.5) %>% 
    addSearchFeatures(targetGroups  = 'name', options = searchFeaturesOptions(zoom=10, openPopup=TRUE))

My Question: Taking inspiration from this famous computer science problem https://en.wikipedia.org/wiki/Four_color_theorem, I would like to randomly color the polygons in such a way that no neighbouring polygons have the same color.

I think that I first need to convert the shapefile/map into a network graph:

library(igraph)
adj <- st_touches(nc, sparse = TRUE)
g <- graph_from_adjacency_matrix(as.matrix(adj))
plot(g)

I am not sure how to continue this problem - currently, I thought of an indirect method where I simply choose many different random colors to decrease the odds of two polygons having the same color, but I am interested in learning about new and creative ways to solve my original problem.

Can someone please show me how to do this?

Thanks!

Answer 1

You can use the MapColoring package. This package uses the DSatur algorithm,which in this case is able to find a minimal (four color) solution to the problem, whereas solutions based on the greedy algorithm do not. In general, DSatur has been shown to derive significantly better vertex colourings than the greedy algorithm.

devtools::install_github("hunzikp/MapColoring")
library(sp)
library(MapColoring)

my.palette = RColorBrewer::brewer.pal(4, 'Set1')
nc$color = my.palette[getColoring(as(nc, 'Spatial'))]

To view the map you can just use the same code you had in your question:

leaflet(data = nc) %>% 
    addTiles(group = "OSM") %>% 
    addProviderTiles(provider = providers$OpenStreetMap) %>% 
    addPolygons(data = nc, weight=1, popup = ~NAME,
                label = ~NAME, group = "name", fillColor = ~color, fillOpacity = 0.5) %>% 
    addSearchFeatures(targetGroups  = 'name', options = searchFeaturesOptions(zoom=10, openPopup=TRUE))

Answer 2

Using tmap:

library(tmap)
tmap_mode(mode = "view") # if you want an interactive leaflet map
tm_shape(nc) + tm_polygons(col = "MAP_COLORS", minimize = TRUE)

In tm_polygons, when col is set to "MAP_COLORS":

polygons will be colored such that adjacent polygons do not get the same color. See the underlying function map_coloring [from tmaptools] for details

From ?map_coloring about the minimize argument:

logical that determines whether algorithm will search for a minimal number of colors

---

Use palette argument in tm_polygons to set preferred colours.

For future reference, related post for non-interactive ggplot map:

Coloring differently adjacent regions on a map with ggplot