I have a shapefile of US counties and high-resolution elevation data that spans the entire contiguous United States. My goal is to calculate a terrain ruggedness index for each county. The functions (that I've been able to find, e.g.
spatialEco::tri) all take raster layers as arguments.
Based on mdsummer's excellent answer and given a boundary shapefile and a raster layer of elevation data, it's easy to calculate zonal statistics:
require(sf) require(tidyverse) # Shapefile of US counties in California calif <- USAboundaries::us_counties("1960-01-01", resolution = "high", states = c("CA")) %>% mutate(county_fips = as.numeric(fips)) %>% select(county_fips, geometry) # Load elevation data (at a low resolution for now) elev <- elevatr::get_elev_raster(as(calif, "Spatial"), z = 2, src = "aws") # Group the elevation raster according to county_fips polymap <- fasterize::fasterize(calif, elev, field = "county_fips") elev[is.na(values(polymap))] <- NA # Zonal statistics # v <- raster::values zonal_stats <- tibble(value = raster::values(elev), county_fips = raster::values(polymap)) %>% group_by(county_fips) %>% summarize(mean_elev = mean(value)) map <- left_join(x = calif, y = zonal_stats, by = "county_fips") plot(map["mean_elev"])
I'm having difficulty seeing how to apply a function that takes a raster layer to each county individually. If I run the following code:
# Terrain Ruggedness Index (entire state) tri.calif <- spatialEco::tri(polymap) plot(tri.calif) tri.calif.crop <- crop(tri.calif, extent(calif)) plot(tri.calif.crop) plot(st_geometry(calif), add = TRUE)
this calculates the TRI across the state using the default cell size of the
but obviously these calculations aren't happening strictly within each county. How do I apply a function (like
tri) that takes a raster layer to the raster that's contained within each county individually?
Once I have that, it's easy enough to calculate the mean TRI across all cells within the county, for example, using the same zonal statistics approach described above?
Once I have that