Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

enter image description here

> 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!

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

1 Answer 1

up vote 4 down vote accepted

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) 
share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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