roptimizationr-raster

How to perform raster calculation (e.g. aspect) on subset of raster based on point intersection in R


I'm working with some raster data in R using the raster 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 SpatialPointsDataFrame 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.

I thought maybe the raster::stackApply 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.

These are my thoughts so far on the for loop, but I'm not sure how best to even do this.

# 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...)
}

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.


Solution

  • Interesting question. Here is a possible approach.

    library(raster)
    
    r <- raster()
    set.seed(0)
    values(r) <- runif(ncell(r), min = 0, max = 1000)
    coords <- cbind(runif(25, min = -100, max = 100),  
                   runif(25, min = -50, max = 50))
    
    x <- rasterize(coords, r)
    f <- focal(x, w=matrix(1, nc=3, nr=3), na.rm=TRUE)             
    rr <- mask(r, f)
    slope <- terrain(rr, "slope")
    
    extract(slope, coords)
    # [1] 0.0019366236 0.0020670699 0.0006305257 0.0025334280 0.0023480935 0.0007527267 0.0002699272 0.0004699626
    # [9] 0.0004869054 0.0025651333 0.0010415805 0.0008574920 0.0010664869 0.0017700297 0.0001666226 0.0008405391
    #[17] 0.0017682167 0.0009854172 0.0015350466 0.0017714466 0.0012994945 0.0016563132 0.0003276584 0.0020499529
    #[25] 0.0006582073
    

    Probably not much efficiency gain, as it still processes all the NA values

    So maybe like this, more along your line of thinking:

    cells <- cellFromXY(r, coords)
    ngbs <- raster::adjacent(r, cells, pairs=TRUE)
    
    slope <- rep(NA, length(cells))
    for (i in 1:length(cells)) {
       ci <- ngbs[ngbs[,1] == cells[i], 2]
       e <- extentFromCells(r, ci)
       x <- crop(r, e)
       slope[i] <- terrain(x, "slope")[5]
    }
    slope
    #[1] 0.0019366236 0.0020670699 0.0006305257 0.0025334280 0.0023480935 0.0007527267 0.0002699272 0.0004699626
    #[9] 0.0004869054 0.0025651333 0.0010415805 0.0008574920 0.0010664869 0.0017700297 0.0001666226 0.0008405391
    #[17] 0.0017682167 0.0009854172 0.0015350466 0.0017714466 0.0012994945 0.0016563132 0.0003276584 0.0020499529
    #[25] 0.0006582073
    

    But I find that brute force

    slope <- terrain(r, "slope")
    extract(slope, coords)
    

    is fastest, 9x faster than my first alternative and 4 times faster than the second alternative