Get colortable from three-band RGB raster (brick) with R

Question

I have a three band RGB Raster Geotiff. I will provide the sample data below, but I'm guessing my problem should be solvable generically. I would like to turn this three band raster into a single band raster with a color table.

The original three bands can each range from 1 to 255, resulting in over 16M (256^3) theoretically possible combinations of Red, Green and Blue. Of these 16M+ Combinations, little more than 11'000 combinations are actually used. Be that as it may, the resulting color table allows only 256 values in total.

Is there an algorithm to reduce the 16M+ possible RGB-values to an appropriate set of 256 values?

I'm guessing that I'd first have to decide on a colour palette with 256 colours and then somehow assign each RGB Value to one of these colours. Before I reinvent the wheel, I thought I'd ask here whether this problem already has been solved. I'd like to use R, the package raster and any other necessary package for this.

Note: Interestingly, the producers of said map also provide a single band geotiff with just 13 (!!) different colours in the colortable, but it doesn't look much different from the three band Raster (see below).

This is the three band RGB raster:

This is the single band raster with a color table:

See the sample data on the swisstopo homepage. In case the Link breaks: https://shop.swisstopo.admin.ch > "Digital national maps > Swiss Map Raster > Swiss Raster Map 50 > Sample Data > Swiss Map Raster 50, Sample data. The three band RGB Raster is in the folder SMR_50KREL_254dpi_uncompressed/SMR50_LV95_KREL_10L_Mosaic.tif. The single band Raster with a colortable is here: SMR_50KOMB_508dpi_LZW/SMR50_LV95_KOMB_Mosaic.tif

Answer 1

There's a function in rgdal for this SGDF2PCT, so here I coerce to SpatialGridDataFrame, build the colour table and rebuild the raster. Note that indexing in raster is assuming [0, 255].

Control the number of colours with ncolors argument.

library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
pct <- rgdal::SGDF2PCT(as(b, "SpatialGridDataFrame"))

r <- setValues(raster(b), pct$idx-1)
colortable(r) <- pct$ct
plot(r)

r now has one layer and the colour table in r@legend@colortable (use colortable(r)).