# How to rasterize SpatialPolygons in R?

I'm trying to extract the bathymetry values of my area of interest from a world bathymetry raster layer using the 'rasterize' function in the {sp} package.

*Edits: I found the 'extract' function which seems to be more what I'm looking for.

This is what I've done so far:

``````> class(subarea0) #This is my area of interest (Eastern Canadian Arctic Sea)
[1] "SpatialPolygons"
attr(,"package")
[1] "sp"

> extent(subarea0)
class       : Extent
xmin        : -82.21997
xmax        : -57.21667
ymin        : 60.2
ymax        : 78.16666

library(marelac)
data("Bathymetry")#World bathymetric data in library (marelac)
names(Bathymetry);class(Bathymetry);str(Bathymetry)
[1] "x" "y" "z"
[1] "list"
List of 3
\$ x: num [1:359] -180 -179 -178 -177 -176 ...
\$ y: num [1:180] -89.5 -88.5 -87.5 -86.5 -85.5 ...
\$ z: num [1:359, 1:180] 2853 2873 2873 2873 2873 ...

raster_bath<-raster(Bathymetry)#Transformed into a raster layer
extent(raster_bath) <- extent(subarea0)#Transform the extend of my raster to the extend of my SpatialPolygons

>ras_sub0<-rasterize(subarea0,raster_bath)#rasterize my SpatialPolygons (*Edits: not the function that I need here, but I am still interested to learn what results mean)
Found 2 region(s) and 10 polygon(s)
> plot(ras_sub0)
``````

``````> ras_sub0
class       : RasterLayer
dimensions  : 180, 359, 64620  (nrow, ncol, ncell)
resolution  : 0.06964709, 0.0998148  (x, y)
extent      : -82.21997, -57.21667, 60.2, 78.16666  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84
values      : in memory
min value   : 1
max value   : 2
layer name  : layer
``````

I don't understand the result. Why am I getting 2 colors for each of my polygons? What do they mean?

How can I ultimately get bathymetry depth contour? Is this something to do with my resolution or changing the dimensions?

*Edits: ok, I've now done the following:

``````v <- extract(raster_bath, subarea0)#Extract data from my Raster Layer for the locations of my SpatialPolygons
``````

v is a list and I am not quiet sure how/under what form to rebind this info with my spatial polygon...

Thank you!

-
When you say 'get bathymetry depth contour' do you mean you want to also generate contours? – Simbamangu May 9 '12 at 7:13

Your line `ras_sub0<-rasterize(subarea0,raster_bath)` is just taking the index number of the polygons and assigning that to the values of the raster.

If you want just the intersection of your polygon and the raster:

``````subarea0_bathy <- intersect(raster_bath, subarea0)
``````

Update: As @GodinA notes, looks like intersect() sometimes doesn't return a raster that has the complete extent of the polygon! To get around this, you can intersect with a slightly larger raster than your original:

``````r2 <- raster() # create a generic raster
extent(r2) <- extent(subarea0) + 1 # add 1 degree of extent (0.5 each side)
r3 <- intersect(raster_bath, r2)
``````
-
Thank you @Simbamangu, I tried the "intersect" command, however when I plot the subarea0_bathy and my spatialpolygon (subarea0), they do not overlap completely, why? To answer your previous question, yes I would like to eventually also have my bathymetry contour. Any suggestions? – GodinA May 9 '12 at 12:10
Yes, looked at the output from the bathymetry data and a Tanzania polygon and I see the same effect. See my update above. For contours, it's as easy as `cont <- contour(r3)` then `lines(cont)` Isn't R great? – Simbamangu May 9 '12 at 12:44
Great, yes R can be fantastic! Thanks @Simbamangu, this works! However, I would like to only plot my Spatialpolygons with my depth contour. Maybe create a SpatialPolygonsDataframe? Using intersect, I get the info for the whole square encompassing my Spatialpolygons. I need to somehow combine my Spatialpolygons with my depth contour, but not sure how do to so. This is why I thought that the function 'extract' was a good start, but I end up with a list of depth that I am not sure how I can rebind with my spatialPolygons... Any ideas? Thanks again! – GodinA May 9 '12 at 13:14
What do you want to end up with? Contours (SpatialLines) that overlap with your subarea0? You might want to start a new question as this is now going beyond the original. – Simbamangu May 9 '12 at 13:37
Yes, I would like to only get my subarea0 with the depth contour for this define area: gis.stackexchange.com/questions/25112/… – GodinA May 9 '12 at 15:14