gdal color tables

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))

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