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

[Alexander Shenkin]

Hi All,

Thanks for everyone's input, especially Forrest Stevens'.  I've put my 
solution up here: https://gist.github.com/ashenkin/7fceb77e78efc33961a8 .

Thanks,
Allie

On 3/11/2016 3:06 PM, Forrest Stevens wrote:
> Using the raster package's zonal() function maybe be better suited for 
> your situation Allie, depending on the size of the underlying 
> elevation raster and how efficient you really need it to be, this 
> approach could be used (tacked on to your reproducible example):
>
> ##Use zonal from the raster package:
> r <- raster(meuse.grid["dist"])
> b <- rasterize(my_blpi, r)
> my_blpi[["a"]] <- zonal(r, b, fun='mean')[,2]
> spplot(my_blpi["a"])
>
> As somewhat of an aside, this requires you to rasterize your transects 
> (the second line in the code above) which if being done over very 
> large areas and fine grained scales can be quite inefficient using the 
> raster package.  In the past I've called the GDAL rasterize utility 
> directly via a system() call to get modest gains in this context.
>
> Hope this helps!
>
> Forrest
>
> On Fri, Mar 11, 2016 at 6:24 AM Edzer Pebesma 
> <edzer.pebesma at uni-muenster.de <mailto:edzer.pebesma at uni-muenster.de>> 
> wrote:
>
>     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 <mailto: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
>
>     _______________________________________________
>     R-sig-Geo mailing list
>     R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
>     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>


	[[alternative HTML version deleted]]