[R-sig-Geo] extracting time series data from a raster brick of AVHRR satellite data

[janverbesselt]

Robert,

I have tried the different options - and based on your example below with
the simulated raster brick everything works fine.
However when applying it to my dataset a few things/errors occur.

This is my dataset:
> b
class       : RasterStack
dimensions  : 17, 13, 262  (nrow, ncol, nlayers)
resolution  : 231.6564, 231.6564  (x, y)
extent      : 305554.7, 308566.3, 5713341, 5717279  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs
min values  : 2874  535 2526 3779 3626 4520 5747 4820  502 4962 ...
max values  : 7671 9126 7660 7181 8723 8217 8513 8581 8979 8854 ...

1/ when trying:
> nc <- writeRaster(b,'modis.nc',overwrite=TRUE) # ,datatype='INT2U'
Error in `dataType<-`(`*tmp*`, value = "FLT4S") :
  Cannot set datatype of a RasterStack

I tried setting the datatype to ,datatype='INT2U' but this gave the same
error. How could this be solved?

2/ When doing the following in a new R session:
modis <- stack("modis.grd")
modis <- readAll(modis)
> inMemory(modis)
[1] FALSE

somehow when loading the file from disk in a new R session and then doing
readAll() is does not load it into memory and when trying modis[1] is
remains slow. How can I make sure that the file loads into memory?

3/ therefore I currently opted for the following:

> system.time(m <- getValues(modis))
   user  system elapsed
  0.347   0.030   0.453

## create a time series
d <- as.vector(m[70,]) # extract cell 70 from the matrix m

Thanks for your help. I really like the 'raster' package so let me know how
I can help testing the package.


On Wed, Jul 27, 2011 at 9:23 PM, Robert Hijmans [via R-sig-geo] <
ml-node+6627415-691589920-304146 at n2.nabble.com> wrote:

> Jan,
>
> > Any advice? Should I do this via Python or can this be done via R but
> more efficient?
>
> Interesting problem, that I had not noticed before. I think this is because
> a loop that can perhaps be omitted or else perhaps needs to move to a C
> routine. I need to investigate that. For now, for speed, either read the
> values into memory (if possible; perhaps a spatial subset, via 'crop'), or
> perhaps use ncdf files (see example below).
>
> > library(raster)
> Loading required package: sp
> raster version 1.9-3 (27-July-2011)
> > b <- brick(nc=10, nr=10, nl=648)
> > b <- setValues(b, matrix(rep(1:100, 648), nc=648))
> > system.time(v <- b[1])
>    user  system elapsed
>       0       0       0
> >
> > d <- writeRaster(b, 'test.grd', overwrite=TRUE)
> > inMemory(d)
> [1] FALSE
> > system.time(v <- d[1])  # too slow
>    user  system elapsed
>   61.54    0.14   61.98
> >
> > d <- readAll(d)
> > inMemory(d)
> [1] TRUE
> > system.time(v <- d[1])  # OK
>    user  system elapsed
>       0       0       0
> >
> > e <- writeRaster(b, 'test.nc', overwrite=TRUE)
> Loading required package: ncdf
> > inMemory(e)
> [1] FALSE
> > system.time(v <- e[1]) # OK
>    user  system elapsed
>    0.01    0.00    0.02
> >
>
>
> Robert
>
>
> ------------------------------
>  If you reply to this email, your message will be added to the discussion
> below:
>
> http://r-sig-geo.2731867.n2.nabble.com/extracting-time-series-data-from-a-raster-brick-of-AVHRR-satellite-data-tp6622055p6627415.html
>  To unsubscribe from extracting time series data from a raster brick of
> AVHRR satellite data, click here<http://r-sig-geo.2731867.n2.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=6622055&code=amFudmVyYmVzc2VsdEBnbWFpbC5jb218NjYyMjA1NXwtMTAxMTIxMDUxNQ==>.
>
>


--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/extracting-time-series-data-from-a-raster-brick-of-AVHRR-satellite-data-tp6622055p6629146.html
Sent from the R-sig-geo mailing list archive at Nabble.com.