Here's my two cents,

The way I see it, you have three options to convert a greyscale image into a colour-mapped version.

1. You can just import your greyscale data into ArcGIS/QGIS and tell it how to render that using the colourmap you want. In QGIS, this option is called "Singleband Pseudocolor".
2. You can create your own colourmap dictionary from matplotlib (or another library) and add that to the tif with `write_colormap()`
3. You can modify the greyscale image to just BE the colours.

``` python
import rasterio
import numpy as np
import matplotlib.cm

cmap = matplotlib.cm.get_cmap('viridis')
data = np.random.uniform(0, 256, (100, 100)).astype(np.uint8)

# Option 1: Just use greyscale, and then interpret it in GIS software
with rasterio.open('greyscale.tif', 'w', width=100, height=100, count=1, dtype=np.uint8) as dst:
    dst.write(data, 1)

# Option 2: Explicitly add a cmap obtained from somewhere else
cmap_dict = {n: np.array(cmap(n))*256 for n in range(256)}
with rasterio.open('mapped.tif', 'w', width=100, height=100, count=1, dtype=np.uint8) as dst:
    dst.write(data, 1)
    dst.write_colormap(1, cmap_dict)

# Option 3: Save the greyscale as a coloured image (not recommended)
colour = (cmap(data)*256).astype(np.uint8) # shaped [100, 100, {RGBA}]
with rasterio.open('coloured.tif', 'w', width=100, height=100, count=3, dtype=np.uint8) as dst:
    for x in range(3):
        dst.write(colour[..., x], x+1)
```

I think you are asking about option 2. The answer, then, is that you can just use matplotlib's colourmap library to generate your colormap for you.

Most GIS software allow you to do Option 1 very easily, anyway. So if it's a one-off, it's probably easier to create/select your colormap using a GUI.

Good luck.