Writing Paletted Rasters

The raster package doesn't have much support for paletted rasters. You can read them in, and the colour palette is stored in a slot and used for plotting. But how can you write a raster with a colour palette?

First lets create a small raster with some integer values:

library(raster)

## Loading required package: sp

## raster 2.0-12 (1-September-2012)

r = raster(matrix(c(1, 1, 2, 0, 2, NA, 3, 4, 4, 5, 0, 5), 3, 4), xmn = -2, xmx = 2, 
    ymn = 54, ymx = 55)
projection(r) = "+init=epsg:4326"
plot(r)

This plot shows we have a raster that is being treated as continuous values. We want to add a colour table. First, we write our raster to a geoTiff using single byte values:

writeRaster(r, "r.tif", datatype = "INT1U", overwrite = TRUE)

## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.1, released 2012/05/15 Path to GDAL shared
## files: /usr/share/gdal/1.9 Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September
## 2009, [PJ_VERSION: 470] Path to PROJ.4 shared files: (autodetected)

## class       : RasterLayer 
## dimensions  : 3, 4, 12  (nrow, ncol, ncell)
## resolution  : 1, 0.3333  (x, y)
## extent      : -2, 2, 54, 55  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/rowlings/Work/R/Raster/Paletted/r.tif 
## names       : r 
## values      : 0, 5  (min, max)

We need a few support routines:

makePalette <- function(colourvector) {
    cmat = cbind(t(col2rgb(colourvector)), 255)
    res = apply(cmat, 1, function(x) {
        sprintf("<Entry c1=\"%s\" c2=\"%s\" c3=\"%s\" c4=\"%s\"/>", x[1], x[2], 
            x[3], x[4])
    })
    res = paste(res, collapse = "\n")
    res
}

makePaletteVRT <- function(raster, colourvector) {
    s = sprintf("<VRTDataset rasterXSize=\"%s\" rasterYSize=\"%s\">\n<VRTRasterBand dataType=\"Byte\" band=\"1\">\n<ColorInterp>Palette</ColorInterp>\n<ColorTable>\n", 
        ncol(raster), nrow(raster))
    p = makePalette(colourvector)
    s = paste0(s, p, "\n</ColorTable>\n</VRTRasterBand>\n</VRTDataset>\n")
    s
}

writePaletteVRT <- function(out, raster, colourvector) {
    s = makePaletteVRT(raster, colourvector)
    cat(s, file = out)
}

These functions create a small file containing the colour palette. GDAL colour palettes start from zero, so we need six colours to cover 0 to 5 in our raster. Let's write that out:

writePaletteVRT("test.vrt", r, rainbow(6))

Now we need to add that palette to the raster file. Currently I don't think rgdal exposes enough of GDAL to do it, so we need a few lines of python. Cut this out and save it as addPalette.py in your working directory or anywhere else in your $PATH - you'll need a python with gdal libraries installed.

#! /usr/bin/env python
#
# take a raster and add a palette
#

#  1. write an integer raster from R like this: writeRaster(r,"rint.tiff",datatype="INT1U")
#
#  2. write a palette file from R like this:  writePaletteVRT("test.vrt",r,rainbow(20))
#
#  3. run this script like this: 
#      addPalette test.vrt rint.tiff
#     and it will modify rint.tiff, adding the palette.

try:
    from osgeo import gdal
except ImportError:
    import gdal

import sys

def Usage():
    print('Usage: addPalette.py palette_file raster_file')
    sys.exit(1)

gdal.AllRegister()
argv = gdal.GeneralCmdLineProcessor( sys.argv )
if argv is None:
    Usage()

if len(argv)!=3:
    Usage()

pal_file = argv[1]
raster_file = argv[2]

src_ds = gdal.Open(raster_file)

pal_ds = gdal.Open(pal_file)
ct = pal_ds.GetRasterBand(1).GetRasterColorTable().Clone()
src_ds.GetRasterBand(1).SetRasterColorTable( ct )

gtiff_driver = gdal.GetDriverByName( 'GTiff' )

gtiff_driver.CreateCopy(raster_file,src_ds)

We can run that from R:

system("./addPalette.py test.vrt r.tif")

and if we read that file back in and plot it, we should see the rainbow palette:

r2 = raster("r.tif")
print(r2)

## class       : RasterLayer 
## dimensions  : 3, 4, 12  (nrow, ncol, ncell)
## resolution  : 1, 0.3333  (x, y)
## extent      : -2, 2, 54, 55  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/rowlings/Work/R/Raster/Paletted/r.tif 
## names       : r 
## values      : 0, 5  (min, max)

plot(r2)

Note that there are some problems with the raster package here - the plot doesn't have the correct aspect ratio even though it has the same extent and projection as the source, and missing data shows up the same as colour zero. If you load this raster into Quantum GIS however, it shows up correctly, with the right aspect ratio and with the missing data cell as transparent.

Note that the colour table is stored in the following slot of the raster:

head(r2@legend@colortable)

## [1] "#FF0000" "#FFFF00" "#00FF00" "#00FFFF" "#0000FF" "#FF00FF"

Obviously these three operations - write integer raster, create colourmap file, add colourmap file to raster - could be all incorporated into a single function, but I'm hoping for better paletted raster functionality in the raster package soon.

You are free to copy and use the code in this document for any purpose, no guarantees of efficacy or safety are given.