5

I want to perform a trend analysis for each pixel over it's time dimension. Therefore I wrote a loop over each pixel doing the calculation but this is extremely slow. I am sure there must be a faster way to perform the computation.

Also the dimensions of the file are 1440,720,480. Is this a doable size for R? The problem is that i cannot find the statistical function implemented in IDL which I otherwise use.

gridfile = 'test.nc'
data = brick(gridfile)

require(fume)

mk = matrix(numeric(0), 720,1440) 

for(i in 1:72) {
  for(j in 1:144) {
  temp = mkTrend(as.vector(data[i,j,]))
  mk[i,j] = temp$`Sen's Slope`
  }
}
8
  • 1
    Do you feel obliged to use the Sen slope estimator? It is hideously inefficient (requiring about 10^5 calculations for each pixel in your case). There are far faster robust methods (some requiring only around 10 calculations per pixel), which would lead to a speedup of three or four orders of magnitude. Consequently you could use your computational budget to learn far more about these data than just making a map of slopes.
    – whuber
    Mar 18, 2014 at 18:53
  • I am looking for significant trends in climatological data and the mann kendall test seems to be pretty standard within this area since the data does not need to fit specific requirements, with the exception of not being autocorrelated. That's why I am trying this adapted mann kendall test in this R package which corrects the effects of autocorrelation. Do you perhaps know of similar trend analysis techniques?
    – Spamiad
    Mar 18, 2014 at 21:49
  • Unfortunately your data are likely to exhibit strong autocorrelation and resorting to a nonparametric test isn't going to cure that problem. Moreover, it appears that mkTrend may be doing far more calculation than you need (the Sen slope estimate requires all 480*479/2 pairs of values to be compared). Be aware, too, of the serious multiple comparisons problem associated with conducting 720*1440 separate (but highly interdependent) tests at once. These considerations all suggest using simpler and more revealing procedures.
    – whuber
    Mar 18, 2014 at 23:25
  • For the multiple comparisons problem I assume the simplest way would be to apply a Bonferroni correction. The Mann-Kendall test applied in my example, using that specific R package, adjusts the p-values according to autocorrelation within the data: A modified Mann-Kendall trend test for autocorrelated data. I will look for alternative methods identifying the trend, although for this type if analysis I can only think of regression models.
    – Spamiad
    Mar 19, 2014 at 9:21
  • 1
    Going back to your original issue, though: it looks like if you avoid mkTrend you might see a speedup of up to several orders of magnitude. There are loads of ways to identify trends: consider exploring Cross Validated for possibilities. For instance, there are robust generalizations of the "three point method" I describe at stats.stackexchange.com/a/35717. E.g., you take a few points near the beginning, a few near the middle, and a few at the end; compute each of their medians (componentwise), and fit a line through the three resulting points. You can even test for linearity.
    – whuber
    Mar 19, 2014 at 14:03

1 Answer 1

3

Package raster provides a function calc that executes a function over all raster pixels, returning the raster with results; it can deal with large rasters (out-of-memory). raster also provides cluster/multi-core functionality, look into function beginCluster; it's meant to work with calc.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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