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

[Forrest Stevens]

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>
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
> > 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
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

	[[alternative HTML version deleted]]