[R-sig-Geo] plotting of large raster dataset fails

Etienne etiennesky at yahoo.com
Tue Mar 29 03:47:32 CEST 2011


Greetings,

I am trying to plot a large raster dataset (2159x925 = 2M pixels) using spplot() and it fails with the following error after about 30 seconds:

> r2=readGDAL("img/TM.burnpix.pnsc.2000-2006.freq.tif")
> image(r2)
> spplot(r2)
Error: segfault from C stack overflow

I can plot the same data fine with the image() function in less than a second.

The main problem comes from the fact that many pixels are NA's (I do not wish to include some in the analysis).  If I set the NAs to zero, I can plot the data but it takes a very long time (about 30 seconds), using 100% cpu and a lot of memory (1GB) - orders of magnitude higher than using the image() function.

Questions:

1)
Is this volume of data supported by the spplot() function, and if so is there any trick to get this working properly and faster?

2)
Can I work with NAs or must I convert them to another value (0,Nan ?)

Many thanks,
Etienne

The dataset has the following structure,:
=====
> str(r2)
Formal class 'SpatialGridDataFrame' [package "sp"] with 6 slots
  ..@ data       :'data.frame':	1997075 obs. of  1 variable:
  .. ..$ band1: num [1:1997075] 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] 290965 7751395
  .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  .. .. ..@ cellsize         : num [1:2] 30 30
  .. .. ..@ cells.dim        : int [1:2] 2159 925
  ..@ grid.index : int(0) 
  ..@ coords     : num [1:2, 1:2] 290965 355704 7751395 7779115
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] 290950 7751380 355719 7779130
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. ..@ projargs: chr " +proj=utm +zone=23 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0"
=====



More information about the R-sig-Geo mailing list