4

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 tri function:

TRI for California

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

  • 1
    One question, why? It does not make sense to calculate the central tendency of TRI for an entire county. You could get at county-level topographic variability by just deriving the zonal variance. The TRI is intended to represent local deviation and is not a metric that is suitable for aggregation to a large unit boundary (this would meet the definition of an ecological fallacy). It functionally ceases to represent a process when using the mean TRI across a large area. Please, take a moment to reconsider your question and associated analysis. – Jeffrey Evans Feb 8 at 21:12
  • @JeffreyEvans Taking the mean TRI of all the cells in an area is a standard technique in the economics literature (for example, this recent paper from a top economics journal). That's why I had this idea in the first place. Could you give a few more details as to why this represents an ecological fallacy? (That's a sincere question; GIS isn't my forte) – Michael A Feb 8 at 21:33
  • (Because, to be clear, I'm 100% open to the possibility that the top journals in my field are barking up the wrong tree. It's happened before) – Michael A Feb 8 at 21:36
  • 2
    Just because something is published does not mean that it is correct. In this case aggregating the TRI to a county is certainly incorrect. The distributional qualities of the metric, in relation to inference, become meaningless. Given the linked journal, bad dogs! I would highly recommend reading up on MAUP, perhaps starting with Cressie's "Change of support and the modifiable areal unit problem" and ecological fallacy in spatial data by reading Wakefield's "Spatial Aggregation and the Ecological Fallacy". – Jeffrey Evans Feb 8 at 23:13
2

Just because something is published does not mean that it is necessarly correct. In this case aggregating the TRI to a county is certainly incorrect. The distributional qualities of the metric, in relation to inference, become meaningless. Given the linked journal, bad dogs! You are functionally taking the mean of a derivative metric that represents localized mean deviation.

I would highly recommend reading up on MAUP, perhaps starting with Cressie's "Change of support and the modifiable areal unit problem" and ecological fallacy in spatial data by reading Wakefield's "Spatial Aggregation and the Ecological Fallacy".

Since the basic idea here is to identify topographic variability within an experimental unit to indicate "ruggedness", one could address the underlying distributions directly. Since highly relieved areas would also be expected to exhibit highly skewed, standard Gaussian moments may not be adequate. You can step out into non parametric statistics such as Median Absolute Deviation from Median (MAD).

Here is an example of what I am getting at and some potential solutions.

Add libraries and data.

library(raster)
library(spatialEco)
library(elevatr)
library(USAboundaries)

counties <- as(us_counties(map_date = "1930-01-01", 
              resolution = "high", states = c("CA")),
              "Spatial")
elev <- get_elev_raster(counties, z=5)

First, let's calculate the pixel-level TRI and calculate the mean for each county. You can see that the variability is not correctly represented, at least not visually.

r.tri <- spatialEco::tri(elev) 
counties@data <- data.frame(counties@data, r.tri = extract(r.tri, 
                            counties, fun=mean))
spplot(counties, "r.tri") 

Now we can calculate the MAD by passing the function directly to raster::extract.

counties@data <- data.frame(counties@data, tri = extract(elev, 
                            counties, fun=tri))
spplot(counties, "rough")

We can also write a global approximation of TRI using the median and the deviation value. This actually looks fairly reasonable and is comparable to MAD. Although, it did pick up Frenso county as very high ruggedness (which spans the southern Sierra's) whereas MAD did not.

tri <- function(x, ...) {
  x <- x[!is.na(x)]
  return( sqrt(sum(((median(x) - x)^2))) )
}

counties@data <- data.frame(counties@data, tri = extract(elev, 
                            counties, fun=tri))
spplot(counties, "tri")
1

Edit and disclaimer: Please see Jeffrey Evan's comment above. Below you will find a computational solution to OP's question, but users should be aware of the consequences of deciding to aggregate something like TRI across broad spatial scales.

This isn't super efficient, but it'll get the job done. If x is your raster and y is your county shape:

CountyTRI = function(mycounty, myraster){
    CoRas = raster::mask(x = myraster, mask = mycounty)
    CoTRI = raster::terrain(CoRas, opt = "TRI")
    CoMean = mean(as.matrix(CoTRI), na.rm=T)
    return(CoMean)
}

Edit: Here's a way to quickly apply the function to a state-level shapefile that has county information. Note that I switched the order of the arguments in the original CountyTRI function so that it can easilty be used in lapply. You'll need to replace NAME with whatever title the name column in calif has.

mycounties = unique(calif$NAME)

countyshapes = lapply(mycounties,
                      FUN = function(x){
                        return(calif[calif$NAME==x,])
                      })
countyshapes

countyTRIs = lapply(countyshapes,
                    FUN = CountyTRI,
                    myraster = elev)
  • Thank you! Two comments, though. I think fun = "mean" should be fun = mean because the raster::calc function actually takes a function as the second argument. Also, what exactly does this function return? It looks like it returns a RasterLayer, but when I convert it to a tibble (using the same fasterize function in my original code) it has many observations per county_fips, whereas I expected it to just have one, since that's the end goal. – Michael A Feb 8 at 17:34
  • (This is probably due to my unfamiliarity with the raster package; it's not meant as a criticism) – Michael A Feb 8 at 17:35
  • 1
    @MichaelA yep, you're right. See updated code so that CountyTRI returns a single value. That value is the mean value of all TRI grid cells within a given county. – JepsonNomad Feb 8 at 19:16
  • When I run CountyTRI(elev, calif), I get a single, scalar value. Shouldn't there be one value for each county? – Michael A Feb 8 at 20:16
  • Running as.matrix on the original result (from your first CountryTRI function) returns a 1329 by 1034 matrix, but there are only 58 counties in California, according to the calif shapefile. I don't know what to make of that object either. – Michael A Feb 8 at 20:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.