¿Cómo rasterizar los polígonos espaciales en R?

Estoy tratando de extraer los valores de batimetría de mi área de interés de una capa ráster de batimetría mundial utilizando la función 'rasterizar' en el paquete {sp}.

* Ediciones: Encontré la función 'extraer' que parece ser más lo que estoy buscando.

Esto es lo que he hecho hasta ahora:

> 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)
> plot(subarea0, add=TRUE)

> 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 

No entiendo el resultado. ¿Por qué obtengo 2 colores para cada uno de mis polígonos? ¿Qué quieren decir?

¿Cómo puedo obtener finalmente el contorno de profundidad de batimetría? ¿Tiene esto algo que ver con mi resolución o con cambiar las dimensiones?

* Ediciones: ok, ahora he hecho lo siguiente:

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

v es una lista y no estoy seguro de cómo / bajo qué forma volver a vincular esta información con mi polígono espacial ...

¡Gracias!

Su línea ras_sub0<-rasterize(subarea0,raster_bath)solo toma el número de índice de los polígonos y lo asigna a los valores del ráster.

Si desea solo la intersección de su polígono y el ráster:

subarea0_bathy <- intersect(raster_bath, subarea0)

Actualización : como señala @GodinA, parece que intersect () a veces no devuelve un ráster que tiene la extensión completa del polígono. Para evitar esto, puede cruzarse con un ráster un poco más grande que el 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) 

Aquí hay otra solución, usando una función básica de subconjunto en raster.

# create a raster with the dimensions that you are interested in.
r <- raster(extent(subarea0)+1) # +1 to increase the area
res(r) <- 0.1 # you can change the resolution here. 0.1 can be about 10x10 km if you are close to the equator.
crs(r) <- projection(subarea0) # transfer the coordinate system to the raster

# Now you can add data to the cells in your raster to mark the ones that fall within your polygon.
r[subarea0,] <- 1 # this an easy way to subset a raster using a SpatialPolygons object

# check your raster
plot(r)