[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