0

First some context: I want to run a time-series analysis function that is only available in R for a continental sized multi-band time-series raster, which is stored in disk in GeoTIFF format, that was created translating a VRT file to a GeoTIFF with GDAL using parameter INTERLEAVE='PIXEL' (That parameter should store all band values for each pixel continuously on the disk an make it easy and fast to read.)

The function only needs the time-series to work, so no spatial dependency. And the file has 380 bands, each of continental size, summing up 420GB on disk space.

The problem: I used the raster package to read the file as a stack, but when I try to get the values for a pixel, R just stucks and appears to be processing, but this operation should be fast. I tried using the [] notation, using the extract function and none worked. (Even if I try to get only 3 bands from a small RGB file, R seems slow.)

Using gdallocationinfo it is extremely fast, and with Python and GDAL module it's also incredibly fast.

So shouldn't it be fast in R too? Am I doing something wrong or missing a function designed for that? Doesn't the Raster package use rgdal internally?

But most importantly, is there a way to do it using rgdal then?

  • Is it tiled? raster::extract is very slow in this case because it reads line-by-line, requiring multiple repeated reads of tiles. That needs to be fixed. You could use the lower level rgdal tools more directly I think, or possibly rgdal2. For R/raster/extract, the only fix to to write to an untiled file afaik. (Happy to explore this but it will have to be tomorrow. ) – mdsumner Jul 21 '16 at 16:26
  • I put details here of what I see here: gist.github.com/mdsumner/8404668a41b47b1dacb5d8b8f326acc1 For a single band I would read the entire thing in and do the extraction, but that probably won't help for your multi-band. – mdsumner Jul 21 '16 at 16:35
  • I'm not sure if it's tiled, but probably is. (Would appreciate a hint about how to find that information :)). Well, I imagined that I'd probably have to use rgdal, but I never worked with it before and am a bit lost. Finally, thank you for pointing rgdal, I will look into it. Looking forward for your exploration tomorrow :) – Hemeligur Jul 21 '16 at 17:52
1

I've tried out a bit and found that:

  1. RasterBricks, as mentioned by RobertH's answer, do work and are more user-friendly and easy to use;
  2. Rgdal methods like readGDAL also work, but with more parameters it's a little bit less user-friendly;

So which option should one use?

According to my tests (on my 420GB GeoTiff with dimensions of 18660x21592 and 374 bands) Rgdal is faster. Maybe due to less overhead of a higher level library such as the Raster-package.

Here are my results, using system.time and replicate:

With brick:

> modis_ndvi_ts_brick <- brick("../data/pa_br_mod13q1_ndvi_250_2000_2016.tif")
> system.time(replicate(100, modis_ndvi_ts_brick[7000,7000]))
   user  system elapsed 
 94.024   5.468  99.562

While with rgdal:

> system.time(replicate(100, readGDAL("../data/pa_br_mod13q1_ndvi_250_2000_2016.tif", offset=c(7000,7000), region.dim=c(1,1))))
#some printed output from readGDAL here
   user  system elapsed 
 88.752   5.400  94.213 

So as you can see Rgdal is slightly faster. For those whose this does not make a difference, I recommend using RasterBrick, it's simpler. But for those whose, like me, are struggling to create a high performance code in which every millisecond matters: Use Rgdal


Note: I'm sorry for the not exactly reproducible data, but one could try it out with other data and post here if their results differ somehow

  • RasterBrick builds on rgdal so that could explain the difference. Then again, the replicate test is not relevant, because you normally would extract many cells in one step, which is a different problem. – Robert Hijmans Jul 22 '16 at 23:06
  • I agree, maybe a test with various randomly chosen cells (or maybe blocks of cells) would be more appropriate. – Hemeligur Jul 23 '16 at 11:37
0

As you have a single file, you should create a RasterBrick. That should make things faster as it could indeed benefit from the by pixel interleave. By creating a RasterStack you create a list of RasterLayers, i.e. you treat each "band" as a separate file.

  • I didn't use a RasterBrick because I read on various places that it loads all the information to memory, which with my 420GB file wouldn't be possible. But it seems that you're suggesting that the RasterBrick is intelligent enough to know how to work with this. I will try and post back here. Thank you – Hemeligur Jul 22 '16 at 18:37
  • I just tried and it works, thank you. I will post an answer regarding the differences of using brick and rgdal, for future comers. – Hemeligur Jul 22 '16 at 20:38
  • "I read on various places that it loads all the information to memory" Really? That is patently wrong. – Robert Hijmans Jul 22 '16 at 23:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.