You have a digital elevation model (DEM), and you want to make an elevation profile between two or more points. You want that elevation profile to reflect an average elevation across a swath, and not just a single point from a line. Here's one way to do it in R. Thanks to Forrest Stevens and other folks from the R-sig-geo list.

make_elevation_profile_swaths.r

# Thanks to Forrest Stevens.  Some of the code here borrowed from him here: https://github.com/ForrestStevens/Scratch/blob/master/swath_slices.R

library(raster)
library(rgdal)
library(sp)
library(rgeos)
library(gtools)
library(ggplot2)
library(plyr)
library(zoo)

# we'll use this data, and use the "dist" column as if it were elevation, since we're thinking about DEMs here (though it's a general algorithm)
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"))

dem = raster(meuse.grid, layer = "dist") # if we want to use raster for large DEMs

create_perp_buffers <- function(x1, y1, x2, y2, grid_dist = 50, slice_width = 50, proj4string, id_prefix = "") {
  # creates buffers perpendicular to a line defined by x1, y1 and x2, y2.  Grid_dist is how far apart to take your swaths along
  #   the line, and slice_width is the width of the swath, all in terms of your projection.  Probably easiest to use meter-based 
  #   projections, and not degrees, etc.  Make sure grid_dist is as big or bigger than your raster grid size.  This is the really
  #   reusable part of this code.
  xdiff <- x2 - x1
  ydiff <- y2 - y1
  lineangle <- atan2(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)
  dist = as.matrix(dist(midpoints), diag=T, upper=T)[,1]
  
  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"=paste0(id_prefix,1:nrow(begin.coord)), begin.coord, end.coord, dist))
  proj4string(sldf) = proj4string
  
  blpi <- gBuffer(sldf, byid=TRUE, id=sldf$lineID, width = grid_dist/2)
  
  return(blpi)
}

# Make fake site data though which we place lines
site_data = data.frame(x = c(179000, 180000, 181000), y = c(330000, 331000, 333000), site = c("site1", "site2", "site3"))
coordinates(site_data) = c("x","y")
proj4string(site_data) = CRS(proj4string(meuse.grid))

# loop through our transects
ele_trans = list()
for (i in 1:(nrow(site_data)-1)) {
  this_site = site_data[i,]
  next_site = site_data[(i+1),]
  buff_trans <- create_perp_buffers(this_site@coords[1], this_site@coords[2], next_site@coords[1], next_site@coords[2], proj4string = proj4string(site_data), id_prefix = as.character(this_site$site))
  dem_rastlyr = rasterize(buff_trans, dem)
  buff_trans[["mean_elevation"]] = zonal(dem, dem_rastlyr, fun = "mean")[,2]
  buff_trans[["min_elevation"]] = zonal(dem, dem_rastlyr, fun = "min")[,2]
  buff_trans[["max_elevation"]] = zonal(dem, dem_rastlyr, fun = "max")[,2]
  buff_trans[["sd_elevation"]] = zonal(dem, dem_rastlyr, fun = "sd")[,2]
  ele_trans[[as.character(this_site$site)]] = buff_trans
  #spplot(buff_trans,"mean_elevation")
  buff_trans$site = NA
  buff_trans@data[1,]$site = as.character(this_site$site)  
  if (i == 1) {
    ele_trans_total = buff_trans
  } else {
    buff_trans$dist = buff_trans$dist + max(ele_trans_total$dist, na.rm = T)
    ele_trans_total = rbind(ele_trans_total, buff_trans)
  }
}

# add final site to last point
  ele_trans_total@data[nrow(ele_trans_total),"site"] = as.character(site_data[nrow(site_data),]$site)

# Take a look at the DEM and the polygons
spplot(ele_trans[[i]], "mean_elevation", sp.layout = ele_trans, xlim = bbox(site_data)[1,], ylim = bbox(site_data)[2,])

image(dem)
plot(ele_trans_total, col = "blue", add = T)

# Smooth elevation data and plot it
ele_profile = ele_trans_total@data

ele_profile = rename(ele_profile, c("dist" = "M_accumulated", "mean_elevation" = "Z", "site" = "plot_code"))

ele_smooth <- function(x, moving_window = 5) {
  # function to smooth a line
  x_smooth <- rollmean(x, k = moving_window, na.pad=T)
  
  # fill the beginning
  for (i in 1:round(moving_window/2)) {
    x_smooth[i] = rollmean(x[1:i], k = i)
  }
  x_smooth[1] = x[1]
  
  # fill the end
  for (i in (moving_window - round(moving_window/2)):1) {
    this_row = length(x) - i
    x_smooth[this_row] = rollmean(x[this_row:length(x)], k = i)[2]
  }
  x_smooth[length(x)] = x[length(x)]
  return(x_smooth)
}

# Smooth out the line
ele_profile$Z_smooth = ele_smooth(ele_profile$Z)

ele_profile$sd_elevation_smooth = ele_smooth(ele_profile$sd_elevation)

ggplot(ele_profile, aes(x=M_accumulated/1000, y=Z_smooth)) + 
  geom_text(aes(label = plot_code), hjust = 0, nudge_x = 0.5, size = 2.5) +
  geom_ribbon(aes(ymax = Z_smooth + 2*sd_elevation_smooth, ymin = Z_smooth - 2*sd_elevation_smooth), fill="gray52", alpha = 0.5) +
  geom_line(aes(y=Z), color = "grey45", size = 0.25) +
  xlab("Distance (km)") + ylab("Elevation (km)") + 
  geom_line(size = .5) +
  geom_point(data=subset(ele_profile, !is.na(plot_code)), aes(x=M_accumulated/1000, y=Z_smooth), fill = "black", size = 1.5, stroke = 0.5, color = "black", shape = 21) +
  theme_bw() +
  theme(panel.grid.major = element_line(linetype = 'dashed'), 
        panel.grid.minor = element_blank(), 
        panel.background = element_rect(colour = 'white'))

# write.csv(ele_profile, file = "elevation_data.csv") # if we want to save our data