[R-sig-Geo] rasterize in parrallel

[John Baumgartner]

I like using gdalUtils for this type of thing. It's often much faster to
work with files on disk using the gdal utilities than to use the raster
package, and may not even warrant parallelisation.

library(maptools)
library(raster)
library(gdalUtils)

# Read poly and assign CRS
polyg <- readShapePoly('polyg.shp')
proj4string(polyg) <- "+init=epsg:2154"

# Create the template raster
grain <- 50
r_template <- raster(polyg, res=c(50, 50))

# For each field (name), write out the template raster (to the working
directory),
# and then burn in the field's values. See ?gdal_rasterize and
# http://www.gdal.org/gdal_rasterize.html for details. Setting output_raster
# to TRUE means that the created rasters are returned to R and can be
# conveniently stacked.
s <- stack(lapply(names(polyg), function(nm) {
  f <- paste0('polyg_', nm, '.tif')
  suppressWarnings(writeRaster(r_template, f))
  gdal_rasterize(src_datasource='polyg.shp', dst_filename=f, a=nm,
                 output_Raster=TRUE)
}))

This took 56 sec on my system.

It's a simple step to then parallelise, which, for me, reduces the time
taken to 16 sec.

library(parallel)
nsc <- 4
cl <- makeCluster(nsc)

# Export required objects to each processes
clusterExport(cl, c('r_template', 'polyg'))

# Load the libraries on each process
clusterEvalQ(cl, sapply(c('raster', 'gdalUtils'), require, char=TRUE))

# Send to processes
system.time(s <- stack(parLapply(cl, names(polyg), function(nm) {
  f <- paste0('polyg_', nm, '.tif')
  suppressWarnings(writeRaster(r_template, f))
  gdal_rasterize(src_datasource='polyg.shp', dst_filename=f, a=nm,
                 output_Raster=TRUE)
})))


Cheers,
John




On Tue, Jul 14, 2015 at 8:03 AM, Jérome Mathieu <jerome.mathieu at upmc.fr>
wrote:

> Dear all,
>
> I need to rasterize the same shapefile according to different fields stored
> in the @data of the shapefile. So I need to loop trough the fields (named
> sc1 to sc27) to rasterize upon each field. With a regular sequential loop
> it takes a very long time (much longer than in ArcMap), so I would like to
> parallelize it, but I didn't succeed so far.
>
> I tried :
>
> # read shapefile
>  library(maptools)
>  polyg<-readShapePoly("D:\\polyg.shp")
>  polyg at proj4string <- CRS("+init=epsg:2154")
>
> library(doParallel)
> cl <- makeCluster(4)
> registerDoParallel(cl)
>
> grain<-50
> nsc<-4
>
> foreach(idxsc = 1 : nsc , .packages="raster", polyg=polyg) %dopar% {
>
>
> eval(parse(text=paste("RRpolyg_",idxsc,"<-raster(polyg,res=c(grain,grain))",sep="")))
> # create raster template
> eval(parse(text=paste("projection(RRpolyg_",idxsc,") <-
> proj4string(polyg)",sep="")))      # projection
>
> eval(parse(text=paste("RRpolyg_",idxsc,"<-rasterize(polyg,RRpolyg_",idxsc,",field=polyg at data
> $sc",idxsc,")",sep="")))
> # rasterization
>
> }
>
> stopCluster(cl)
>
> and I get the following error :
>
> Error in { :
>   task 1 failed - "unable to find an inherited method for function 'raster'
> for signature '"numeric"'"
>
>
> the data are here
>
> http://we.tl/7qqegfLoBA
>
>
> Any help would be greatly appreciated !
>
> Jerome
>
>
> > sessionInfo()
> R version 3.1.3 (2015-03-09)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
> Running under: Windows 7 x64 (build 7601) Service Pack 1
>
> locale:
> [1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
> LC_MONETARY=French_France.1252 LC_NUMERIC=C
> LC_TIME=French_France.1252
>
> attached base packages:
> [1] parallel  stats     graphics  grDevices utils     datasets  methods
> base
>
> other attached packages:
>  [1] rgeos_0.3-11     rgdal_0.9-2      maptools_0.8-36  lattice_0.20-31
> foreign_0.8-63   raster_2.3-40    sp_1.1-0         doParallel_1.0.8
> iterators_1.0.7  foreach_1.4.2
>
> loaded via a namespace (and not attached):
> [1] codetools_0.2-11 compiler_3.1.3   grid_3.1.3       tools_3.1.3
>
>
>
>
>
>
>
> --
> Jérôme Mathieu
> Université Pierre & Marie Curie
> Institute of Ecology and Environmental Science (Paris)
>
> Bât. A - 7ème Etage, porte 715
> 7 quai St.-Bernard
> F-75252 Paris Cedex 05
>
> tel: 01 44 27 34 22
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]