[R-sig-Geo] gdal color tables

Agustin Lobo alobolistas at gmail.com
Fri Jan 4 11:50:08 CET 2013 

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.

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

More information about the R-sig-Geo mailing list