I am trying to perform a calculation on every cell of a spatial raster.
In reality the raster represents a digital elevation model and is very large (1280000000 cells), I am using tide height to calculate the proportion that each cell is emerged (exposed to the air).
Example:
Tide data
df <- c(3.879, 4.078, 4.211, 4.252, 4.204, 4.077, 3.872, 3.588, 3.259,
2.883, 2.48, 2.065, 1.635, 1.199, 0.766999999999999, 0.339,
-0.0840000000000005,
-0.503, -0.906, -1.284, -1.649, -1.998, -2.326, -2.603, -2.801,
-2.959, -3.108, -3.237, -3.329, -3.353, -3.343, -3.303, -3.199,
-3.041, -2.803, -2.503, -2.173, -1.789, -1.348, -0.869000000000001,
-0.373, 0.141999999999999, 0.657999999999999, 1.207, 1.728, 2.226,
2.683, 3.055, 3.393, 3.655, 3.841, 3.956, 3.988, 3.938, 3.816,
3.63, 3.365, 3.047, 2.69, 2.292, 1.871, 1.433, 0.981999999999999,
0.524, 0.0759999999999996, -0.367, -0.805000000000001, -1.226,
-1.637, -2.036, -2.422, -2.741, -2.956, -3.137, -3.322, -3.481,
-3.593, -3.662, -3.727, -3.791, -3.79, -3.707, -3.557, -3.356,
-3.077, -2.732, -2.354, -1.962, -1.515, -1.035, -0.515000000000001,
0.00599999999999934, 0.532999999999999, 1.05, 1.563, 2.032, 2.462,
2.794, 3.098, 3.313)
Raster
require(raster)
r1 <- raster(matrix(seq(-4, 1.5, 0.5), ncol = 3))
Function to calculate exposure
expo <- list()
for(i in 1:ncell(r1)){
depth <- df - r1[i]
expo[[i]] <- length(depth[depth < 0])/length(depth)
}
Convert back to raster
r2 <- raster(matrix(unlist(expo), ncol = ncol(r1), byrow = T))
This takes a long time on a large raster, I'm wondering if anyone can help speed it up. I've tried to write a function to use with raster::calc but have not got it to work.
Thanks
Here is an approach you can try
I renamed df to v, just to make clear that it is a numeric vector, and not (and should not be) a data.frame.
v <- df
lv <- length(v)
f <- function(i) {
depth <- v - rep(i, lv)
mean(depth < 0)
}
The trick with debugging calc is to test the function with one cell. e.g.
f(r1[2])
Now use it
x <- calc(r1, f)
Or on one line as:
x <- calc(r1, function(i) mean((v-rep(i,lv)) < 0) )
To speed up more, you can try
library(compiler)
ff <- cmpfun(f)
x <- calc(r1, ff)
Or
library(Rcpp)
fcpp <- cppFunction('double flood(double x) {
// [[Rcpp::plugins(cpp11)]]
if (std::isnan(x)) return(NAN);
std::vector<double> v = {3.879, 4.078, 4.211, 4.252, 4.204, 4.077, 3.872, 3.588, 3.259, 2.883, 2.48, 2.065, 1.635, 1.199, 0.766999999999999, 0.339, -0.0840000000000005, -0.503, -0.906, -1.284, -1.649, -1.998, -2.326, -2.603, -2.801, -2.959, -3.108, -3.237, -3.329, -3.353, -3.343, -3.303, -3.199, -3.041, -2.803, -2.503, -2.173, -1.789, -1.348, -0.869000000000001, -0.373, 0.141999999999999, 0.657999999999999, 1.207, 1.728, 2.226, 2.683, 3.055, 3.393, 3.655, 3.841, 3.956, 3.988, 3.938, 3.816, 3.63, 3.365, 3.047, 2.69, 2.292, 1.871, 1.433, 0.981999999999999, 0.524, 0.0759999999999996, -0.367, -0.805000000000001, -1.226, -1.637, -2.036, -2.422, -2.741, -2.956, -3.137, -3.322, -3.481, -3.593, -3.662, -3.727, -3.791, -3.79, -3.707, -3.557, -3.356, -3.077, -2.732, -2.354, -1.962, -1.515, -1.035, -0.515000000000001, 0.00599999999999934, 0.532999999999999, 1.05, 1.563, 2.032, 2.462, 2.794, 3.098, 3.313};
unsigned lv = v.size();
unsigned depth = 0;
for (size_t i=0; i<lv; i++) {
depth += ((v[i] - x) < 0);
}
return (double(depth) / lv);
}')
x <- calc(r1, fcpp)
If r1 is not very large, you can perhaps speed it up further with
r2 <- setValues(r1, sapply(values(r1), fcpp))
Or better still:
library(Rcpp)
fcpp2 <- cppFunction('std::vector<double> flood(std::vector<double> x) {
// [[Rcpp::plugins(cpp11)]]
unsigned sizex = x.size();
std::vector<double> out(sizex);
std::vector<double> v = {3.879, 4.078, 4.211, 4.252, 4.204, 4.077, 3.872, 3.588, 3.259, 2.883, 2.48, 2.065, 1.635, 1.199, 0.766999999999999, 0.339, -0.0840000000000005, -0.503, -0.906, -1.284, -1.649, -1.998, -2.326, -2.603, -2.801, -2.959, -3.108, -3.237, -3.329, -3.353, -3.343, -3.303, -3.199, -3.041, -2.803, -2.503, -2.173, -1.789, -1.348, -0.869000000000001, -0.373, 0.141999999999999, 0.657999999999999, 1.207, 1.728, 2.226, 2.683, 3.055, 3.393, 3.655, 3.841, 3.956, 3.988, 3.938, 3.816, 3.63, 3.365, 3.047, 2.69, 2.292, 1.871, 1.433, 0.981999999999999, 0.524, 0.0759999999999996, -0.367, -0.805000000000001, -1.226, -1.637, -2.036, -2.422, -2.741, -2.956, -3.137, -3.322, -3.481, -3.593, -3.662, -3.727, -3.791, -3.79, -3.707, -3.557, -3.356, -3.077, -2.732, -2.354, -1.962, -1.515, -1.035, -0.515000000000001, 0.00599999999999934, 0.532999999999999, 1.05, 1.563, 2.032, 2.462, 2.794, 3.098, 3.313};
unsigned sizev = v.size();
for (size_t j=0; j<sizex; j++) {
if (std::isnan(x[j])) {
out[j] = NAN;
} else {
unsigned depth = 0;
for (size_t i=0; i<sizev; i++) {
depth += ((v[i] - x[j]) < 0);
}
out[j] = (double(depth) / sizev);
}
}
return(out);
}')
r2 <- setValues(r1, fcpp2(values(r1)))