Libraries needed

library(raster)
library(stars)
library(sf)
library(foreach)
library(doParallel)
library(iterators)
library(mapview)

Intro

To have a first glimpse on parallel computation in R I suggest the chapter on parallel computation of the “R Programming for Data Science” book. https://bookdown.org/rdpeng/rprogdatascience/parallel-computation.html
BLAS : For linux users (but also for windows with a little more effort) it could be usefull to check if the BLAS (Basic Linear Algebra Subprograms) parallelization could be optimized on the machine. Because the default libraries have not been properly optimized, you could see large increases in speed, for code that is dependent on linear algebra computations, by switching to a different BLAS library. The traditional open source solutions for this have been ATLAS, GotoBLAS (no longer developed), and more recently OpenBLAS (GotoBLAS2 fork). This post provides some tips: http://brettklamer.com/diversions/statistical/faster-blas-in-r/
to install OpenBLAS e ATLAS on Ubuntu run:
install OpenBLAS
sudo apt-get install libopenblas-base
install ATLAS
sudo apt-get install libatlas3-base liblapack3
To switch between the different blas libraries use:
sudo update-alternatives –config libblas.so.3-x86_64-linux-gnu
To check the currently used LAPACK and BLAS within an R session, run:

if(.Platform$OS.type!="windows"){
  # These functions will not work on windows OS
  La_library()
  extSoftVersion()["BLAS"]
}
##                                              BLAS 
## "/usr/lib/x86_64-linux-gnu/openblas/libblas.so.3"

Scope

Extract common statistics of the values of a raster inside certain polygons (zonal statistics) using parallel computation to speed up the process.
The foreach package is an attempt to make parallel computing much simpler by providing a parallel for-loop construct for R. It requires a parallel system (backend) like doParallel or doMPI. We will only use doParallel.
Create clusters:
First we create the clusters and then we register them to the {doParallel} package. On windows the only option is the so-called “SNOW functionality”, but on linux the so-called “Multi Core” functionality is also availabe.

if(.Platform$OS.type=="windows"){ # If we are on windows
  # Snow-like functionality
  cl<-makeCluster(8) # adjust according to the number of available cores (CPU)
  registerDoParallel(cl)
  
}else{ # If we are on unix
  # Multi Core functionality
  registerDoParallel(cores = 8)
  
}


getDoParWorkers()
## [1] 8
getDoParName()
## [1] "doParallelMC"
getDoParVersion()
## [1] "1.0.15"

Read data

We read in some cadastral parcels for the province of Bozen/Bolzano as sf polygon object. The original file is available here: http://geocatalogo.retecivica.bz.it/geokatalog/

polygon=readRDS("polygons.RDS")

We read the global radiation (year sum) for the province of Bozen/Bolzano as a stars object. The original file is available here: http://geocatalogo.retecivica.bz.it/geokatalog/

rad=read_stars("EURAC_RADGLOB_YEAR_25m.tif")

let’s plot the solar radiation and the polygons

map=mapview(as(rad,"Raster"),layer.name="Solar Radiation")+mapview(polygon)
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 26152192 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  26152192 '
#map

to save memory on the website we will just show a screenshot of the actual mapview map