<p>Based on the replies, I tested two options: the use of <a href="https://cran.r-project.org/package=tdigest" rel="nofollow noreferrer"><code>tdigest</code> package</a> and the built-in parallelization routine from <code>terra</code> package (<code>cores</code> parameter). The t-Digest construction algorithm, by <a href="https://arxiv.org/abs/1902.04023v1" rel="nofollow noreferrer">Dunning et al., (2019)</a>, uses a variant of 1-dimensional k-means clustering to produce a very compact data structure that allows accurate estimation of quantiles. I recommend to use the <code>tquantile</code> function, which reduced by third the processing time with the tested dataset.</p>
<p>For those who were thinking about <code>foreach</code> parallelization, there are no simple solution to run foreach loop with terra objects. For such tasks, I'm still using the good old raster package. <a href="https://stackoverflow.com/questions/67445883/terra-package-returns-error-when-try-to-run-parallel-operations">It is a planned update but not on short term - see here</a>. More details below.</p>
<h1>Toy dataset</h1>
<pre><code>library(terra)
# load elevation coming with the terra pakage
r &lt;- rast( system.file(&quot;ex/elev.tif&quot;, package=&quot;terra&quot;) )

plot(r)

# number of iteration
n_it &lt;- 20

# With `stats::quantile()` function
start_time &lt;- Sys.time()
for (i in 1:n_it){
  ra &lt;- aggregate(r, 2 , fun = function(x) quantile(x, probs = .5, na.rm = T))
}
end_time &lt;- Sys.time()
print(end_time-start_time)
</code></pre>
<blockquote>
<p>Time difference of 6.013551 secs</p>
</blockquote>
<h1>With <code>tdigest::tquantile()</code></h1>
<pre><code>library(tdigest)

start_time_tdigest &lt;- Sys.time()
for (i in 1:n_it){
  ra_tdigest &lt;- aggregate(r, 2 , fun = function(x) tquantile(tdigest(na.omit(x)), probs = .5))
  }
end_time_tdigest &lt;- Sys.time()
print(end_time_tdigest-start_time_tdigest)
</code></pre>
<blockquote>
<p>Time difference of 1.922526 secs</p>
</blockquote>
<p>As suspected by Martin, the use of the <code>cores</code> parameter in the <code>terra:aggregate</code> function did not improve the processing time:</p>
<h1><code>stats::quantile()</code> + parallelization</h1>
<pre><code>start_time_parallel &lt;- Sys.time() 
for (i in 1:n_it){
ra_tdigest_parallel &lt;- aggregate(r, 2 , fun = function(x) quantile(x, probs = .5, na.rm = T), cores = 2) 
} 
end_time_parallel &lt;- Sys.time() 
print(end_time_parallel-start_time_parallel)
</code></pre>
<blockquote>
<p>Time difference of 8.537751 secs</p>
</blockquote>
<h1><code>tdigest::tquantile()</code> + parallelization</h1>
<pre><code>tdigest_quantil_terra &lt;- function(x) {   
require(tdigest)   
tquantile(tdigest(na.omit(x)), probs = .5) }
        
start_time_tdigest_parallel &lt;- Sys.time() for (i in 1:n_it){
ra_tdigest_parallel &lt;- aggregate(r, 2 , 
fun = function(x, ff) ff(x), cores = 2 , 
ff = tdigest_quantil_terra)     
} 

end_time_tdigest_parallel &lt;- Sys.time() 
print(end_time_tdigest_parallel-start_time_tdigest_parallel)
</code></pre>
<blockquote>
<p>Time difference of 7.520231 secs</p>
</blockquote>
<p>In a nutshell:</p>
<blockquote>
<p>1                tdigest 1.922526 secs</p>
<p>2          base_quantile 6.013551 secs</p>
<p>3       tdigest_parallel 7.520231 secs</p>
<p>4 base_quantile_parallel 8.537751 secs</p>
</blockquote>
Terra R - Speed up aggregate() of raster data with custom function
<p>I would like to use <code>aggregate</code> function from the <code>terra</code> <code>R</code> package to aggregate raster with a quantiles approach as aggregation function. Here below, I used the <code>quantile</code> function from <code>R base</code> to compute 50th percentile (i.e. the median) using a raster from the local package directory. I chose the 50th percentile for comparison with median but my goal is indeed to compute other quantile(s)...</p>
<pre><code>library(terra)
# load elevation coming with the terra pakage
r &lt;- rast( system.file(&quot;ex/elev.tif&quot;, package=&quot;terra&quot;) )

plot(r)

# number of iteration
n_it &lt;- 20

# with a custom function
start_time &lt;- Sys.time()
for (i in 1:n_it){
ra &lt;- aggregate(r, 2 , fun = function(x) quantile(x, probs = .5, na.rm = T))
}
end_time &lt;- Sys.time()
</code></pre>
<p>It took my computer approx. 6 secs to do it 20 times.</p>
<pre><code>print(end_time-start_time)
</code></pre>
<blockquote>
<p>Time difference of 6.052727 secs</p>
</blockquote>
<p>WHen I ran the same <code>aggregate</code> run with the median built-in function, it took approx. 40 times less time to perform the very same 20 iterations!</p>
<pre><code># with a built-in function
start_time &lt;- Sys.time()
for (i in 1:n_it){
  ra &lt;- aggregate(r, 2 , fun = median)
}
end_time &lt;- Sys.time()
print(end_time-start_time)
</code></pre>
<blockquote>
<p>Time difference of 0.1456101 secs</p>
</blockquote>
<p>As I would like to compute other percentiles than the 50th, could someone provide some advises to speed up <code>aggregate</code> when using custom functions?</p>