[R-sig-Geo] extract mean raster value from polygon

[Alexander Shenkin]

Hello All,

I've been working to be able to make elevation profiles from a DEM along 
a swath, rather than just a line (thanks to Forrest Stevens for the help 
so far).  To that end, I've made a function, create_perp_buffers, that 
creates polygons perpendicular to, and along, a transect (see graphic 
below).

Now, I would like to be able to intersect each polygon with the 
SpatialGridDataFrame DEM, and get an average elevation for each 
polygon.  I've tried using over(), but it just hangs forever (my actual 
DEM is quite large, a mosaic from SRTM).  Can anyone suggest a good way 
to do this?

Thanks,
Allie






create_perp_buffers <- function(x1, y1, x2, y2, grid_dist, slice_width, 
proj4string = "+init=epsg:28992 
+towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812") {
     xdiff <- x2 - x1
     ydiff <- y2 - y1
     lineangle <- atan(ydiff/xdiff)
     dx = cos(lineangle)
     dy = sin(lineangle)
     left_angle = lineangle - 90
     right_angle = lineangle + 90
     midpoints = data.frame(x = seq(from = x1, to = x2, by = 
dx*grid_dist), y = seq(from = y1, to = y2, by = dy*grid_dist))
     begin.coord = data.frame("x" = cos(left_angle)*slice_width + 
midpoints$x, "y" = sin(left_angle)*slice_width + midpoints$y)
     end.coord = data.frame("x" = cos(right_angle)*slice_width + 
midpoints$x, "y" = sin(right_angle)*slice_width + midpoints$y)

     l <- vector("list", nrow(begin.coord))

     for (i in seq_along(l)) {
         l[[i]] <- Lines(list(Line(rbind(begin.coord[i,], 
end.coord[i,]))), as.character(i))
     }

     sl <- SpatialLines(l)

     names(begin.coord) <- c("begin_x", "begin_y")
     names(end.coord) <- c("end_x", "end_y")
     sldf <- SpatialLinesDataFrame(sl, 
data.frame("lineID"=1:nrow(begin.coord), begin.coord, end.coord))
     proj4string(sldf) = proj4string

     blpi <- gBuffer(sldf, byid=TRUE, id=sldf$lineID, width = grid_dist/2)

     return(blpi)
}

my_blpi <- create_perp_buffers(180000, 331500, 181000, 332500, 100, 100)

data(meuse.grid)
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE
proj4string(meuse.grid) <- 
CRS(paste("+init=epsg:28992","+towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812"))
plot(meuse.grid, axes = T)
plot(my_blpi, col="red", add=T)

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20160311/228a7824/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: igeciice.png
Type: image/png
Size: 8315 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20160311/228a7824/attachment.png>