9

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

  • 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 '16 at 7:43
  • I put summary in body of question. – Diogo Caribé Oct 6 '16 at 14:53
  • Not summary, just print - the info I asked for us not tgere – mdsumner Oct 6 '16 at 20:15
  • Sorry, I put the print. – Diogo Caribé Oct 7 '16 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 '16 at 7:23
17

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:

  • Very good answer! – Diogo Caribé Oct 7 '16 at 9:54
  • Do you not just set n as the number of cores on the machine? – Sam Jan 25 '17 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. – Guzmán Jan 25 '17 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 '17 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 – Diogo Caribé Apr 19 '18 at 12:49

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.