Rasterize multiple SpatialPolygonsDataframe attributes

Question

I need to 'rasterize' multiple attributes of a SpatialPolygonsDataframe. Then, I will multiply the values in each raster layer by the values of a given vector (e.g. 'gp'). I am using the wrld_simpl dataset from maptools for a start. wrld_simpl has an attribute called POP2005 and I will add more attrributes corresponding to population estimates for the years (2010:2100). I'm trying to write a loop or function so that I don't have to manually rasterize each attribute, then calculate their new values independently, and then repeat all these steps once and again.

data(wrld_simpl)  
gp <- seq(1,246)  

myraster <- raster(nrow = 572, ncol = 1440, xmn = -180, xmx = 180, ymn= -58, ymx = 85)

Because the rasterize() function does not allow me to select a specific attribute from the spatialPolygonsDataframe, I have used raster:::.polygonsToRaster() instead (which does exactly the same) to convert my attributes to a raster layer.

rastergp <- raster:::.polygonsToRaster(wrld_simpl, myraster, field = wrld_simpl$POP2005)  
stackraster <- stack(rastergp, gp)  
estimation <- calc(stackraster, fun = function(x) x[1]*x[2])

Does anyone have any suggestions about how to proceed? Thanks a lot!

Spacedman · Accepted Answer · 2012-07-23 16:41:50Z

3

Ummm you should cut down your script to the bare minimum to illustrate what you are trying to do. There's some stuff in there that isn't relevant and is distracting and getting in the way of us figuring it out.

What I think you are trying to do (and hard to tell when you use an undocumented hidden function raster:::.polygonsToRaster) is get the values of the attributes of polygons at a grid of locations? You can do this by overlaying a grid of points (rather than a raster) on your polygons.

I'll do it here for a smaller grid - your resolution will take 100 times longer:

 > myraster <- raster(nrow = 57, ncol = 144, xmn = -180, xmx = 180, ymn= -58, ymx = 85)
 > pts = SpatialPoints(xyFromCell(myraster,1:prod(dim(myraster))),proj4string=CRS(proj4string(wrld_simpl)))
 > overGrid = pts %over% wrld_simpl

Now overGrid is just a standard data frame with the values from wlrd_simpl at all the grid points. Its in the same order as your raster. Hence:

> POP2010 = myraster
> POP2010[] = overGrid$POP2010
> plot(POP2010)

will map the POP2010 raster. Compare:

> spplot(wrld_simpl,"POP2010")

So that's done all the overlay business. You can then make a stack from the columns of overGrid.

answered Jul 23, 2012 at 16:41

Spacedman

94.7k12 gold badges148 silver badges231 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

Climate Vet Over a year ago

Thanks for your input Spacedman. I have reduced the script to a minimum to illustrate what I need to do. What you propose gives basically the same result of what I am doing now. But my main problem is that I am creating one file at the time (from a total of 91) and I want to create a loop not to have to run the script 91 times. Any input in that matter will be greatly appreciated.

Collectives™ on Stack Overflow

Rasterize multiple SpatialPolygonsDataframe attributes

1 Answer 1

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related