Use terra to calculate relative similarity of raster values between areas inside and outside of a group of polygons in R

Question

This question builds on a helpful solution provided for calculating uniqueness across a SpatRaster -- Using terra in R to calculate & map similarity and uniqueness across cells of a large GDM-based raster

Here, I want to calculate the similarity of each grid cell (across a study extent) to those found in a set of polygons, while accounting for uniqueness (uniqueness = average similarity of a given grid cell in relation to the entire study extent).

I provide an example and a possible solution. My solution runs without error.

• Looking for feedback on whether my code logic is sound.
• Also, why the routine only works with dissimilarity values (not similarity). Code for both is shown.

---

library(terra)

#### Assemble Sample Data ####

# create multilayer sample raster

set.seed(1)
r <- rast(ncol=100, nrow=100, xmin=-150, xmax=-80, ymin=20, ymax=60,
          nlyr=5, vals=runif(10000*5))
r[[1:3]][10:20] <- NA
r[100:120] <- NA

# create a polygon
lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(id=1, part=1, lon, lat)
pols1 <- vect(lonlat, type="polygons", crs="+proj=longlat +datum=WGS84")

# create a second polygon
lon <- c(-140.8, -135.2, -134.9, -138.9, -137.2, -138.4, -140.0)
lat <- c(51.3, 52.9, 52.4, 49.8, 47.6, 48.3, 47.6)
lonlat <- cbind(id=1, part=1, lon, lat)
pols2 <- vect(lonlat, type="polygons", crs="+proj=longlat +datum=WGS84")

# combine polygons into a single SpatVector
pa <- vect(c(pols1, pols2))

---

#### ANALYSES ####

# Use an arbitrary value for intercept to set up function
gdmRastMod <- list(intercept = 2)

#### Assemble dissimilarity data for entire study extent ####

# Calculate mean value of each layer in raster
m <- global(r, "mean", na.rm=T)

# Calculate sum of differences b/w raster grid cell values and mean raster grid cell values  
edist <- sum(r - unlist(m))

# Create SpatRaster of dissimilarity for entire study extent;
dissim <- 1 - exp(-1 * (gdmRastMod$intercept + edist))

# Alternative -- SpatRaster of similarity
# sim <- exp(-1 * (mod + edist))
# Note -- using similarity doesn't work; no errors, just nonsensical results 

#### Assemble dissimilarity data for areas in polygons (e.g., could be protected areas - PAs) ####

# Crop raster to PA extent
pa_rast <- crop(r, pa, mask = TRUE)

# calculate mean values in cropped out area of raster
m_pa <- global(pa_rast, "mean", na.rm=T)

# Calculate sum of differences b/w raster grid cell values and mean raster grid cell values in PAs
edist_pa <- sum(r - unlist(m_pa))

# create SpatRaster of dissimilarity for PAs 
dissim_pa <- 1 - exp(-1 * (gdmRastMod$intercept + edist_pa))

#### Calculate relative dissimilarity between PAs & study extent ####
rel_sim <- dissim_pa / dissim