I found the solution on this awesome website
addColorTable <- function(inRstName, outRstName, rat.df){
library(rgdal)
r<- readGDAL(inRstName)
rat.df$color<- as.character(rat.df$color)
rat.df$attribute<- as.character(rat.df$attribute)
outRst <- writeGDAL(r, outRstName, type="Byte",
colorTable=list(rat.df$color),
catNames=list(rat.df$attribute), mvFlag=11L)
return(raster(outRst))
}
library(rgdal)
library(raster)
# create dummy data set
r <- raster(nrow=10, ncol=10)
r[] <- 0
r[51:100] <- 1
r[3:6, 1:5] <- 2
r[1, 1] <- 3
writeRaster(r,'dummy_raster.tif',overwrite=T)
# This defines the values, the color and the attribute
valT <- c(0,1,2,3)
colT <- c("#FF0000", "#FF9900" ,"#99FF00","#0000FF")
attT <- c('Forest','Water body','City','Cropland')
rat.df <- data.frame(value=valT,color=colT,attribute=attT)
# apply the magic function
rnew <- addColorTable('dummy_raster.tif', 'dummy_raster_with_symbology.tif', rat.df)
# plot the results and tadaaaa
spplot(rnew, col.regions=colT)