Raster parallel process in R

May 10, 2021 · 5 min read · raster parallel snow foreach ·

By Fidel Maureira Sotomayor

Install libraries if there are no already in your computer.

 1packages = c('raster','terra', 'snow', 'parallel', 'foreach','doParallel','rasterVis', 'mapview','ggplot2', 'lubridate','pryr')
 2out = sapply(packages ,  function(p){
 3  print(p)
 4  if(require(p,character.only = TRUE)==FALSE) {
 5    install.packages(p)}
 6   
 7    require(p,character.only = TRUE)  })
 8
 9
10
11my_function = function(r, na.rm=TRUE, fun=mean) {
12  fun(getValues(r), na.rm=na.rm)
13}

Download files, alternative is possible use the function raster::getData. More details on https://worldclim.org/data/worldclim21.html

1temp <- tempfile()
2tempd <- tempdir()
3#download worldclim climate data
4download.file('https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_10m_tavg.zip',temp, mode="wb")
5#raster::getData()
6unzip(temp, exdir=tempd)
7#dir(tempd) # explore the temp directory

Explore the data

1file_names = list.files(tempd,'wc2.1_10m_tavg')
2tmean <- stack(file.path(tempd, file_names))
3mapview(tmean[[7]])

## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 2332800 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  2332800 '

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : 104
## projected point(s) not finite

## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : 104
## projected point(s) not finite

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

Paraleling

Create cluster and load functions. Alternative you can create your own scrip file with function and just load this file.

1cl <- makeCluster(4, type='SOCK', outfile='log.txt')
2registerDoParallel(cl)
3out=clusterEvalQ(cl,library(raster))
4env= where('my_function')
5clusterExport(cl, list('my_function'),envir = env)
6#out =clusterEvalQ(cl, source('my_functions.R'))

We going to use two way of paralleling in R with the snow/parallel and foreach libraries.

1system.time({
2  clusterExport(cl,list('tmean'))
3  mean_planet = parLapply(cl,c(1:12), function(l){
4    my_function(tmean[[l]])
5  })
6  mean_planet = unlist(mean_planet)
7})

##    user  system elapsed 
##    0.00    0.00    1.51

1system.time({
2  mean_planet2 = foreach(l = c(1:12),
3                         .combine = rbind,
4                         .errorhandling = 'remove') %dopar%
5    my_function(tmean[[l]])
6})

##    user  system elapsed 
##    0.04    0.02    1.37

1results = data.frame(layer_name = names(tmean),
2                     month =  month(1:12,label=TRUE, abbr = TRUE),
3                     parLapply_out =mean_planet,
4                     foreach_out = mean_planet2)

The first plot show the global mean temperature for each month. The second figure show that the two paralleling outcome the same results

1ggplot(results,aes(x=month,y=parLapply_out, col='parLapply_out'))+
2  geom_point()+labs(y='Mean temperature (°C)', col='', x='')+
3  geom_point(aes(y=foreach_out,col = 'foreach_out'))+theme_minimal()

Split large raster

We download the raster and make a split overlapped extent. The overlapped extent is only for the use of focus function, any other process should use a non-overlapped extent! Also is possible use the function focus in raster. but it is slower function.

1file_raster = 'https://www.dropbox.com/s/onp3f75pab7z3j9/intro-1508529851.jpg?dl=1'
2
3r_raster = stack(file_raster)
4e =extent(r_raster)
5extent_r =list(r1 =extent(0, floor(e@xmax/2)+2, 0, floor(e@ymax/2)+2),
6               r2 = extent(ceiling (e@xmax/2), e@xmax, floor (e@ymax/2), e@ymax),
7               r3 =  extent(0, floor (e@xmax/2)+2,  floor (e@ymax/2), e@ymax),
8               r4 =  extent(ceiling (e@xmax/2), e@xmax,  0, floor (e@ymax/2)+2))

We call the parLapply from a lapply environment, so we need export elements from this environment (env).

 1tmean_smooth = lapply(c(1:3) , function(l){
 2  env= where('l')
 3  clusterExport(cl, list('extent_r','r_raster','l'),envir=env)
 4  r_parts = parLapply(cl,seq_along(extent_r),function(part){
 5    r = r_raster[[l]]
 6    r_part <- crop(r, extent_r[[part]])
 7    #focal(r_part,w=matrix(1,3,3), fun='mean',na.rm=TRUE)
 8    terra::focal(r_part,w=matrix(1,3,3), fun='mean',na.rm=TRUE)
 9    #r_part
10  })
11  
12  m1 <- mosaic(r_parts[[1]], r_parts[[2]],r_parts[[3]],r_parts[[4]], fun=mean)
13  
14  m1
15})
16stopCluster(cl)
17#r_smooth = stack(tmean_smooth[[1]],tmean_smooth[[2]],tmean_smooth[[3]])
18r_smooth = do.call(stack,tmean_smooth)
19
20plotRGB(r_smooth)
21abline(v=390,h=219,col='green')

Versus original

1plotRGB(r_raster)

##Visualization

more info in https://oscarperpinan.github.io/rastervis/

1data(volcano)
2volcano = raster(volcano)
3library(viridis)

## Loading required package: viridisLite

1rasterVis::gplot(volcano) + 
2  geom_tile(aes(fill = (value)))  +
3  viridis::scale_fill_viridis(direction = -1, 
4                              na.value='#FFFFFF00') +
5  labs(fill='Altitud (m) ', x='E',y='N')+
6  theme_bw()

1levelplot(volcano,contour=TRUE)

1xyplot(wc2.1_10m_tavg_03~wc2.1_10m_tavg_07|cut(x, 4), data = tmean, auto.key = list(space='right'))