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

[sj]

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