[R-sig-Geo] rasterize in parrallel

[John Baumgartner]

(Note that the lines reading in the polygon and assigning the CRS are
actually redundant, since with gdal_rasterize we refer directly to the
shapefile's file path, and the .shp already had the appropriate CRS
assigned.)

On Tue, Jul 14, 2015 at 10:02 AM, John Baumgartner <johnbaums at gmail.com>
wrote:

> 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]]