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

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!

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))