Dear Roger
Thanks for your detailed reply.
>> I would like to use R from my bash console with the following
>> command:
>>
>> /path/to/R --vanilla --slave --args < /data/myscript.R
>>
>> this script invokes a range of functions and tools, among others to
>> process a bunch of raster images in a loop using GDAL. The thing
>> is that
>> releasing this command in the bash results in the following error:
>
> The possibilities are almost endless. Without having any idea of the
> classes of the arguments to your function, it is hard to know why the
> object passed to writeGDAL() may be corrupt.
>
> Only write a script when you have a good practical understanding of
> all
> the steps you are taking. Your script indicates that this is not
> the case.
> You have not even explained that what you seem to be trying to do
> is to
> take a SpatialGridDataFrame with one column, and to convert it into a
> four-column SpatialGridDataFrame with red, green, blue, and alpha
> columns,
> based on a classification of the input variable given in a numeric
> vector
> breaks.
Yes, the description of what the script intends to do is missing.
Sorry, my fault. Your description of the script is correct. And yes,
it was the first time I did anything with GDAL.
>> ERROR:
>> *************************************
>>
>> Error in GDAL.close(tds.out) :
>> GDAL Error 1: TIFFReadDirectory:/tmp/RtmpBgrSyt/file2ef6a5b2:
>> cannot handle zero scanline size
>> Calls: plot.georef -> GDAL.close -> .Call
>> Execution halted
The weird thing is, that the script works as expected if called
interactively from the R console. It gives the same error, but the
resulting TIFF images are ok. The execution is only halted if R is
called from the command line (with the script name as argument).
>> CODE:
>> *************************************
>>
>> plot.georef=function(grid.sp,graph.file,col.palette,breaks){
>
> grid.sp is what? A SpatialGridDataFrame?
Yes, this SpatialGridDataFrame results from an interpolated
temperature field. writeAsciiGrid(grid.sp,asc.file) results in the
expected text file with a correct georef.
> str(grid.sp)
Formal class 'SpatialGridDataFrame' [package "sp"] with 6 slots
..@ data :'data.frame': 84000 obs. of 2 variables:
.. ..$ temperature_mean.pred: num [1:84000] 5.4 5.4 5.4 5.4 5.5
5.5 5.5 5.5 5.5 5.5 ...
.. ..$ temperature_mean.var : num [1:84000] NA NA NA NA NA NA NA
NA NA NA ...
..@ grid :Formal class 'GridTopology' [package "sp"] with 3
slots
.. .. ..@ cellcentre.offset: Named num [1:2] 655050 278050
.. .. .. ..- attr(*, "names")= chr [1:2] "x_coord" "y_coord"
.. .. ..@ cellsize : Named num [1:2] 100 100
.. .. .. ..- attr(*, "names")= chr [1:2] "x_coord" "y_coord"
.. .. ..@ cells.dim : Named int [1:2] 350 240
.. .. .. ..- attr(*, "names")= chr [1:2] "x_coord" "y_coord"
..@ grid.index : int(0)
..@ coords : num [1:2, 1:2] 655050 689950 278050 301950
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : NULL
.. .. ..$ : chr [1:2] "x_coord" "y_coord"
..@ bbox : num [1:2, 1:2] 655000 278000 690000 302000
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "x_coord" "y_coord"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
.. .. ..@ projargs: chr "+init=world:CH1903"
> graph.file is what? A file name - single element character vector?
> str(graph.file)
chr "/data/graph_output/temp_205_100m_2008031117.tiff"
> col.palette is what? A character vector?
> str(col.palette)
chr [1:26] "#97FF02" "#5BFA05" "#00C632" "#17D9ED" "#28CCEE" ...
> breaks is what? A numeric vector?
> str(breaks)
num [1:27] -100 -18 -16 -14 -12 -10 -8 -6 -4 -2 ...
>> # get desired picture type
>> base.name=strsplit(graph.file,"/")[[1]][length(strsplit
>> (graph.file,"/")[[1]])]
>
> Do you know that base.name is what you think it is? Did you
> consider using
> the basename() function?
Yes, basename() should be used. The above is a relict from an old
script, and does the same as basename().
>> ft=0
>> ft=sum(c(ft,grep(".tiff",base.name)))
>
> See ?grep - "." means something.
>
>> if (ft==1){ file.typ="TIFF" }
>> ft=0
>> ft=sum(c(ft,grep(".jpg",base.name)))
>> if (ft==1){ file.typ="JPG" }
>> ft=0
>> ft=sum(c(ft,grep(".png",base.name)))
>> if (ft==1){ file.typ="PNG" }
This part of the script actually works without problems. I know that
the determination of the file type from the given file name is not
done in a very elaborate way. Any suggestions?
>> # create RGB bands
>> names(grid.sp)=c("band1","band2")
>
> How many columns did grid.sp have when it came into the function,
> do you
> know that there were two? If not, this will not create them.
Yes, it always has two.
>> grid.sp$band3=grid.sp$band2
>> grid.sp$band4=grid.sp$band2
>
> If band2 wasn't created above, this won't create band3 or band4.
Yes, but as band2 is always present, this is not an issue.
>> col.class=as.vector(cut(grid.sp$band1,br = breaks,
>> labels=col.palette))
>
> Curious use of cut:
>
> col.class <- col.palette[findInterval(grid.sp$band1, vec=breaks,
> all.inside=TRUE)]
>
> seems more natural, and avoids creating a factor (which is what cut
> () is
> for).
I'll keep this in mind. But doesn't findInterval() need sorted input?
For this specific problem I would now use vec2RGB as suggested below.
>> col.class[is.na(col.class)]="#FFFFFF"
>
> Do you know how many are NAs?
No, this differs for the individual cases.
>> col.class.rgb=col2rgb(col.class)
>
> Assuming that you have full control of col.class, this should work.
It does:
> str(col.class.rgb)
int [1:3, 1:84000] 214 226 238 214 226 238 214 226 238 214 ...
- attr(*, "dimnames")=List of 2
..$ : chr [1:3] "red" "green" "blue"
..$ : NULL
>> grid.sp$band1=col.class.rgb["red",]
>
> Here you are overwriting the input object - do you need to?
> Wouldn't it be
> cleaner to start a fresh object?
As I do not need the data from the input object any more, I tried to
save memory by reusing this object (as these can become fairly
large). Could this be done better in a different way?
>> grid.sp$band2=col.class.rgb["green",]
>> grid.sp$band3=col.class.rgb["blue",]
>> grid.sp$band4[]=256
>
> Should the alpha band run 0-255 or 1-256? col2rgb() uses 0-255 (and
> can
> also add an alpha column).
I actually determined 256 by trial and error (the alpha channel
option from col2rgb() gives 255). With 256, writeGDAL() produced an
image with a transparent alpha channel. Assumed that the higher
number represents complete transparency.
The script results in the following SpatialGridDataFrame:
> str(grid.sp)
Formal class 'SpatialGridDataFrame' [package "sp"] with 6 slots
..@ data :'data.frame': 84000 obs. of 4 variables:
.. ..$ band1: int [1:84000] 214 214 214 214 214 214 214 214 214
214 ...
.. ..$ band2: int [1:84000] 226 226 226 226 226 226 226 226 226
226 ...
.. ..$ band3: int [1:84000] 238 238 238 238 238 238 238 238 238
238 ...
.. ..$ band4: num [1:84000] 256 256 256 256 256 256 256 256 256
256 ...
..@ grid :Formal class 'GridTopology' [package "sp"] with 3
slots
.. .. ..@ cellcentre.offset: Named num [1:2] 655050 278050
.. .. .. ..- attr(*, "names")= chr [1:2] "x_coord" "y_coord"
.. .. ..@ cellsize : Named num [1:2] 100 100
.. .. .. ..- attr(*, "names")= chr [1:2] "x_coord" "y_coord"
.. .. ..@ cells.dim : Named int [1:2] 350 240
.. .. .. ..- attr(*, "names")= chr [1:2] "x_coord" "y_coord"
..@ grid.index : int(0)
..@ coords : num [1:2, 1:2] 655050 689950 278050 301950
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : NULL
.. .. ..$ : chr [1:2] "x_coord" "y_coord"
..@ bbox : num [1:2, 1:2] 655000 278000 690000 302000
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "x_coord" "y_coord"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
.. .. ..@ projargs: chr "+init=world:CH1903"
>> #writeGDAL(grid.sp, graph.file,drivername = "GTiff",type =
>> "band",options="ALPHA=YES")
>
> The two versions are equivalent.
I know. I first suspected the error in writeGDAL() and therefore
tried to achieve the same result by directly using the "underlying"
function.
> Most likely, the cause of your problems
> is type = "band", which is not a valid type, being passed to
> underlying
> GDAL functions as "Unknown". My guess would be that you mean "Byte".
>
> The full list of GDAL data types is on www.gdal.org (down now for
> me, so
> no specific link), and from the R prompt as:
>
> rgdal:::.GDALDataTypes
That's the one! When I checked out writeGDAL I obviously
misinterpreted the type option. As it worked from the command line, I
did not suspect this option to be wrong.
> Finally, wouldn't it be easier to use the vec2RGB() function in
> rgdal, see
> the example using the meuse.grid data set at the foot of the
> example on
> the help page?
Yes, could be used instead of the "misused" cut (gives the same
result). Another reinvention from my side ...
Thanks a lot for your help and time.
Best regards,
Simon
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
