[R-sig-Geo] gdal color tables

Thu Jan 3 20:20:43 CET 2013

On Thu, 3 Jan 2013, Agustin Lobo wrote:

> Happy 2013 to everyone!
> Thanks Barry, very useful page.
> Nevertheless, I think the equivalent to zlim in plot() should be
> included also, because most of the time
> we use standard palettes that must be applied to the particular range
> of interest of the raster(s) (note
> real raster data often has outliers), i.e.:
> plot(r,col=matlab.like(64),zlim=c(20,50))
>

This is exactly why I referred to classInt, which returns the indices 
that you must have to attach a color table to a file. With a color table, 
the indices (zero-based) point to rows in the color table, so using 
findIntervals() or similar as used in classInt is the obvious route to 
provide standard palettes.

Roger

> The same would be useful for makePalette(). Currentely, my raster is
> almost all black if use
> writePaletteVRT("test.vrt", rin, matlab.like(64))
> system("python addPalette.py test.vrt test.tif")
>
> By now i'll use plot() to write a jpeg with the appropriate zlim and
> use attachpct.py to apply
> the palette to all rasters, or extract the tuned palette from the
> output of plot() and use your functions
> as they are now. By some reason I'm having problems with raster on my
> home computer and
> cannot use plot(), but will try
> tomorrow from the office. I'll let you know.
>
> Agus
>
> On Thu, Dec 27, 2012 at 11:35 PM, Barry Rowlingson
> <b.rowlingson at lancaster.ac.uk> 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.
>>
>> 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
>>
>> _______________________________________________
>> 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