[R-sig-Geo] gdal color tables
Agustin Lobo
alobolistas at gmail.com
Fri Jan 4 13:50:50 CET 2013
For the case of the quantitative 8bit images, it's actually easier, just tune
the palette to zlim and then apply Barry's functions. I've made tuneZpalette:
tuneZpalette <- function(palette,ncol,zlim){
zlim <- zlim+c(-1,1)
micolv <- classIntervals(zlim[1]:zlim[2], n=ncol,style="equal")
micolv <- findColours(micolv,palette(ncol))
micolv <- c(rep(micolv[1],zlim[1]),micolv,rep(micolv[length(micolv)],255-zlim[2]))
micolv
}
Then the following R plot:
require(colorRamps)
require(raster)
rin <- raster("test2R.tif")
plot(rin,zlim=c(99,176) ,col=matlab.like(64),main=" ",axes=FALSE)
and the tif file created by the following code:
micolv <- tuneZpalette(palette=matlab.like,ncol=64,zlim=c(100,175))
writePaletteVRT("test2R.vrt", rin, micolv)
system("python addPalette.py test2R.vrt test2R.tif")
have the same colors, distributed within the interval 100,175
(and with 0-99 as deepest blue and 176-255 as deepest red)
Agus
On Fri, Jan 4, 2013 at 11:50 AM, Agustin Lobo <alobolistas at gmail.com> 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.
>
> 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