0

I have raster layer of days in which the temperature exceeds the threshold of 0 degrees, and a shapefile of the states of United States. I would like to extract the frequency of the pixel values per ecoregion in terra in R. So far I have tried this:

library(terra)

x = rast("D:/temperature/data.tif")
states = vect("D:/studyarea/US_states.shp")
states = project(states, x)
z = extract(x, states, na.rm = T, fun = freq(x))

and

z = extract(freq(x), states, na.rm = T)

But it did not work. How can I extract the frequency of the raster values per state?

1
  • What does "did not work" mean? Did it all come back zeroes? Error messages? Are your data correctly spatially referenced? Can you show a map or summary info of your data? We don't have your data so we can't replicate your problem.
    – Spacedman
    Commented Feb 13, 2023 at 11:12

3 Answers 3

2

Example data

library(terra)
r <- round(rast(system.file("ex/elev.tif", package="terra"))/100)
names(r) <- "days"
v <- vect(system.file("ex/lux.shp", package="terra"))

You can use crosstab for this after rasterizing the polygons

x <- rasterize(v, r, "NAME_2", wopt=list(names="states"))
ct <- crosstab(c(x, r))
ct
#                  days
#states               1   2   3   4   5
#  Capellen           0   0 260  70   0
#  Clervaux           0   0   3 139 419
#  Diekirch           0  43 195 129  27
#  Echternach         0  35 208  81   0
#  Esch-sur-Alzette   0   8 369  57   0
#  Grevenmacher       3  79 273  24   0
#  Luxembourg         0  14 321  88   0
#  Mersch             0  44 265 111   0
#  Redange            0   0 198 145 123
#  Remich             7 116  95   3   0
#  Vianden            0   7  46  45  32
#  Wiltz              0   0  48 279 146

You can create a "long" format data.frame like this

b <- as.data.frame(a)
head(b)
#            states days Freq
#1         Capellen    1    0
#2         Clervaux    1    0
#3         Diekirch    1    0
#4       Echternach    1    0
#5 Esch-sur-Alzette    1    0
#6     Grevenmacher    1    3

Alternatively, you can use extract

e <- extract(r, v, table, ID=FALSE)[[1]]

e is a list in which element matches a geometry in v. The table for the first geometry is

e[[1]]

  3   4   5 
  3 139 419 

But the nice thing about aldo_tapia's answer is that you can more easily combine the results for each layer (by calling do.call(rbind, output) after the loop. With the above you can get there like this:

z <- lapply(1:length(e), \(i) cbind(geom=i, as.data.frame(e[[i]])))
out <- do.call(rbind, z)
out$geom <- v$NAME_2[out$geom]

head(out)
#      geom Var1 Freq
#1 Clervaux    3    3
#2 Clervaux    4  139
#3 Clervaux    5  419
#4 Diekirch    2   43
#5 Diekirch    3  195
#6 Diekirch    4  129
1

You can use a for loop for this purpose, saving the results in a list:

output <- list()

for(i in seq_along(states)){
  x_temp <- crop(x, states[i,])
  output[[i]] <- freq(x_temp)
}
1
  • 3
    For this to be correct you need to also mask x_temp with states[i,]. You can use crop(x, states[i,], mask=TRUE) Commented Feb 13, 2023 at 2:00
0

This can also be done with exactextractr. Although exactextractr does not have a freq operation, we can get it by multiplying the results of frac and count.

Using the example data from Robert Hijmans' solution:

library(terra)
library(sf)
library(exactextractr)
library(dplyr)

r <- round(rast(system.file("ex/elev.tif", package="terra"))/100)
names(r) <- "days"
v <- st_read(system.file("ex/lux.shp", package="terra"))

exact_extract(r, st_as_sf(v), c('count', 'frac'), append_cols = 'NAME_2') %>%
    mutate(across(starts_with('frac'), function(x) x * count))  %>%
    rename_with(function(n) sub('frac', 'freq', n))

Result:

             NAME_2    count   freq_1     freq_2    freq_3    freq_4    freq_5
1          Clervaux 553.2812 0.000000   0.000000   2.58476 137.85263 412.84379
2          Diekirch 392.1883 0.000000  42.314376 192.78626 128.33565  28.75196
3           Redange 463.6169 0.000000   0.000000 198.47038 143.40036 121.74611
4           Vianden 129.0564 0.000000   6.783227  46.34947  45.16423  30.75945
5             Wiltz 472.7357 0.000000   0.000000  48.08434 276.51430 148.13704
6        Echternach 327.6882 0.000000  35.998649 210.01369  81.67591   0.00000
7            Remich 218.1701 6.276011 114.403333  94.49070   3.00000   0.00000
8      Grevenmacher 373.0518 1.997318  77.842841 270.82947  22.38219   0.00000
9          Capellen 330.4629 0.000000   0.000000 257.97941  72.48348   0.00000
10 Esch-sur-Alzette 432.5632 0.000000   7.804500 367.45097  57.30769   0.00000
11       Luxembourg 424.8329 0.000000  13.882485 321.89943  89.05096   0.00000
12           Mersch 419.2414 0.000000  44.480128 265.39091 109.37033   0.00000

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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