Интерполяция больших растровых рядов в R

Question

У меня есть 36-летние месячные температурные оценки сетки/растра, которые я бы хотел преобразовать в ежедневные оценки. На данный момент я устанавливаю ежемесячные оценки в середине месяца и простую линейную интерполяцию. Для этого я пытаюсь использовать растровые :: calc и stats :: approx, описанные в этом question. Однако, при этом, я получаю следующее сообщение об ошибке:

Error in is.infinite(v) : default method not implemented for type 'list'

Ниже приведен код, который мы надеемся, обеспечивает имитацию рода воссоздать проблему. Я думаю, что проблема заключается в том, как NA обрабатываются, потому что бит конца q_interp (где ни один растр не установлен в NA). Тем не менее, я не совсем уверен, что делать с этой информацией.

library(raster) 

#The parameters of the problem 
num_days = 9861 
months_num = 324 
num_na = 191780 

#generate baseline rasters 
r <- raster(nrows=360, ncols=720); 
values(r) <- NA 
x <- sapply(1:months_num, function(...) setValues(r, runif(ncell(r)))) 

#make them a stack 
s = stack(x) 

#define what x coordinates the rasters refer to (e.g. loosely convert monthly to daily). Probably not the most elegant solution in the world. 
num_day_month = c(31,28,31,30,31,30,31,31,30,31,30,31) 
days = as.character(seq(as.Date('1989/01/01'), as.Date('2015/12/31'), by = 'day')) 
months = as.character(seq(as.Date('1989/01/01'), as.Date('2015/12/01'), by = 'month')) 
months = substr(months, 1,nchar(months)-3) 
mid_points = as.vector(lapply(months, function(x) grep(x,days,value =T)[round(length(grep(x,days,value =T))/2)])) 
mp_loc = days %in% mid_points 
#output is the monthly mid points on the daily scale 
mp_day_locs = (1:length(days))[mp_loc] 

#make some of the cells NA throughout the whole span. In the actual dataset, the NAs generally represent oceans. 
s[sample(ncell(s), num_na)] = NA 

#a function to interpolate 
interp_row <- function(base_indexes, value_vector, return_indexes, rule_num =2) { 
    nnn = length(value_vector) 
    if (any(is.na(value_vector))) { 
    return(rep(NA, nnn)) 
    } else { 
    return(approx(x = base_indexes, y= value_vector, xout = return_indexes, rule=rule_num)$y) 
    } 
} 

#this is the function call that causes the error to be thrown 
s_interp = calc(s, function(y) interp_row(base_indexes = mp_day_locs, value_vector = y, return_indexes = 1:length(days),rule_num = 2)) 

#Now make a without NAs-- seems to work 
#generate baseline rasters 
r <- raster(nrows=360, ncols=720); 
values(r) <- NA 
x <- sapply(1:months_num, function(...) setValues(r, runif(ncell(r)))) 
#make them a stack 
q = stack(x) 
q_interp = calc(q, function(y) interp_row(base_indexes = mp_day_locs, value_vector = y, return_indexes = 1:length(days),rule_num = 2)) 

Answer 1

Проблема (насколько я могу видеть) - это вектор длины возврата, созданный, если любое значение равно NA. В вашем случае длина его такая же, как входной вектор:

# length = 324 
nnn = length(value_vector) 
return(rep(NA, nnn)) 

Однако возвращение вектор создан в правильном случае (без NA) намного больше (суточные значения):

#length = 9861 
return(approx(x = base_indexes, y= value_vector, xout = return_indexes, rule=rule_num)$y) 

Что касается как я знаю, два случая возвращения должны иметь одинаковую длину. Попробуйте изменить свою функцию следующим путем установки rep(NA, length(return_indexes)):

interp_row <- function(base_indexes, value_vector, return_indexes, rule_num =2) { 
    nnn = length(value_vector) 
    if (any(is.na(value_vector))) { 
    return(rep(NA, length(return_indexes))) 
    } else { 
    return(approx(x = base_indexes, y= value_vector, xout = return_indexes, rule=rule_num)$y) 
    } 
} 

Примечание: Ваш код кажется довольно медленно. Одним из быстрых решений является использование функции clusterR() для ускорения вашего кода. Эта функция позволяет использовать многоядерную обработку некоторых растровых функций, таких как calc. Тип ?clusterR или посмотрите here.