[R-sig-Geo] gdal error 1 (execution halted)

[Roger Bivand]

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