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