3

I'm working on spatial polygon data frame (spdf) dataset. This dataset contains a time series of monthly climate data. What I'm trying to achieve is to convert the spdf to raster stack with 1000m resolution for further statistical analysis. I've writen an R code but is is very slow and it took ages to convert one column.I would appreciate if any of you could suggest tips to make it faster.

hru<-readOGR("E:\\Tade\\HRU\\ubn_merge.shp",layer="ubn_merge") # spatial polygon
spdf<-merge(hru,spdf.2000,by.x="HRU",by.y="HRU",all.x=T,sort=F) # spdf nrow=565 ncol=375
# convert sp to raster
hru.ras<-raster(hru,resolution=1000) # raster hru shape to 1km
for (i in 1:length(spdf){
  et.ras<-rasterize(spdf,hru.ras,field=paste("m",1,sep="")) # rasterize
  et.stack<-stack(et.stack,et.ras)
}

Thanks

Tade
  • 43
  • 5
  • 1
    In my experience you're in for a world of pain using R and the `rasterize()` function to do this task. From my testing the conversion from vector to raster is much more manageable when using the `gdal_rasterize` utilty directly. I call it using a system call, e.g. `system(paste0("gdal_rasterize -burn 1 -l census_buffer ", census_buffer_path, " /tmp/census_mask.tif"))` etc. which obviously assumed the `gdal` utilities are in your path. – Forrest R. Stevens Nov 04 '15 at 21:04
  • @ForrestR.Stevens thank you. – Tade Nov 05 '15 at 09:05

2 Answers2

3

As Forrest says (and you have experienced), rasterize is a bit slow, but you can do much better than what you have now. You do not need to use a loop at all:

r <- raster(spdf, resolution=1000)
et.ras <-rasterize(spdf, r, field=paste0("m",1:ncol(spdf)))

This will create a single RasterLayer with a Raster Attribute Table. To create a RasterStack, do:

s <- deratify(et.ras)

If you were to use a loop, use rasterize only once, the get the polygon IDs, and then use subs for the actual variables of interest.

And of course this saves you the pain of external dependencies.

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Thanks. The speed improved dramatically ( few hours to few minutes). – Tade Nov 05 '15 at 09:07
  • As always Robert, you provide a lot of food for thought (in such a good way). I'd not thought of the approach you outlined, and indeed I think it would improve performance in some of my use cases over issuing multiple calls to `gdal_rasterize`. Thank you! – Forrest R. Stevens Nov 05 '15 at 14:59
2

As of March 2020, you can now consider this new package.

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

fasterize()

Works exactly as 

rasterize()

But 100-1000 times faster. I think this can help a lots of people.

Paul Roub
  • 36,322
  • 27
  • 84
  • 93
Antoine C.
  • 21
  • 2