Какой самый эффективный способ использовать movingFun в больших временных рядах растров? - PullRequest
Новая
       38

Какой самый эффективный способ использовать movingFun в больших временных рядах растров?

0 голосов
/ 06 июня 2018

Мне нужно сгладить большой временной ряд, и я использую функцию movingFun из пакета 'растр'.Я протестировал несколько вариантов на основе предыдущих сообщений (см. Мои варианты ниже).Первые 2 работают, но работают очень медленно при использовании реальных данных (все временные ряды MOD13Q1 для всей Австралии).Поэтому я попытался вариант 3 и потерпел неудачу.Я был бы признателен, если бы кто-то мог помочь указать, что не так в этой функции.У меня есть доступ к памяти, я использую сервер RStudio Server с оперативной памятью 700 ГБ, но я не уверен, что будет лучшим подходом для выполнения этой работы.Заранее спасибо.

a) с помощью movingFun и overlay 

library(raster)
r <- raster(ncol=10, nrow=10)
r[] <- runif(ncell(r))
s <- brick(r,r*r,r+2,r^5,r*3,r*5)
ptm <- proc.time()
v <- overlay(s, fun=function(x) movingFun(x, fun=mean, n=3, na.rm=TRUE, circular=TRUE)) #works
proc.time() - ptm

   user  system elapsed 
  0.140   0.016   0.982

b) с созданием функции и с использованием clusterR.Я думал, что это будет быстрее, чем (а). 

fun1 = function(x) {overlay(x, fun=function(x) movingFun(x, fun=mean, n=6, na.rm=TRUE, circular=TRUE))}

beginCluster(4)
ptm <- proc.time()
v = clusterR(s, fun1, progress = "text")
proc.time() - ptm
endCluster()
   user  system elapsed 
  0.124   0.012   4.069

в) Я нашел этот документ , написанный Робертом Дж. Хеймансом, и я попытался (и не смог) написать функцию какописано в виньетках.Я не могу полностью выполнить все шаги в этой функции, поэтому терпит неудачу.

smooth.fun <- function(x, filename='', smooth_n ='',...) { #x could be a stack or list of rasters
  out <- brick(x)
  big <- ! canProcessInMemory(out, 3)
  filename <- trim(filename)
  if (big & filename == '') {
    filename <- rasterTmpFile()
  }
  if (filename != '') {
    out <- writeStart(out, filename, ...)
    todisk <- TRUE
  } else {
    vv <- matrix(ncol=nrow(out), nrow=ncol(out))
    todisk <- FALSE
  }

  bs <- blockSize(out)
  pb <- pbCreate(bs$n)

  if (todisk) {
    for (i in 1:bs$n) {
      v <- getValues(out, row=bs$row[i], nrows=bs$nrows[i] )
      v <- movingFun(v, fun=mean, n=smooth_n, na.rm=TRUE, circular=TRUE)
      out <- writeValues(out, v, bs$row[i])
      pbStep(pb, i)
    }
    out <- writeStop(out)
  } else {
    for (i in 1:bs$n) {
      v <- getValues(out, row=bs$row[i], nrows=bs$nrows[i] )
      v <- movingFun(v, fun=mean, n=smooth_n, na.rm=TRUE, circular=TRUE)
      cols <- bs$row[i]:(bs$row[i]+bs$nrows[i]-1)
      vv[,cols] <- matrix(v, nrow=out@ncols)
      pbStep(pb, i)
    }
    out <- setValues(out, as.vector(vv))
  }
  pbClose(pb)
  return(out)
}

s <- smooth.fun(s, filename='test.tif', smooth_n = 6, format='GTiff', overwrite=TRUE)

 Error in .local(.Object, ...) : 
  `/path-to-dir/test.tif' does not exist in the file system,
and is not recognised as a supported dataset name.

1 Ответ

0 голосов
/ 08 июня 2018

Это решение, которое я нашел, благодаря моему коллеге.Он рассчитывает каждый год (из 23 файлов) за 20 минут.Возможно, что-то и улучшится, но на этом этапе я рад, что могу выполнять работу всего за 20 минут в год.

Так что здесь я работаю 5 лет одновременно, используя пакет foreach.Затем цикл for создает массив из 6 файлов одновременно;помните, что мне нужно было 3-месячное скользящее окно, в 16-дневном наборе данных MOD13Q1, это 6 файлов.Затем цикл вычисляет средние значения в массиве, используя ColMeans, создает пустой растр, присваивает средние значения новому растру и сохраняет его.Обратите внимание, что мы воссоздали параметр circular функции movingFun.Итак, среднее число 1-й даты делается на основе последних дат того же года. 

require(raster)
require(rgdal)
library(foreach)
library(doParallel)

rasterOptions(maxmemory = 3e10, chunksize = 2e10)

ip <- "directory-with-grids"
op <- "directory-where-resuls-are-being-saved"

years = c(2000:2017)   

k <- 6    # moving window size
k2 <- floor((k-1)/2)
slot <- 0

# determine clusters
cl <- makeCluster(5, outfile = "") # make worker prints visible
registerDoParallel(cl)

foreach(j = 1:length(years), .packages=c("raster")) %dopar% {
  ip1 = paste(ip, years[j],sep='/')
  ndvi.files <- list.files(ip1, pattern = 'ndvi.*tif$',full.names = T) 
  nfiles <- length(ndvi.files)

  for (n in (1-(k-1)):nfiles) {
    i <- (n + k2 - 1) %% nfiles + 1
    print(ndvi.files[i])
    r <- raster(ndvi.files[i])
    if (slot == 0) {
      win <- matrix(data = NA, nrow = k, ncol = r@nrows * r@ncols)
    }
    slot <- slot %% k + 1
    win[slot,] <- getValues(r)
    if (n > 0) {
      o <- raster(extent(c(xx,xx,xx ,xx))); res(o)=c(xx,xx) # your extent and resolution
      crs(o) <-'+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0'
      o[] <- colMeans(win)
      o[o<0] <- NA
      # write out m as the nth raster
      fname = paste(names(r),'smoothed',sep='_')
      out.dir =  file.path(op, paste(years[j], sep='/'))
      dir.create(out.dir,showWarnings = FALSE)
      out.path = file.path(out.dir, fname)
      writeRaster(o, out.path, format="Geotiff", overwrite=T,  datatype='INT2S')
    }
  }
}

stopCluster(cl)
...