(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 <[email protected]> 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 <[email protected]> > 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@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@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 >> [email protected] >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
