[R-sig-Geo] gdal color tables

[Barry Rowlingson]

Here's a solution:

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
}


Then:

 write a raster tiff

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


Then use makePalette to create the colortable lines:

cat(makePalette(iom at legend@colortable)) # replace
iom at legend@colortable with your colour table vector

That spits out a bunch of <Entry> lines. Stick them in your .vrt file thus:

<VRTDataset rasterXSize="413" rasterYSize="397">
  <VRTRasterBand dataType="Byte" band="1">
    <ColorInterp>Palette</ColorInterp>
<SimpleSource>
<SourceFilename relativeToVRT="1">test.tif</SourceFilename>
</SimpleSource>
    <ColorTable>
<Entry c1="0" c2="0" c3="0" c4="255"/>
<Entry c1="230" c2="0" c3="77" c4="255"/>
<Entry c1="255" c2="0" c3="0" c4="255"/>
<Entry c1="204" c2="77" c3="242" c4="255"/>
<Entry c1="204" c2="0" c3="0" c4="255"/>
[etc]
   </ColorTable>
  </VRTRasterBand>
</VRTDataset>

Note you have to fill in the rasterXSize and rasterYSize, and the
source filename. Also you'll need to get the geotransform of the
original and stick that in otherwise your raster is geolocated at
(1:Nrows, 1:Ncolumns).

If I do all that, I get a raster created by R that I can read into
QGis and is coloured according to my colour scheme, as long as I open
*test.vrt* in QGis and *not* test.tif.

I'm sorry this isn't (yet) a one-shot solution, and its a bit of a
construction set. Like I said, two minute job, twenty minute job to do
properly... I've spent ten minutes on it :)

Barry




On Sun, Dec 23, 2012 at 9:25 PM, Rowlingson, Barry
<b.rowlingson at lancaster.ac.uk> wrote:
> On Sun, Dec 23, 2012 at 8:23 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>
>> There is no support in raster or sp objects for symbology, and unless
>> someone extends the classes to accommodate symbology on a dynamic
>> by-attribute basis, it isn't going to happen.
>>
>
>  I'm not sure that's true - I can read in a raster from a geoTIFF and
> get a colour table from it, then plot will show it with the right
> colours - here's a Land Use tiff:
>
>> lu
> class       : RasterLayer
> dimensions  : 397, 413, 163961  (nrow, ncol, ncell)
> resolution  : 100, 100  (x, y)
> extent      : 3356200, 3397500, 3533200, 3572900  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
> +ellps=GRS80 +units=m +no_defs
> data source : /data/rowlings/MapLibrary/Europe/UK/IsleOfMan/iomcorine.tiff
> names       : iomcorine
> values      : 0, 255  (min, max)
>
>> lu at legend@colortable[1:10]
>  [1] "#000000" "#E6004D" "#FF0000" "#CC4DF2" "#CC0000" "#E6CCCC" "#E6CCE6"
>  [8] "#A600CC" "#A64D00" "#FF4DFF"
>
>
>  Support for this is minimal in package:raster, and I think what Agus
> is trying to do (correct me if wrong) is create a raster with a colour
> table such as this from R. If I just use writeRaster I get a numeric
> raster with no palette.
>
> I think Robert has talked about better support for paletted rasters in the past.
>
> Barry
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo