[R-sig-Geo] gdal color tables

[Roger Bivand]

On Thu, 27 Dec 2012, Barry Rowlingson wrote:

> How to create colour-paletted rasters using R and Python:
>
> http://www.rpubs.com/geospacedman/rasterColourPalettes
>
> It would be great to be able to do this solely from R using either
> raster or rgdal, and to have the couple of bugs with paletted rasters
> in package:raster mentioned in the Rpubs doc fixed. I shall try and
> raise this with Robert H when everyone comes out of hibernation in
> January.

Barry:

As of SVN revision 422 on R-forge, writeGDAL() can write a properly 
configured ColorTable to a Byte GTiff (other drivers untried, probably 
limited to Byte, so #colors 256), by band. It can add CategoryNames by 
band; these and the colors seem to be read by ArcGIS 10.1; the color table 
is read by GRASS 6.4.2, but not the category names. See the examples at 
foot of ?writeGDAL. I haven't so far checked out how one might pass the 
ColorTable through from raster (from the object metadata?), and the 
creation of the table isn't very intuitive (rather like your 
makePalette()) - the look-up from the raster should be zero-based (value 0 
gets the first color, etc.). I'd be interested in feedback with your 
examples.

Best wishes,

Roger


>
> On Sun, Dec 23, 2012 at 9:51 PM, Rowlingson, Barry
> <b.rowlingson at lancaster.ac.uk> wrote:
>> 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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no