CODEWITHSUNDEEP
r
Climate Vet
Climate Vet

Reputation: 1

Rasterize multiple SpatialPolygonsDataframe attributes

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!

Upvotes: 0

Views: 2019

Answers (1)

Spacedman
Spacedman

Reputation: 94277

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.

Upvotes: 3

Related Questions