Search code examples
rgeospatialr-rasterterraintopography

Calculate Topographic Position Index for coordinates stored in a dataframe in R


I'm trying to calculate Topographic Position Index (TPI) for 177 points of interest. I have their coordinates stored in a data.frame and elevation in a raster of 7.5 arc sec spatial resolution. And the TPIs I'm calculating is basically: the elevation of point of interest minus the average elevation of its surrounding cells, then the intermediate result is divided by the spatial resolution of the raster. (resolution(dem)) to account for differences in the spatial scale of the DEM and the TPI values.

And since studies usually calculate two TPIs (a small scale + a large scale), I'm also using two windows, where in the small one the surrounding 55 cells are used, and in the large one the surrounding 1010 cells are used.

I am getting this error message

Error in .focal_fun(v, w, as.integer(c(tr$nrows[1] + addr, nc)), runfun,  :    
   Evaluation error: could not find function "resolution"

Code:

library(raster)
library(sp)
library(terra)
library(haven)

# Make a dataframe with longitudes and latitudes
df <- data.frame(lon = coords$longitude, lat = coords$latitude)

# Convert the dataframe to a SpatialPointsDataFrame
coordinates(df) <- c("lon", "lat")
proj4string(df) <- CRS("+proj=longlat +datum=WGS84")

# Extract the elevation values at the points
elev <- extract(dem, df)

# Define the scales for TPI calculation
scales <- c(5, 10)

# Loop over the scales and calculate TPI
tpi_list <- list()
for (scale in scales) {
  # Define the size of the moving window
  win_size <- scale * 5
  
  # Calculate TPI
  tpi <- focal(dem, w = matrix(1, win_size, win_size), fun = function(x) {
    (elev - mean(x)) / resolution(dem) * 5
  })

  # Extract TPI values at the points
  tpi_vals <- extract(tpi, df)
  
  # Store the TPI values in a list
  tpi_list[[as.character(scale)]] <- tpi_vals
}

#Error in .focal_fun(v, w, as.integer(c(tr$nrows[1] + addr, nc)), runfun,  :    Evaluation error: could not find function "resolution"

# Combine the TPI values for different scales into a dataframe
tpi_df <- data.frame(tpi_list, row.names = rownames(df))

Solution

  • The error you get is clear:

    Evaluation error: could not find function "resolution"
    

    You are using a function resolution, but R does not know about that function. It does not exist in the current workspace. I suppose you were looking for res.

    Here is a working example.

    library(terra)
    dem <- rast(system.file("ex/elev.tif", package="terra"))
    df <- data.frame(lon=c(5.9, 6.0, 6.2), lat=c(49.9, 49.6, 49.7))
    
    tpifun <- \(x, f) x[f] - mean(x[-f], na.rm=TRUE)
    
    scales <- c(5, 11)
    tpilst <- vector("list", length(scales))
    for (i in seq_along(scales)) {
        win_size <- scales[i] * 5
        mid <- ceiling(win_size^2 / 2)
        tpi <- focal(dem, w=win_size, fun=tpifun, f=mid, wopt=list(names="tpi"))
        tpilst[[i]] <- data.frame(scale=scales[i], extract(tpi, df))
    }
    
    tpi <- do.call(rbind, tpilst)
    tpi$tpi <- tpi$tpi / (mean(res(dem)) * 5)
    
    tpi
    #  scale ID        tpi
    #1     5  1 -1330.6184
    #2     5  2   340.1538
    #3     5  3  -135.3077
    #4    11  1  -585.7952
    #5    11  2   255.3344
    #6    11  3   292.0155
    

    A couple of things:

    res returns two numbers, the x and y resolution. In the example above I take the mean.

    What you were doing in the function supplied to focal is not possible. You supplied a data.frame with elevation data for a few points. How can focal understand what that is all about? Instead, you can compute the TPI for each cell and extract these values.

    You cannot use a value of 10 for scale because the weights matrix must have odd size. Otherwise it is not clear how it should be centered on the focal cell.

    You say that if scale is 5, the "surrounding 55 cells are used". But that is not the case. The number of surrounding cells used is 624.

    scale <- 5
    (scale * 5)^2 - 1
    #[1] 624