05_SpatialInteractionVectorRaster.Rmd

# Spatial Interactions of Vector and Raster Data {#int-RV}

```{r chap5_setup, echo = FALSE, results = "hide"}
library(knitr)
knitr::opts_chunk$set(
  echo = TRUE,
  cache = FALSE,
  # cache = TRUE,
  comment = NA,
  message = FALSE,
  warning = FALSE,
  tidy = FALSE,
  cache.lazy = FALSE
)

suppressMessages(library(here))
opts_knit$set(root.dir = here())
```

```{r , eval = FALSE, echo = FALSE}
setwd(here())
```

```{r, include = FALSE, cache = FALSE}
#--- load packages ---#
library(data.table)
library(exactextractr)
library(prism)
library(sf)
library(terra)
library(stars)
library(raster)
library(tidyverse)
library(DT)
library(tictoc)
library(tmap)
library(parallel)
library(maps)
```

```{r figure_setup, echo = FALSE}
theme_update(
  axis.title.x = element_text(size = 12, angle = 0, hjust = .5, vjust = -0.3, face = "plain", family = "Times"),
  axis.title.y = element_text(size = 12, angle = 90, hjust = .5, vjust = .9, face = "plain", family = "Times"),
  axis.text.x = element_text(size = 10, angle = 0, hjust = .5, vjust = 1.5, face = "plain", family = "Times"),
  axis.text.y = element_text(size = 10, angle = 0, hjust = 1, vjust = 0, face = "plain", family = "Times"),
  axis.ticks = element_line(size = 0.3, linetype = "solid"),
  # axis.ticks = element_blank(),
  axis.ticks.length = unit(.15, "cm"),
  # axis.ticks.margin = unit(.1,'cm'),
  # axis.text = element_text(margin=unit(.1,'cm')),

  #--- legend ---#
  legend.text = element_text(size = 10, angle = 0, hjust = 0, vjust = 0, face = "plain", family = "Times"),
  legend.title = element_text(size = 10, angle = 0, hjust = 0, vjust = 0, face = "plain", family = "Times"),
  legend.key.size = unit(0.5, "cm"),

  #--- strip (for faceting) ---#
  strip.text = element_text(size = 10, family = "Times"),

  #--- plot title ---#
  plot.title = element_text(family = "Times", face = "bold", size = 12),

  #--- margin ---#
  # plot.margin = margin(0, 0, 0, 0, "cm"),

  #--- panel ---#
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  panel.border = element_rect(fill = NA)
)
```

```{r ch05-define_theme, echo = F}
theme_set(theme_bw())

theme_for_map <- theme(
  axis.ticks = element_blank(),
  axis.text = element_blank(),
  axis.line = element_blank(),
  panel.border = element_blank(),
  panel.grid.major = element_line(color = "transparent"),
  panel.grid.minor = element_line(color = "transparent"),
  panel.background = element_blank(),
  plot.background = element_rect(fill = "transparent", color = "transparent")
)
```

## Before you start {-}

In this chapter we learn the spatial interactions of a vector and raster dataset. We first look at how to crop (spatially subset) a raster dataset based on the geographic extent of a vector dataset. We then cover how to extract values from raster data for points and polygons. To be precise, here is what we mean by raster data extraction and what it does for points and polygons data:

+ **Points**: For each of the points, find which raster cell it is located within, and assign the value of the cell to the point.  
 
+ **Polygons**: For each of the polygons, identify all the raster cells that intersect with the polygon, and assign a vector of the cell values to the polygon

This is probably the most important operation economists run on raster datasets. 

We will show how we can use `terra::extract()` for both cases. But, we will also see that for polygons, `exact_extract()` from the `exactextractr` package is often considerably faster than `terra::extract()`.

Finally, you will see conversions between `Raster`$^*$ (`raster` package) objects and `SpatRaster` object (`terra` package) because of the incompatibility of object classes across the key packages. I believe that these hassles will go away soon when they start supporting each other.

### Direction for replication {-}

**Datasets**

All the datasets that you need to import are available [here](https://www.dropbox.com/sh/l84zfidaxmjrti9/AAB1GrDRoIlidJ3_zArMN24ua?dl=0). In this chapter, the path to files is set relative to my own working directory (which is hidden). To run the codes without having to mess with paths to the files, follow these steps:

+ set a folder (any folder) as the working directory using `setwd()`
+ create a folder called "Data" inside the folder designated as the working directory (if you have created a "Data" folder previously, skip this step)
+ download the pertinent datasets from [here](https://www.dropbox.com/sh/l84zfidaxmjrti9/AAB1GrDRoIlidJ3_zArMN24ua?dl=0) 
+ place all the files in the downloaded folder in the "Data" folder

**Packages**

Run the following code to install or load (if already installed) the `pacman` package, and then install or load (if already installed) the listed package inside the `pacman::p_load()` function.

```{r Chap5_packages, eval = F}
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  terra, # handle raster data
  raster, # handle raster data
  exactextractr, # fast extractions
  sf, # vector data operations
  dplyr, # data wrangling
  tidyr, # data wrangling
  data.table, # data wrangling
  prism, # download PRISM data
  tictoc, # timing codes
  tigris, # to get county sf
  tmap # for mapping
)
```

## Cropping to the Area of Interest {#raster-crop}

Here we use PRISM maximum temperature (tmax) data as a raster dataset and Kansas county boundaries as a vector dataset. 

Let's download the tmax data for July 1, 2018 (Figure \@ref(fig:prism-tmax-map)).

```{r prism_download, cache = F, results = "hide"}
#--- set the path to the folder to which you save the downloaded PRISM data ---#
# This code sets the current working directory as the designated folder
options(prism.path = "Data")

#--- download PRISM precipitation data ---#
get_prism_dailys(
  type = "tmax",
  date = "2018-07-01",
  keepZip = FALSE
)

#--- the file name of the PRISM data just downloaded ---#
prism_file <- "Data/PRISM_tmax_stable_4kmD2_20180701_bil/PRISM_tmax_stable_4kmD2_20180701_bil.bil"

#--- read in the prism data ---#
prism_tmax_0701_sr <- rast(prism_file)
```

```{r prism-tmax-map, echo = F, fig.cap = "Map of PRISM tmax data on July 1, 2018"}
tm_shape(as(prism_tmax_0701_sr, "Raster")) +
  tm_raster() +
  tm_layout(frame = NA)
```

We now get Kansas county border data from the `tigris` package (Figure \@ref(fig:ks-county-map)) as `sf`. 

```{r get_KS_county, cache = F, results = "hide"}
#--- Kansas boundary (sf) ---#
KS_county_sf <-
  #--- get Kansas county boundary ---
  tigris::counties(state = "Kansas", cb = TRUE) %>%
  #--- sp to sf ---#
  st_as_sf() %>%
  #--- transform using the CRS of the PRISM tmax data  ---#
  st_transform(terra::crs(prism_tmax_0701_sr))
```

```{r ks-county-map, echo = FALSE, fig.cap = "Kansas county boundaries"}
#--- gen map ---#
tm_shape(KS_county_sf) +
  tm_polygons() +
  tm_layout(frame = NA)
```

---

Sometimes, it is convenient to crop a raster layer to the specific area of interest so that you do not have to carry around unnecessary parts of the raster layer. Moreover, it takes less time to extract values from a raster layer when the size of the raster layer is smaller. You can crop a raster layer by using `terra::crop()`. It works like this:

```{r crop_syntax, eval = FALSE}
#--- syntax (this does not run) ---#
terra::crop(SpatRaster, sf)
```

So, in this case, this does the job.

```{r crop_prism_to_KS, cache = F}
#--- crop the entire PRISM to its KS portion---#
prism_tmax_0701_KS_sr <- terra::crop(prism_tmax_0701_sr, KS_county_sf)
```

The figure below (Figure \@ref(fig:prism-ks-viz)) shows the PRISM tmax raster data cropped to the geographic extent of Kansas. Notice that the cropped raster layer extends beyond the outer boundary of Kansas state boundary (it is a bit hard to see, but look at the upper right corner).  

```{r prism-ks-viz, fig.cap = "PRISM tmax raster data cropped to the geographic extent of Kansas", echo = F}
tm_shape(as(prism_tmax_0701_KS_sr, "Raster")) +
  tm_raster() +
  tm_shape(KS_county_sf) +
  tm_borders(col = "blue")
```

## Extracting Values from Raster Layers for Vector Data

In this section, we will learn how to extract information from raster layers for spatial units represented as vector data (points and polygons). For the demonstrations in this section, we use the following datasets:

+ Raster: PRISM tmax data cropped to Kansas state border for 07/01/2018 (obtained in \@ref(raster-crop)) and 07/02/2018 (downloaded below)
+ Polygons: Kansas county boundaries (obtained in \@ref(raster-crop))
+ Points: Irrigation wells in Kansas (imported below) 

### Simple visual illustrations of raster data extraction

**Extracting to Points**

Figure \@ref(fig:illustrate-points) shows visually what we mean by "extract raster values to points."

```{r illustrate-points, fig.cap = "Visual illustration of raster data extraction for points data", echo = F}
set.seed(378533)

#--- create polygons ---#
polygon <- st_polygon(list(
  rbind(c(0, 0), c(8, 0), c(8, 8), c(0, 8), c(0, 0))
))

raster_like_cells <- st_make_grid(polygon, n = c(8, 8)) %>%
  st_as_sf() %>%
  mutate(value = sample(1:64, 64))

stars_cells <- st_rasterize(raster_like_cells, nx = 8, ny = 8)

cell_centroids <- st_centroid(raster_like_cells) %>%
  st_as_sf()

#--------------------------
# Create points for which values are extracted
#--------------------------
#--- points ---#
point_1 <- st_point(c(2.4, 2.2))
point_2 <- st_point(c(6.7, 1.8))
point_3 <- st_point(c(4.2, 7.1))

#--- combine the points to make a single  sf of points ---#
points <- list(point_1, point_2, point_3) %>%
  st_sfc() %>%
  st_as_sf() %>%
  mutate(point_name = c("Point 1", "Point 2", "Point 3"))

#--------------------------
# Create maps
#--------------------------
ggplot() +
  geom_stars(data = stars_cells, alpha = 0.5) +
  scale_fill_distiller(name = "Value", palette = "Spectral") +
  geom_sf_text(data = raster_like_cells, aes(label = value)) +
  geom_sf(data = points, aes(shape = point_name), size = 2) +
  scale_shape(name = "Points") +
  theme_void() +
  theme_for_map
```

In the figure, we have grids (cells) and each grid holds a value (presented at the center). There are three points. We will find which grid the points fall inside and get the associated values and assign them to the points. In this example, Points 1, 2, and 3 will have 50, 4, 54, respectively,

**Extracting to Polygons**

Figure \@ref(fig:polygons-extact-viz) shows visually what we mean by "extract raster values to polygons."

```{r polygons-extact-viz, fig.cap = "Visual illustration of raster data extraction for polygons data", echo = F}
#--------------------------
# Create a polygon for which values are extracted
#--------------------------
polygon_extract <- st_polygon(list(
  rbind(c(1.5, 2), c(6, 2.3), c(7, 6.5), c(2, 5), c(1.5, 2))
))

polygons_extract_viz <-
  ggplot() +
  geom_stars(data = stars_cells, alpha = 0.5) +
  scale_fill_distiller(name = "Value", palette = "Spectral") +
  geom_sf(data = polygon_extract, fill = "gray", alpha = 0.5) +
  geom_sf(data = cell_centroids, color = "black", size = 0.8) +
  geom_sf_text(data = raster_like_cells, aes(label = value), nudge_x = -0.25, nudge_y = 0.25) +
  theme_void() +
  theme_for_map

polygons_extract_viz
```

There is a polygon overlaid on top of the cells along with their centroids represented by black dots. Extracting raster values to a polygon means finding raster cells that intersect with the polygons and get the value of all those cells and assigns them to the polygon. As you can see some cells are completely inside the polygon, while others are only partially overlapping with the polygon. Depending on the function you use and its options, you regard different cells as being spatially related to the polygons. For example, by default, `terra::extract()` will extract only the cells whose centroid is inside the polygon. But, you can add an option to include the cells that are only partially intersected with the polygon. In such a case, you can also get the fraction of the cell what is overlapped with the polygon, which enables us to find area-weighted values later. We will discuss the details these below.

**PRISM tmax data for 07/02/2018**

```{r download_07022018, cache = F, results = "hide"}
#--- download PRISM precipitation data ---#
get_prism_dailys(
  type = "tmax",
  date = "2018-07-02",
  keepZip = FALSE
)

#--- the file name of the PRISM data just downloaded ---#
prism_file <- "Data/PRISM_tmax_stable_4kmD2_20180702_bil/PRISM_tmax_stable_4kmD2_20180702_bil.bil"

#--- read in the prism data and crop it to Kansas state border ---#
prism_tmax_0702_KS_sr <- rast(prism_file) %>%
  terra::crop(KS_county_sf)
```

**Irrigation wells in Kansas:**

```{r import_KS_wells}
#--- read in the KS points data ---#
(
  KS_wells <- readRDS("Data/Chap_5_wells_KS.rds")
)
```

---

Here is how the wells are spatially distributed over the PRISM grids and Kansas county borders (Figure \@ref(fig:tmax-prism-wells)):

```{r tmax-prism-wells, fig.cap = "Map of Kansas county borders, irrigation wells, and PRISM tmax", echo = F}
tm_shape(KS_county_sf) +
  tm_polygons(alpha = 0) +
  tm_shape(KS_wells) +
  tm_symbols(size = 0.02) +
  tm_layout(
    frame = FALSE,
    legend.outside = TRUE,
    legend.outside.position = "bottom"
  )
```

### Extracting to Points

You can extract values from raster layers to points using `terra::extract()`. `terra::extract()` finds which raster cell each of the points is located within and assigns the value of the cell to the point. 

```{r eval = F}
#--- syntax (this does not run) ---#
terra::extract(raster, points)
```

Let's extract tmax values from the PRISM tmax layer (`prism_tmax_0701_KS_sr`) to the irrigation wells:

```{r error = TRUE}
tmax_from_prism <- terra::extract(prism_tmax_0701_KS_sr, KS_wells)
```

Oh, okay. So, the problem is that the points data (`KS_wells`) is an `sf` object, and `terra::extract()` does not like it (I thought this would be fixed by now ...). Let's turn KS_wells into an `SpatVect` object and then apply `terra::extract()`:

```{r extract_tmax}
tmax_from_prism <- terra::extract(prism_tmax_0701_KS_sr, vect(KS_wells))

#--- take a look ---#
head(tmax_from_prism)
```

The resulting object is a `data.frame`, where the `ID` variable represents the order of the observations in the points data and the second column represents that values extracted for the points from the raster cells. So, you can assign the extracted values as a new variable of the points data as follows: 

```{r}
KS_wells$tmax_07_01 <- tmax_from_prism[, -1]
```   

Extracting values from a multi-layer `SpatRaster` works the same way. Here, we combine `prism_tmax_0701_KS_sr` and `prism_tmax_0702_KS_sr` to create a multi-layer `SpatRaster` and then extract values from them.

```{r extract_tmax_run_stack}
#--- create a multi-layer SpatRaster ---#
prism_tmax_stack <- c(prism_tmax_0701_KS_sr, prism_tmax_0702_KS_sr)

#--- extract tmax values ---#
tmax_from_prism_stack <- terra::extract(prism_tmax_stack, vect(KS_wells))

#--- take a look ---#
head(tmax_from_prism_stack)
```

We now have two columns that hold values extracted from the raster cells: the 2nd column for the 1st raster layer and the 3rd column for the 2nd raster layer in `prism_tmax_stack`.

**Side note:**

Interestingly, `terra::extract()` can work with `sf` as long as the raster object is of type `Raster`$^*$. In the code below, `prism_tmax_stack` is converted to `RasterStack` before `terra::extract()` does its job and it works just fine.

```{r extract-raster-sf-code, eval = F}
tmax_from_prism <- terra::extract(as(prism_tmax_stack, "Raster"), KS_wells)
```

```{r check-extracted-values-points}
#--- take a look ---#
head(tmax_from_prism)

#--- check the class ---#
class(tmax_from_prism)
```

As you can see, the resulting object is a `matrix` instead of a `data.frame`. Now, let's see if we can use a `SpatVect` instead of an `sf`.

```{r extract-raster-sv, error = TRUE}
tmax_from_prism <-
  terra::extract(
    as(prism_tmax_stack, "Raster"),
    vect(KS_wells)
  )
```

So, it does not work. It turns out that `terra::extract()` works for the combinations of `SpatRaster`-`SpatVect` or `Raster`$^*$-`sf`/`sp`.

### Extracting to Polygons (`terra` way)

You can use `terra::extract()` for extracting raster cell values for polygons as well. For each of the polygons, it will identify all the raster cells whose center lies inside the polygon and assign the vector of values of the cells to the polygon. Let's first convert `KS_county_sf` (an `sf` object) to a `SpatVector`.

```{r convert-sv}
#--- Kansas boundary (SpatVector) ---#
KS_county_sv <- vect(KS_county_sf)
```

<!-- (You can change this to the cells that intersect with polygons using the `touch = TRUE` option).  -->

Now, let's extract tmax values for each of the KS counties.  

```{r terra-extract-polygon, eval = F}
#--- extract values from the raster for each county ---#
tmax_by_county <- terra::extract(prism_tmax_0701_KS_sr, KS_county_sv)
```

```{r terra-extract-polygon-show, echo = F}
# saveRDS(tmax_by_county, "Data/tmax-by-county-terra.rds")
tmax_by_county <- readRDS("Data/tmax-by-county-terra.rds")
```

```{r check-the-extracted-values}
#--- check the class ---#
class(tmax_by_county)

#--- take a look ---#
head(tmax_by_county)

#--- take a look ---#
tail(tmax_by_county)
```

So, `terra::extract()` returns a `data.frame`, where `ID` values represent the corresponding row number in the polygons data. For example, the observations with `ID == n` are for the **n**th polygon. Using this information, you can easily merge the extraction results to the polygons data. Suppose you are interested in the mean of the tmax values of the intersecting cells for each of the polygons, then you can do the following:

```{r c05-mean-merge}
#--- get mean tmax ---#
mean_tmax <-
  tmax_by_county %>%
  group_by(ID) %>%
  summarize(tmax = mean(PRISM_tmax_stable_4kmD2_20180701_bil))

(
  KS_county_sf <-
    #--- back to sf ---#
    st_as_sf(KS_county_sv) %>%
    #--- define ID ---#
    mutate(ID := seq_len(nrow(.))) %>%
    #--- merge by ID ---#
    left_join(., mean_tmax, by = "ID")
)
```

Instead of finding the mean after applying `terra::extract()` as done above, you can do that within `terra::extract()` using the `fun` option. 

```{r terra-extract-polygon-mean, eval = F}
#--- extract values from the raster for each county ---#
tmax_by_county <-
  terra::extract(
    prism_tmax_0701_KS_sr,
    KS_county_sv,
    fun = mean
  )
```

```{r terra-extract-polygon-mean-show, echo = F}
# saveRDS(tmax_by_county, "Data/tmax-by-county-mean.rds")
tmax_by_county <- readRDS("Data/tmax-by-county-mean.rds")
```

```{r take-a-look-e-polygons-sum}
#--- take a look ---#
head(tmax_by_county)
```

You can apply other summary functions like `min()`, `max()`, `sum()`.

---

Extracting values from a multi-layer raster data works exactly the same way except that data processing after the value extraction is slightly more complicated. 

```{r terra_exatrac_from_stack_run}
#--- extract from a multi-layer raster object ---#
tmax_by_county_from_stack <-
  terra::extract(
    prism_tmax_stack,
    KS_county_sv
  )

#--- take a look ---#
head(tmax_by_county_from_stack)
```

Similar to the single-layer case, the resulting object is a `data.frame` and you have two columns corresponding to each of the layers in the two-layer raster object. 

--- 

Sometimes, you would like to have how much of the raster cells are intersecting with the intersecting polygon to find area-weighted summary later. In that case, you can add `exact = TRUE` option to `terra::extract()`.

```{r extract-with-weights}
#--- extract from a multi-layer raster object ---#
tmax_by_county_from_stack <-
  terra::extract(
    prism_tmax_stack,
    KS_county_sv,
    exact = TRUE
  )

#--- take a look ---#
head(tmax_by_county_from_stack)
```

As you can see, now you have `fraction` column to the resulting `data.frame` and you can find area-weighted summary of the extracted values like this:

```{r area-weighted-extract}
tmax_by_county_from_stack %>%
  group_by(ID) %>%
  summarize(
    tmax_0701 = sum(fraction * PRISM_tmax_stable_4kmD2_20180701_bil) / sum(fraction),
    tmax_0702 = sum(fraction * PRISM_tmax_stable_4kmD2_20180702_bil) / sum(fraction)
  )
```

### Extracting to Polygons (`exactextractr` way)

`exact_extract()` function from the `exactextractr` package is a faster alternative than `terra::extract()` for a large raster data as we confirm later (`exact_extract()` does not work with points data at the moment).^[See [here](https://github.com/isciences/exactextract) for how it does extraction tasks differently from other major GIS software.] `exact_extract()` also provides a coverage fraction value for each of the cell-polygon intersections. The syntax of `exact_extract()` is very much similar to `terra::extract()`. 

```{r eval = FALSE}
#--- syntax (this does not run) ---#
exact_extract(raster, polygons)
```

`exact_extract()` can accept both `SpatRaster` and `Raster`$^*$ objects as the raster object. However, while it accepts `sf` as the polygons object, but it does not accept `SpatVect`. 

Let's get tmax values from the PRISM raster layer for Kansas county polygons, the following does the job: 

```{r exact_extract}
library("exactextractr")

#--- extract values from the raster for each county ---#
tmax_by_county <-
  exact_extract(
    prism_tmax_0701_KS_sr,
    KS_county_sf,
    #--- this is for not displaying progress bar ---#
    progress = FALSE
  )
```

The resulting object is a list of `data.frame`s where **i**th element of the list is for **i**th polygon of the polygons `sf` object.

```{r show_the_results}
#--- take a look at the first 6 rows of the first two list elements ---#
tmax_by_county[1:2] %>% lapply(function(x) head(x))
```

For each element of the list, you see `value` and `coverage_fraction`. `value` is the tmax value of the intersecting raster cells, and `coverage_fraction` is the fraction of the intersecting area relative to the full raster grid, which can help find coverage-weighted summary of the extracted values (like `fraction` variable when you use `terra::extract()` with `exact = TRUE`). 

We can take advantage of `dplyr::bind_rows()` to combine the list of the datasets into a single `data.frame`. In doing so, you can use `.id` option to create a new identifier column that links each row to its original data.


```{r combine_after_ee, cache = F}
(
  #--- combine ---#
  tmax_combined <- bind_rows(tmax_by_county, .id = "id") %>%
    as_tibble()
)
```

`data.table` users can use `rbindlist()` with the `idcol` option like below:

```{r rbindlist-dt}
rbindlist(tmax_by_county, idcol = "id")
```

Note that `id` value of **i** represents **i**th element of the list, which in turn corresponds to **i**th polygon in the polygons sf data. Let's summarize the data by `id` and then merge it back to the polygons `sf` data. Here, we calculate coverage-weighted mean of tmax.

```{r transform_after_ee, cache = F}
#--- weighted mean ---#
(
  tmax_by_id <- tmax_combined %>%
    #--- convert from character to numeric  ---#
    mutate(id = as.numeric(id)) %>%
    #--- group summary ---#
    group_by(id) %>%
    summarise(tmax_aw = sum(value * coverage_fraction) / sum(coverage_fraction))
)

#--- merge ---#
KS_county_sf %>%
  mutate(id := seq_len(nrow(.))) %>%
  left_join(., tmax_by_id, by = "id") %>%
  dplyr::select(id, tmax_aw)
```

---

Extracting values from `RasterStack` works in exactly the same manner as `RasterLayer`.

```{r exatrac_from_stack_run, eval = F}
tmax_by_county_stack <-
  exact_extract(
    prism_tmax_stack,
    KS_county_sf,
    progress = F
  )

#--- take a look at the first 6 lines of the first element---#
tmax_by_county_stack[[1]] %>% head()
```

```{r exatrac_from_stack_run-show, echo = F}
# saveRDS(tmax_by_county_stack, "Data/tmax_by_county_stack.rds")
tmax_by_county_stack <- readRDS("Data/tmax_by_county_stack.rds")

#--- take a look at the first 6 lines of the first element---#
tmax_by_county_stack[[1]] %>% head()
```

As you can see above, `exact_extract()` appends additional columns for additional layers.

```{r combine_them}
#--- combine them ---#
tmax_all_combined <- bind_rows(tmax_by_county_stack, .id = "id")

#--- take a look ---#
head(tmax_all_combined)
```

In order to find the coverage-weighted tmax by date, you can first pivot it to a long format using `dplyr::pivot_longer()`.

```{r pivot_to_longer}
#--- pivot to a longer format ---#
(
  tmax_long <- pivot_longer(
    tmax_all_combined,
    -c(id, coverage_fraction),
    names_to = "date",
    values_to = "tmax"
  )
)
```

And then find coverage-weighted tmax by date:

```{r dplyr_way_tmax_cov}
(
  tmax_long %>%
    group_by(id, date) %>%
    summarize(tmax = sum(tmax * coverage_fraction) / sum(coverage_fraction))
)
```

For `data.table` users, this does the same:

```{r datatable_way_tmax_cov}
(
  tmax_all_combined %>%
    data.table() %>%
    melt(id.var = c("id", "coverage_fraction")) %>%
    .[, .(tmax = sum(value * coverage_fraction) / sum(coverage_fraction)), by = .(id, variable)]
)
```

## Extraction speed comparison {#extract-speed}

Here we compare the extraction speed of `raster::extract()`, `terra::extract()`, and `exact_extract()`. 

### Points: `terra::extract()` and `raster::extract()`

`terra::extract()` uses C++ as the backend. Therefore, it is considerably faster than `raster::extract()`.

```{r points_extraction_comp}
#--- terra ---#
KS_wells_sv <- vect(KS_wells)
tic()
temp <- terra::extract(prism_tmax_0701_KS_sr, KS_wells_sv)
toc()

#--- raster ---#
prism_tmax_0701_KS_lr <- as(prism_tmax_0701_KS_sr, "Raster")
tic()
temp <- raster::extract(prism_tmax_0701_KS_lr, KS_wells)
toc()
```

As you can see, `terra::extract()` is much faster. The time differential between the two packages can be substantial as the raster data becomes larger. Until recently, `raster::extract()` was much much slower. However, recent updates to the package made significant improvement in extraction speed.

### Polygons: `exact_extract()`, `terra::extract()`, and `raster::extract()`

`terra::extract()` is faster than `exact_extract()` for a relatively small raster data. Let's time them and see the difference.  

```{r extract-from-polygons-comp, cahce = F}
#--- terra::extract ---#
tic()
terra_extract_temp <-
  terra::extract(
    prism_tmax_0701_KS_sr,
    KS_county_sv
  )
toc()

#--- exact_extract ---#
tic()
exact_extract_temp <-
  exact_extract(
    prism_tmax_0701_KS_sr,
    KS_county_sf,
    progress = FALSE
  )
toc()

#--- raster::extract ---#
tic()
raster_extract_temp <-
  raster::extract(
    prism_tmax_0701_KS_lr,
    KS_county_sf
  )
toc()
```

As you can see, `raster::extract()` is the slowest. `terra::extract()` is faster than `exact_extract()`. However, once the raster data becomes larger (or spatially finer), then `exact_extact()` starts to shine. 

---

Let's disaggregate the prism data by a factor of 10 to create a much larger raster data.^[We did not introduce this function as it is very rare that you need this function in research projects.]

```{r prism_disaggregate}
#--- disaggregate ---#
(
  prism_tmax_0701_KS_sr_10 <- terra::disagg(prism_tmax_0701_KS_sr, fact = 10)
)
```

The disaggregated PRISM data now has 10 times more rows and columns.   

```{r dimensions}
#--- original ---#
dim(prism_tmax_0701_KS_sr)

#--- disaggregated ---#
dim(prism_tmax_0701_KS_sr_10)
``` 

---

Now, let's compare `terra::extrct()` and `exact_extrct()` using the disaggregated data.

```{r extract_from_polygons_comp_2}
#--- terra extract ---#
tic()
terra_extract_temp <-
  terra::extract(
    prism_tmax_0701_KS_sr_10,
    KS_county_sv
  )
toc()

#--- exact extract ---#
tic()
exact_extract_temp <-
  exact_extract(
    prism_tmax_0701_KS_sr_10,
    KS_county_sf,
    progress = FALSE
  )
toc()
```

As you can see, `exact_extract()` is considerably faster. The difference in time becomes even more pronounced as the size of the raster data becomes larger and the number of polygons are greater. The time difference of several seconds seem nothing, but imagine processing PRISM files for the entire US over 20 years, then you would appreciate the speed of `exact_extract()`. 

### Single-layer vs multi-layer

Pretend that you have five dates of PRISM tmax data (here we repeat the same file five times) and would like to extract values from all of them. Extracting values from a multi-layer raster objects (`RasterStack` for `raster` package) takes less time than extracting values from the individual layers one at a time. This can be observed below.   

---

**`terra::extract()`**

```{r write_raster_prism_10, echo = F}
# terra::writeRaster(prism_tmax_0701_KS_sr_10, "Data/prism_tmax_0701_KS_sr_10.tif", overwrite = TRUE)

prism_tmax_0701_KS_sr_10 <- rast("Data/prism_tmax_0701_KS_sr_10.tif")
```

```{r compare_speed_te}
#--- extract from 5 layers one at a time ---#
tic()
temp <- terra::extract(prism_tmax_0701_KS_sr_10, KS_county_sv)
temp <- terra::extract(prism_tmax_0701_KS_sr_10, KS_county_sv)
temp <- terra::extract(prism_tmax_0701_KS_sr_10, KS_county_sv)
temp <- terra::extract(prism_tmax_0701_KS_sr_10, KS_county_sv)
temp <- terra::extract(prism_tmax_0701_KS_sr_10, KS_county_sv)
toc()

#--- extract from a 5-layer SpatRaster ---#
tic()
prism_tmax_ml_5 <-
  c(
    prism_tmax_0701_KS_sr_10,
    prism_tmax_0701_KS_sr_10,
    prism_tmax_0701_KS_sr_10,
    prism_tmax_0701_KS_sr_10,
    prism_tmax_0701_KS_sr_10
  )
temp <- terra::extract(prism_tmax_ml_5, KS_county_sv)
toc()
```

---

**`exact_extract()`**

```{r compare_speed_ee}
#--- extract from 5 layers one at a time ---#
tic()
temp <- exact_extract(prism_tmax_0701_KS_sr_10, KS_county_sf, progress = FALSE)
temp <- exact_extract(prism_tmax_0701_KS_sr_10, KS_county_sf, progress = FALSE)
temp <- exact_extract(prism_tmax_0701_KS_sr_10, KS_county_sf, progress = FALSE)
temp <- exact_extract(prism_tmax_0701_KS_sr_10, KS_county_sf, progress = FALSE)
temp <- exact_extract(prism_tmax_0701_KS_sr_10, KS_county_sf, progress = FALSE)
toc()

#--- extract from from a 5-layer SpatRaster ---#
tic()
prism_tmax_stack_5 <-
  c(
    prism_tmax_0701_KS_sr_10,
    prism_tmax_0701_KS_sr_10,
    prism_tmax_0701_KS_sr_10,
    prism_tmax_0701_KS_sr_10,
    prism_tmax_0701_KS_sr_10
  )

set.names(prism_tmax_stack_5, paste0("layer_", 1:5))

temp <-
  exact_extract(
    prism_tmax_stack_5,
    KS_county_sf,
    progress = FALSE
  )
toc()
```

The reduction in computation time for both methods makes sense. Since both layers have exactly the same geographic extent and resolution, finding the polygons-cells correspondence is done once and then it can be used repeatedly across the layers for the multi-layer `SparRaster` and `RasterStack`. This clearly suggests that when you are processing many layers of the same spatial resolution and extent, you should first stack them and then extract values at the same time instead of processing them one by one as long as your memory allows you to do so. 

<!-- There is much more to discuss about the computation speed of raster data extraction for polygons. For those who are interested in this topic, go to Chapter \@ref(EE). -->