I need to somehow combine my Spatialpolygons with my depth contour, but not sure how do to so. Maybe create a SpatialPolygonsDataframe?

This is what I've done so far:

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

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

data("Bathymetry")#World bathymetric data in library (marelac)
[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 (not sure if this is the best way to achieve what I want)
    extent(raster_bath) <- extent(subarea0)#Transform the extend of my raster to the extend of my SpatialPolygons
v <- extract(raster_bath, subarea0)#Extract data from my Raster Layer for the locations of my SpatialPolygons

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


I apologize in advance, looks like I've been having a hard time to explain clearly what I was hoping to get. I don't want to summarize my depth info for each of my polygon but rather provide a detailed bathymetry within my area of interest. I am definitely not sure if what I did up to now is the way to go, but ultimately what I would like to achieve is to end-up with a plot of my polygon with my depth information. Something like this enter image description here (thanks to Simbamangu), but without the depth info outside of my polygon. Basically, I just want the bathymetry info comprise in my overall area at the finest scale possible. Should I first transform my SpatialPolygons into one and then into some grid format to then be able to clip with my raster layer?

Thank you!

  • Could you be more precise about what you mean by "combine"? What is the objective of this operation? What formats are the inputs in? What format is needed for output? What software are you willing to use? And how does this question differ from the extremely similar one (involving the same code) you posted at gis.stackexchange.com/questions/25049 ? – whuber May 9 '12 at 15:17
  • definitely helps if you put what software you are using in your title and in a tag. – Stev_k May 9 '12 at 15:18
  • 1
    @whuber, I was asked to start a new question... I am simply trying to get the depth contour for my area of interest (and only my area of interest, nothing else around) (subarea0-Spatialpolygons object in R) from a RasterLayer (raster_bath). When I used 'intersect' as proposed in my previous question, I do get the depth contour, but not only within the boundaries of subarea0. I get a bounding box with some depth contour falling outside of my subarea0. Hopefully this makes a bit more sense. Ultimately, I would like to be able to plot using ggplot2 my depth contour + polygon. Thank you. – GodinA May 9 '12 at 15:59
  • In the first question you asked how to extract bathymetry values for your area of interest - I understand that the next step (this question) is to GENERATE contours for the area of interest, then CLIP them for your area of interest so they don't go beyond the boundaries of subarea0? – Simbamangu May 9 '12 at 17:05
  • Yes @Simbamangu, this is exactly what I would like :) – GodinA May 9 '12 at 18:22

Slightly different approach to first question ... you'll need to create SpatialLines (or a SpatialLinesDataFrame) and work with that.

library(maptools) # Contour tools
library(marelac) # bathymetry data
library(raster) # only used for displaying Bathymetry here
library(rgeos) # for intersections

cont <- contourLines(Bathymetry$x, Bathymetry$y, Bathymetry$z) # generate contours
c.lines <- ContourLines2SLDF(cont) # create SpatialLines DF

#Have a look:
b.rast <- raster(Bathymetry) 
plot(b.rast) # plot world bathymetry
lines(c.lines) # our SpatialLines overlap and look OK

# Clip to your boundary
subarea0 <- readShapePoly("subarea0.shp") 
subarea0.lines <- gIntersection(c.lines, subarea0) # clip the contours to subarea0
  • @GodinA - curious to know if this worked? – Simbamangu May 21 '12 at 10:19
  • sorry for my late reply. I was away from my emails. Thank you very much for your answer. This indeed worked and was exactly what I was looking for. Thank you again! – GodinA Jun 8 '12 at 14:54

v is a list. There is one list element for each polygon. Each list element can have many values (one for each raster cell that it covers). Therefore, you cannot simply combine the data to the SpatialPolygons (and you are not saying why you want to do that anyway). To get a single value for each polygon, you can use extract with a function (mean would be the natural choice)

vv <- extract(raster_bath, subarea0, fun=mean)

or use lapply with v, something like this:

vvv <- lapply(v, mean, na.rm=TRUE) 

turn that into a data.frame and use function SpatialPolygonsDataFrame to create that type of object.

  • Thanks for your help @RobertH, I added details in my question, hopefully, this clarifies a little what I am hoping to do. – GodinA May 10 '12 at 12:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.