extracting time series data from a raster brick of AVHRR satellite data

classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

extracting time series data from a raster brick of AVHRR satellite data

janverbesselt
hi guys!

I am trying to optimise the extracting of a time series of NDVI data for
specific locations from a .hdr (envi file) containing 648 images. At the
moment I am doing the following:

> data <- readGDAL(file_in, offset=cOffset, region.dim=cRegion)
GIMMS_8208_FLOAT_BIL_australia has GDAL driver ENVI
and has 476 rows and 570 columns
>
> system.time(b <- brick(file_in))
   user  system elapsed
  0.171   0.001   0.173
> nlayers(b)
[1] 648
> # layerNames(b)
> xy <- cbind(148.15125,-35.65043)
> sp <- SpatialPoints(xy)
> system.time(data <- extract(b, sp))
   user  system elapsed
150.605   1.022 152.752
> ndvi <- ts(as.vector(data[1,]),start=c(1982,1),frequency=24)

This works correctly but the extract() is extremely slow, i.e. 152seconds
for one pixel.
Is there a faster/better way of doing this?

Alternatively I could use:
cOffset <- c(0,0)
cRegion <- c(1,1)

# read file and info
fileinfo <- GDALinfo(file_in)
bands <- fileinfo[[3]]
data <- readGDAL(file_in, offset=cOffset, region.dim=cRegion)
which is fast but then I still need to define the specific location (in
lat/long) when using readGDAL().

Thanks for your help.
Best,
Jan


> sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] C/en_US.UTF-8/C/C/C/C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] bfast_1.3         forecast_2.19     fracdiff_1.3-2    tseries_0.10-25
quadprog_1.5-4    MASS_7.3-14       strucchange_1.4-4
 [8] sandwich_2.2-7    zoo_1.7-2         raster_1.9-1      rgdal_0.6-33
 sp_0.9-84

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: extracting time series data from a raster brick of AVHRR satellite data

Martin
This post was updated on .
what I do to extract a pixel is:

>row=40
>col=58

>test = mybrick[row,col]

>s=ts(test, start=2001, end=c(2003, 36), frequency=36)


this is extremely fast, but you need to know the row and col of the pixel...
Martin Brandt

Department of Geography and
Regional Research
UZA II, Althanstr. 14
1090 Wien, AUSTRIA
http://geooekologie.univie.ac.at
Reply | Threaded
Open this post in threaded view
|

Re: extracting time series data from a raster brick of AVHRR satellite data

Martin
sorry, it has to be:

>row=40
>col=58

>test = mybrick[row,col]

btw: does anyone know how to get the row, col information in QGIS or GRASS?
Martin Brandt

Department of Geography and
Regional Research
UZA II, Althanstr. 14
1090 Wien, AUSTRIA
http://geooekologie.univie.ac.at
Reply | Threaded
Open this post in threaded view
|

Re: extracting time series data from a raster brick of AVHRR satellite data

janverbesselt
I have done the following:

1. read a small subset of the modis image time series (2000-2011, i.e. 262 images) in and save as a raster stack.
2. below read in the .grd file as a brick with 262 layers and dim = 17,13 containing integer values
3. now I am trying to extract time series per pixel for further processing, this however is very slow.
It takes more than 30sec for one time series. I have tried modis[row,col] but this is not faster.

Any advice? Should I do this via Python or can this be done via R but more efficient?


> ## read in Brick
> modis <- brick("modis.grd")
> nlayers(modis)
[1] 262
> dim(modis)
[1]  17  13 262
>
> (system.time(test <- extract(modis,1)))
   user  system elapsed
 35.985   0.084  36.194
> (system.time(test <- modis[1,1]))
   user  system elapsed
 35.847   0.076  36.051

Reply | Threaded
Open this post in threaded view
|

Re: extracting time series data from a raster brick of AVHRR satellite data

Robert Hijmans
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
Reply | Threaded
Open this post in threaded view
|

Re: 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] <[hidden email]> 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



To unsubscribe from extracting time series data from a raster brick of AVHRR satellite data, click here.

Reply | Threaded
Open this post in threaded view
|

Re: extracting time series data from a raster brick of AVHRR satellite data

janverbesselt
In reply to this post by Robert Hijmans
I have also made the 'modis.grd' data available via:
## Public Link
download.file("http://dl.dropbox.com/u/8150174/raster/modis.grd","m.grd")
download.file("http://dl.dropbox.com/u/8150174/raster/modis.gri","m.gri")

More info here:

Thanks, Jan
Reply | Threaded
Open this post in threaded view
|

Re: extracting time series data from a raster brick of AVHRR satellite data

Robert Hijmans
In reply to this post by janverbesselt


On Thu, Jul 28, 2011 at 1:53 AM, janverbesselt [via R-sig-geo] <[hidden email]> wrote:
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?


Ah, bug, happens when values of b are in memory, not when they are on disk. Fixed in devel version. Thanks, Robert

 

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?


Oh, that not right either (readAll with a RasterStack does not set the inMemory flag correctly).

modis <- brick("modis.grd")
modis <- readAll(modis) 

 
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.


Precisely by sending these type of reports. 
Thanks, Robert.

 

On Wed, Jul 27, 2011 at 9:23 PM, Robert Hijmans [via R-sig-geo] <[hidden email]> 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



To unsubscribe from extracting time series data from a raster brick of AVHRR satellite data, click here.




To unsubscribe from extracting time series data from a raster brick of AVHRR satellite data, click here.