12

I am converting vector to raster in R. However the process was too long. Is there possibility to put the script into multithread or GPU processing in order to do it more faster?

My script to rasterized vector.

r.raster = raster()
extent(r.raster) = extent(setor) #definindo o extent do raster
res(r.raster) = 10 #definindo o tamanho do pixel
setor.r = rasterize(setor, r.raster, 'dens_imov')

r.raster

class : RasterLayer dimensions : 9636, 11476, 110582736 (nrow, ncol, ncell) resolution : 10, 10 (x, y) extent : 505755, 620515, 8555432, 8651792 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

setor

class : SpatialPolygonsDataFrame features : 5419 extent : 505755, 620515.4, 8555429, 8651792 (xmin, xmax, ymin, ymax) coord. ref. : +proj=utm +zone=24 +south +ellps=GRS80 +units=m +no_defs variables : 6 names : ID,CD_GEOCODI, TIPO, dens_imov,area_m,domicilios1 min values : 35464, 290110605000001, RURAL, 0.00000003,100004,1.0000 max values : 58468, 293320820000042, URBANO, 0.54581673,99996, 99.0000

Print of setor enter image description here

7
  • Can you post summaries of setor and r.raster? I'd like to have some idea of the number of objects in setor and the dimensions of r.raster. just print them is fine
    – mdsumner
    Oct 6, 2016 at 7:43
  • I put summary in body of question. Oct 6, 2016 at 14:53
  • Not summary, just print - the info I asked for us not tgere
    – mdsumner
    Oct 6, 2016 at 20:15
  • Sorry, I put the print. Oct 7, 2016 at 0:55
  • Ah, disappointed I didn't think of this until I saw the print-out - make sure the raster's projection matches the polygons, it doesn't at the moment - try r <- raster(setor); res(r) <- 10; setor.r = rasterize(setor, r, 'dens_imov') - but also try, setting res(r) <- 250 first so you get an idea of how long the high-res version will take
    – mdsumner
    Oct 7, 2016 at 7:23

2 Answers 2

19

I tried to "parallelize" the function rasterize using the R package parallel in this way:

  1. split the SpatialPolygonsDataFrame object in n parts
  2. rasterize every part separately
  3. merge all the parts into one raster

In my computer, the parallelized rasterize function took 2.75 times less than the no-parallelized rasterize function.

Note: the code below download a polygon shapefile (~26.2 MB) from the web. You can use any SpatialPolygonDataFrame object. This is only an example.

Load libraries and example data:

# Load libraries
library('raster')
library('rgdal')

# Load a SpatialPolygonsDataFrame example
# Load Brazil administrative level 2 shapefile
BRA_adm2 <- raster::getData(country = "BRA", level = 2)

# Convert NAMES level 2 to factor 
BRA_adm2$NAME_2 <- as.factor(BRA_adm2$NAME_2)

# Plot BRA_adm2
plot(BRA_adm2)
box()

# Define RasterLayer object
r.raster <- raster()

# Define raster extent
extent(r.raster) <- extent(BRA_adm2)

# Define pixel size
res(r.raster) <- 0.1

BrazilSPDF

Figure 1: Brazil SpatialPolygonsDataFrame plot

Simple thread example

# Simple thread -----------------------------------------------------------

# Rasterize
system.time(BRA_adm2.r <- rasterize(BRA_adm2, r.raster, 'NAME_2'))

Time in my laptop:

# Output:
# user  system elapsed 
# 23.883    0.010   23.891

Multithread thread example

# Multithread -------------------------------------------------------------

# Load 'parallel' package for support Parallel computation in R
library('parallel')

# Calculate the number of cores
no_cores <- detectCores() - 1

# Number of polygons features in SPDF
features <- 1:nrow(BRA_adm2[,])

# Split features in n parts
n <- 50
parts <- split(features, cut(features, n))

# Initiate cluster (after loading all the necessary object to R environment: BRA_adm2, parts, r.raster, n)
cl <- makeCluster(no_cores, type = "FORK")
print(cl)

# Parallelize rasterize function
system.time(rParts <- parLapply(cl = cl, X = 1:n, fun = function(x) rasterize(BRA_adm2[parts[[x]],], r.raster, 'NAME_2')))

# Finish
stopCluster(cl)

# Merge all raster parts
rMerge <- do.call(merge, rParts)

# Plot raster
plot(rMerge)

BrazilRaster

Figure 2: Brazil Raster plot

Time in my laptop:

# Output:
# user  system elapsed 
# 0.203   0.033   8.688 

More info about parallelization in R:

6
  • Very good answer! Oct 7, 2016 at 9:54
  • Do you not just set n as the number of cores on the machine?
    – Sam
    Jan 25, 2017 at 15:28
  • @Sam I think it should work without problem but I don't know if it's better or not! I assumed that if I splitted the features in n parts equal to the number of cores maybe one of this parts could be easier to process and the core that processed it would be without use! However, if you have more parts than cores when one core finish processing one part it would take other part. But certainly, I'm not sure! Any help on this issue would be appreciated.
    – Guz
    Jan 25, 2017 at 15:43
  • i'm going to run some tests tonight. On a small shapefile (roughly 25km by 25km), rasterized to 50m, there is a tiny improvement in using n = 2,4 or 8 against n = 20, 30 or up to 50. I'll sub in a very large shapefile tonight and rasterize to 25m. Single core processing is 10hrs so we will see what different values of n do!! (n=50 is just under 1 hr)
    – Sam
    Jan 25, 2017 at 15:55
  • @Guzmán I am runnig the code again. However, it returened some error and a don't know why. Can you help me? Error in checkForRemoteErrors(val) : 7 nodes produced errors; first error: object 'BRA_adm2' not found Apr 19, 2018 at 12:49
0

try fasterize, no more need to parallelize

https://cran.r-project.org/web/packages/fasterize/vignettes/using-fasterize.html

1
  • I should note that I am getting very reasonable speeds with rasterizing vector data using the terra package. I have not benchmarked it but, the performance seems fairly comparable to fasterize, to the point that I needed to operationalize a rasterization process, for creating sample grids and masks, in a large scale analysis and just ended up using terra, after comparing with fasterize. There may be a measurable difference when rasterizing a large number of features and would be interesting to test. Apr 19, 2021 at 17:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.