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

Edzer Pebesma edzer.pebesma at uni-muenster.de
Fri Mar 11 12:09:09 CET 2016


Allie, thanks for the reproducible example, which runs after loading
packages sp and rgeos. You can get the aggregate & plot directly by

a = aggregate(meuse.grid["dist"], my_blpi, mean)
spplot(a)

Package raster may have routines that do the same thing faster for large
raster data. In raster, ``aggregate'' has a very different meaning.


On 11/03/16 11:07, Alexander Shenkin wrote:
> 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)
> 
> 
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/
Spatial Statistics Society http://www.spatialstatistics.info

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20160311/67e7eb90/attachment.bin>


More information about the R-sig-Geo mailing list