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!

@ tospig: thank you for your reply! I had never considered this restriction before! I am open to any number of colors! E.g. Let's say 6? Red, Blue, Green, Yellow, Purple, Orange? But I am open to any thing! Thank you so much! — stats_noob, Apr 20 '23 at 22:19
This is not an "R" question. I recommend you do a search for pages discussing simple (so to speak) implementations of 4-color map assignments. — Carl Witthoft, Apr 20 '23 at 23:46
igraph has `greedy_vertex_coloring()`. It will not in general be able to do it with only four colours. You will of course need to construct the adjacency graph of the polygons. — Szabolcs, Apr 21 '23 at 10:19
@ Szabolcs: thank you so much for your reply! I can't find this function in R? — stats_noob, Apr 21 '23 at 15:08
> library(igraph) > g <- make_graph("petersen") > greedy_vertex_coloring(g) Error in greedy_vertex_coloring(g) : could not find function "greedy_vertex_coloring" — stats_noob, Apr 21 '23 at 15:08
I see you added a bounty. Did you try using `greedy_vertex_coloring()`? Did you encounter any specific issue that you couldn't solve? — Szabolcs, Apr 23 '23 at 14:53
@ Szabolcs: thank you for your reply! I tried this function over here: https://stackoverflow.com/questions/76075267/r-assigning-colors-to-graph-nodes — stats_noob, Apr 23 '23 at 15:06

dww · Accepted Answer · 2023-04-29T18:52:00.557

2

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

edited Apr 29 '23 at 18:52

answered Apr 23 '23 at 17:17

dww

30,425
5
68
111

@ dww: thank you so much for your answer! Is the resulting map an interactive leaflet map? thanks! – stats_noob Apr 23 '23 at 17:20
sure - I plotted it using the same leaflet code you had in the Q. I can add that to the answer if you like, but its just a copy and paste from your example – dww Apr 23 '23 at 17:23
1

Regarding the claim that this "will give you the minimal (four color) solution": the MapColoring package claims to use the DSatur algorithm, which does not guarantee an optimal solution, not even for planar graphs. I know that the package also claims to use a minimal number of colours, but if it relies on DSatur, this is simply not true. – Szabolcs Apr 24 '23 at 05:55
BTW igraph 0.10 contains an accurate implementation of DSatur (MapColoring's implementation has some minor deviations from the paper). It will be available in `greedy_vertex_coloring()` once R/igraph is updated to use 0.10, hopefully this summer. – Szabolcs Apr 24 '23 at 06:09
1

Thanks for your comments @Szabolcs. I updated the answer to be more accurate. – dww Apr 24 '23 at 09:16
Thanks for adding the clarification. One more comment I have here is that the current text compares DSatur to "the greedy algorithm". In fact, "greedy" means that we colour vertices one by one in some order, never modifying the already assigned colours. There are many way to choose this order, thus there are many different greedy algorithms. DSatur itself is greedy too! It's just one of the many greedy algorithms. This is why DSatur is just one of the available methods in the `igraph_vertex_coloring_greedy()` function of C/igraph 0.10 (coming soon to R/igraph). – Szabolcs Apr 24 '23 at 10:07
Thus it doesn't make sense to compare DSatur to "the greedy algorithm". In R/igraph 1.4, `greedy_vertex_coloring()` implements a single ordering heuristic: the next vertex we colour is the one with the most already coloured neighbours. DSatur chooses the one with the most neighbouring colours. Both have some extra subtleties on how to break ties. There are several other possible ordering heuristics as well, e.g. NetworkX implements 7 different ones at the moment: https://networkx.org/documentation/stable/reference/algorithms/generated/networkx.algorithms.coloring.greedy_color.html – Szabolcs Apr 24 '23 at 10:12
@ dww: I ran your code (after loading tmaptools library)... but how do I view the map? – stats_noob Apr 29 '23 at 03:18
You can just use your original code to view the map. I added that to the answer. – dww Apr 29 '23 at 18:47

Henrik · Answer 2 · 2023-04-23T17:54:44.247

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

edited Apr 23 '23 at 17:54

answered Apr 23 '23 at 17:26

Henrik

65,555
14
143
159

I tried tmaptools::map_coloring already but the algorithm failed to converge on 4 colors (it used a fifth color for one county). Can you check this works as expected, as I think it depends on the same function. – dww Apr 23 '23 at 17:27
@dww The code above works on OP's `nc`. I haven't hardcoded the number of colors. Cheers – Henrik Apr 23 '23 at 17:33
I tried it out, and as you can see in [this image](https://imgur.com/a/RL5tH5z) there is an extra 5th color for one county – dww Apr 23 '23 at 17:42
1

@dww Indeed, but also note [OP's comment](https://stackoverflow.com/questions/76068597/color-polygons-in-a-map-so-that-adjacent-polygons-have-different-colors/76086263?noredirect=1#comment134155817_76068597): _I am open to any number of colors! E.g. Let's say 6? Red, Blue, Green, Yellow, Purple, Orange? But I am open to any thing!_ . I don't claim that my code results in four colours ;) It's just _a_ (convenient) way to solve the problem. Cheers – Henrik Apr 23 '23 at 17:46
@ Henrik: Thank you so much for your answer! I really like this feature in tmaptools that it allows you to click on any polygon and see the value of all columns for the row corresponding to that polygon from the shapefile. very nice! – stats_noob Apr 29 '23 at 15:49
@stats_noob That's not a feature of `tmaptool`, but rather the `tmap` defaults when rendering a `leaflet` plot; as I wrote in my answer `tmap_mode(mode = "view")` generates an interactive `leaflet` widget. Then check the `popup.vars` argument of `tm_fill` (which is called from `tm_polygons`): "names of data variables that are shown in the popups in `"view"` mode. [...] If they are not specified, all variables are shown." – Henrik Apr 29 '23 at 20:13

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

2 Answers2

Linked