Как переклассифицировать пространственный растр на основе кадра данных - PullRequest
Новая
       10

Как переклассифицировать пространственный растр на основе кадра данных

0 голосов
/ 24 сентября 2018

Я пытаюсь выполнить расчет для каждой ячейки пространственного растра.

На самом деле растр представляет собой цифровую модель рельефа и очень большой (1280000000 ячеек), я использую высоту прилива, чтобы рассчитать пропорцию появления каждой ячейки (под воздействием воздуха).

Пример:

Данные прилива 

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)

Растр 

require(raster)

r1 <- raster(matrix(seq(-4, 1.5, 0.5), ncol = 3))

Функция для расчета экспозиции 

expo <- list()

for(i in 1:ncell(r1)){
    depth <- df - r1[i]
    expo[[i]] <- length(depth[depth < 0])/length(depth)
    }

Преобразование обратно в растр 

r2 <- raster(matrix(unlist(expo), ncol = ncol(r1), byrow = T))

Это занимает много времени на большом растре, мне интересно, кто-нибудь может помочь ускорить его.Я пытался написать функцию для использования с raster :: calc, но не получил ее на работу.

Спасибо

1 Ответ

0 голосов
/ 24 сентября 2018

Вот подход, который вы можете попробовать

Я переименовал df в v, просто чтобы прояснить, что это числовой вектор, а не (и не должен быть) data.frame.

v <- df
lv <- length(v)

f <- function(i) {
    depth <- v - rep(i, lv)
    mean(depth < 0) 
}

Хитрость при отладке calc заключается в проверке функции с одной ячейкой.например, 

f(r1[2])

Теперь используйте его 

x <- calc(r1, f)

Или в одной строке как: 

x <- calc(r1, function(i) mean((v-rep(i,lv)) < 0) )

Чтобы ускорить больше, вы можете попробовать 

library(compiler)
ff <- cmpfun(f)
x <- calc(r1, ff)

Или 

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)

Если r1 не очень велико, возможно, вы можете ускорить его с помощью 

r2 <- setValues(r1, sapply(values(r1),  fcpp))

Или еще лучше: 

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)))
...