title: How to perform raster calculation (e.g. aspect) on subset of raster based on point intersection in R
body: <p>I'm working with some raster data in R using the <code>raster</code> package. I want to calculate and extract some geographic information (e.g., slope, aspect) from the raster, but only at specific points (I also have some data as a <code>SpatialPointsDataFrame</code> at which I want to calculate slope/aspect/etc.). I'm doing this for several high-resolution rasters, and it seems like a poor use of resources to calculate this for every raster cell when I only need maybe 5-10% of them. </p>
  
  <p>I thought maybe the <code>raster::stackApply</code> function might work,  but that seems to perform calculations on subsets of a rasterBrick rather than calculations on subsets of a single raster based on point locations (please correct me if I'm wrong). I also thought I could do a for loop, where I extract the surrounding cells nearest each point of interest, and iteratively calculate slope/aspect that way. That seems clunky, and I was hoping for a more elegant or built-in solution, but it should work. </p>
  
  <p>These are my thoughts so far on the for loop, but I'm not sure how best to even do this. </p>
  
  <pre class="lang-r prettyprint-override"><code># Attach packages
  library(rgdal)
  library(raster)
  
  # Generate example raster data
  r = raster()
  set.seed(0)
  values(r) = runif(ncell(r), min = 0, max = 1000)
  
  # Generate example point data
  df.sp = SpatialPoints(
    coords = cbind(runif(25, min = -100, max = 100),  
                   runif(25, min = -50, max = 50)), 
    proj4string = crs(r))
  
  # Iterate on each row of SpatialPoints
  for (i in 1:nrow(df.sp)) {
  
    # Find cell index of current SpatialPoint
    cell.idx = raster::extract(r, df.sp[i,], cellnumbers = TRUE)[1]
  
    # Find indices of cells surrounding point of interest
    neighbors.idx = raster::adjacent(r, cell.idx, directions = 16)
  
    # Get DEM values for cell and surrounding cells
    vals.local = r[c(cell.idx, neighbors.idx[,2])]
  
    # Somehow convert this back to an appropriate georeferenced matrix
    #r.local = ...
  
    # Perform geometric calculations on local raster
    #r.stack = terrain(r.local, opt = c('slope', 'aspect'))
  
    # Remaining data extraction, etc. (I can take it from here...)
  }
  
  </code></pre>
  
  <p>In summary, I need a method to calculate slope and aspect from a DEM raster only at specific points as given by a SpatialPoints object. If you know of a pre-built or more elegant solution, great! If not, some help finishing the for loop (how to best extract a neighborhood of surrounding cells and run calculations on that) would be most appreciated as well.</p>