gdal error 1 (execution halted)
Dear Simon,
On Mon, 17 Mar 2008, sj wrote:
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).
Scripted R fails by design on error. Interactive R also fails, but hands control back to the prompt, instead of terminating R, because the interactive user is expected to sort the error out. The reason the output seemed to be OK is that the default "Unknown" type happened to make a 4-band RGBA file - by the way, we don't know whether this was at the stage of creating the GDAL dataset, or saving it, because the file gets written twice, once to a temporary location, next to its destination.
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.
No, findInterval() is an excellent and very fast function for this purpose. Only vec= needs to be sorted, but that is the same in cut() and not an unreasonable requirement.
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?
When debugging, it is never a good idea to overwrite. Further, R uses advanced memory management and a garbage collector, so you are not guaranteed immediate reductions in memory usage by overwriting - it will depend on when garbage collection happens.
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.
A string value not in rgdal:::.GDALDataTypes results in the type being set to "Unknown", which then has different behaviour fro different drivers. In the next rgdal release, users of writeGDAL will be forced to give a valid type (but not users of the underlying GDAL functions. Roger
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
Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, 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