Ускорьте растр :: извлечение с весами в г

Question

Я хочу извлечь точное среднее значение растровых значений из экстента области, определенной многоугольником в r.Это работает, используя raster::extract с опцией weights=TRUE.Однако эта операция становится чрезмерно медленной при работе с большими растрами, и функция не распараллеливается, поэтому beginCluster() ... endCluster() не ускоряет процесс.Мне нужно извлечь значения для ряда растров, примеры которых приведены здесь как r, r10 и r100.Есть ли способ ускорить это в r, или есть альтернативный способ сделать это в GDAL?

r <- raster(nrow=1000, ncol=1000, vals=sample(seq(0,0.8,0.01),1000000,replace=TRUE))
r10 <- aggregate(r, fact=10)
r100 <- aggregate(r, fact=100)

v = Polygons(list(Polygon(cbind(c(-100,100,80,-120), c(-70,0,70,0)))), ID = "a")
v = SpatialPolygons(list(v))

plot(r)
plot(r10)
plot(r100)
plot(v, add=T)

system.time({
precise.mean <- raster::extract(r100, v, method="simple",weights=T, normalizeWeights=T, fun=mean)    
})
user  system elapsed 
0.251   0.000   0.253 
> precise.mean
      [,1]
[1,] 0.3994278    


system.time({
  precise.mean <- raster::extract(r10, v, method="simple",weights=T, normalizeWeights=T, fun=mean)    
})

user  system elapsed 
7.447   0.000   7.446 

precise.mean
      [,1]
[1,] 0.3995429

Answer 1

В конце я решил проблему, используя gdalUtils, работающий непосредственно с жестким диском.

Я использовал команду gdalwarp(), чтобы уменьшить разрешение растра до r10, 100. Затем gdalwarp(), чтобы увеличитьразрешение результирующего растра до исходного разрешения r.Затем gdalwarp() с cutline= "v.shp", crop_to_cutline =T для маскировки растра к вектору v.А затем gdalinfo() в сочетании с subset(x(grep("Mean=",x))) для извлечения средних значений.Все это было упаковано в цикл foreach() %dopar% для обработки ряда растров и разрешения.

Несмотря на сложность и, вероятно, не такую ​​точность, как extract::raster, он справился со своей задачей.

Answer 2

На самом деле он должен работать быстрее, если вы сначала вызовете beginCluster (затем функция занимается распараллеливанием).Еще лучше было бы использовать версию 2.7-14, которая имеет намного более быструю реализацию.В настоящее время он находится на рассмотрении в CRAN, но вы также можете получить его здесь: https://github.com/rspatial/raster