[R-sig-Geo] gdal color tables

Roger Bivand Roger.Bivand at nhh.no
Fri Jan 4 12:00:26 CET 2013


On Fri, 4 Jan 2013, Agustin Lobo wrote:

> Yes, using classIntervals() is not hard, but note that the VRT table
> does not have a values field. This is not
> a problem for unsigned integer8 tif files (i.e., for landcover raster
> maps) but it is in the case of raster
> files of float values (i.e., many of my aerial thermograms).  I presume
> gdal tif files make a linear interpolation between the min to the max
> value in this case , but I really cannot dig
> into it myself right now.

The GDAL implementation of color tables only works internally for Byte 
valued rasters - try assigning to say a Float32, and GDAL refuses to add 
the table. I don't know how this works in the GDAL/Python implementation, 
but it cannot do more than the GDAL base functionality. There are internal 
color ramp functions in GDAL, but there are no break points or 
quantization rules (that you find in GRASS, which supports color tables 
for floating point rasters).

If you treat classIntervals() as a way of creating quantization rules, you 
can store the floating point values in one file, and the quantized 0-255 
(or fewer) values in a Byte file.

I do see references to quantization in GDAL for specific drivers, such as 
PNG, but not as a general facility - does anyone know whether such a 
facility exists (to add quantized color tables and rules to arbitrary 
raster bands)?

Roger

>
> But I also have a set of thermograms as 8bit images, so for those I
> can hack Barry's functions with classIntervals()
>
> For the thermograms in 16bits integers and floating values, what I can
> do is to modify Barry's functions to write a QGIS qml , very similar
> to the table in the VRT file but
> including a values field.
>
> Agus
>
> On Thu, Jan 3, 2013 at 8:20 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>> 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
>>
>> _______________________________________________
>> 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



More information about the R-sig-Geo mailing list