11-ras-vec.Rmd

# Raster-vector interactions {#raster-vector}

## Prerequisites {-}

- This chapter requires the following packages:

```{r, message = FALSE} 
library(sf)
library(raster)
library(tidyverse)
library(spData)
library(spDataLarge)
```

## Introduction

This section focuses on interactions between raster and data vector geographic data models, introduced in Chapter \@ref(spatial-class).
It includes four main techniques:
raster cropping and masking using vector objects (section \@ref(raster-cropping));
extracting raster values using different types of vector data (section \@ref(raster-extraction));
and raster-vector conversion (sections \@ref(rasterization) and \@ref(spatial-vectorization)).
<!-- operations are not symmetrical, for example: -->
<!-- - raster clipping - no vector counterpart -->
<!-- - raster extraction is connected to some methods used in vectorization and rasterization -->
<!-- - etc. -->
The above concepts are demonstrated using data used in previous chapters to understand their potential real-world applications.

## Raster cropping

Many geographic data projects involve integrating data from many different sources, such as remote sensing images (rasters) and administrative boundaries (vectors).
Often the extent of input raster datasets is larger than the area of interest.
In this case raster **cropping** and **masking** are useful for unifying the spatial extent of input data.
Both operations reduce object memory use and associated computational resources for subsequent analysis steps, and may be a necessary preprocessing step before creating attractive maps involving raster data.

We will use two objects to illustrate raster cropping:

- A `raster` object `srtm` representing elevation (meters above sea level) in Southwestern Utah.
- A vector (`sf`) object `zion` representing Zion National Park.

Both target and cropping objects must have the same projection.
The following code chunk therefore not only loads the datasets, from the **spDataLarge** package installed in Chapter \@ref(spatial-class).
It also reprojects `zion` (see section \@ref(reproj-geo-data) for more on reprojection):

```{r, results='hide'}
srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = read_sf(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion = st_transform(zion, projection(srtm))
```

We will use `crop()` from the **raster** package to crop the `srtm` raster.
`crop()` reduces the rectangular extent of the object passed to its first argument based on the extent of the object passed to its second argument, as demonstrated in the command below (which generates Figure \@ref(fig:cropmask):B --- note the smaller extent of the raster background):

```{r}
srtm_cropped = crop(srtm, as(zion, "Spatial"))
```

Related to `crop()` is the **raster** function `mask()`, which sets values outside of the bounds a the object passed to its second argument to `NA`.
The following command therefore masks every cell outside of the the Zion National Park boundaries (Figure \@ref(fig:cropmask):C):

```{r}
srtm_masked = mask(srtm, zion)
```

Changing the settings of `mask()` yields in different results.
Setting `maskvalue = 0`, for example, would set all pixels outside the national park to 0.
Setting `inverse = TRUE` would mask everything *inside* the bounds of the park (see `?mask` for details).

```{r cropmask, echo = FALSE, fig.cap="Illustration of raster cropping (center) and raster masking (right)."}
# TODO: split into reproducible script, e.g. in code/09-cropmask.R
library(tmap)
library(grid)
library(rcartocolor)
terrain_colors = carto_pal(7, "TealRose")

pz1 = tm_shape(srtm) + 
  tm_raster(palette = terrain_colors, title = "Elevation (m)", 
            legend.show = TRUE, auto.palette.mapping = FALSE, style = "cont") + 
  tm_shape(zion) +
  tm_borders(lwd = 2) + 
  tm_layout(legend.frame = TRUE, legend.position = c("right", "top"))

pz2 = tm_shape(srtm_cropped) + 
  tm_raster(palette = terrain_colors, title = "Elevation (m)", 
            legend.show = TRUE, auto.palette.mapping = FALSE, style = "cont") + 
  tm_shape(zion) +
  tm_borders(lwd = 2) + 
  tm_layout(legend.frame = TRUE, legend.position = c("right", "top"))

pz3 = tm_shape(srtm_masked) + 
  tm_raster(palette = terrain_colors, title = "Elevation (m)",
            legend.show = TRUE, auto.palette.mapping = FALSE, style = "cont") + 
  tm_shape(zion) +
  tm_borders(lwd = 2) + 
  tm_layout(legend.frame = TRUE, legend.position = c("right", "top"))

grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 3, heights = unit(c(0.25, 5), "null"))))
grid.text("A. Original", vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
grid.text("B. Cropped", vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
grid.text("C. Masked", vp = viewport(layout.pos.row = 1, layout.pos.col = 3))
print(pz1, vp = viewport(layout.pos.row = 2, layout.pos.col = 1))
print(pz2, vp = viewport(layout.pos.row = 2, layout.pos.col = 2))
print(pz3, vp = viewport(layout.pos.row = 2, layout.pos.col = 3))
```

## Raster extraction

Raster extraction is the process of identifying and returning the values associated with a 'target' raster at specific locations, based on a (typically vector) geographic 'selector' object.
The results depend on the type of selector used (points, lines or polygons) and arguments passed to the `raster::extract()` function, which we use to demonstrate raster extraction.
The reverse of raster extraction --- assigning raster cell values based on vector objects --- is rasterization, described in section \@ref(rasterization).

The simplest example is extracting the value of a raster cell at specific **points**.
For this purpose we will use `zion_points`, which contain a sample of 30 locations with the Zion National Park (Figure \@ref(fig:pointextr)). 
<!-- They could represent places where soils properties were measured and we want to know what is the elevation of each point. -->
The following command extracts elevation values from `srtm` and assigns the resulting vector to a new column (`elevation`) in the `zion_points` dataset: 

```{r}
zion_points$elevation = raster::extract(srtm, zion_points)
```

```{r, echo=FALSE, eval=FALSE}
# Aim: demonstrate buffer arg in raster extract
elev_b1 = raster::extract(srtm, zion_points, buffer = 1000)
```

The `buffer` argument can be used to specify a buffer radius (in meters) around each point.
The result of `raster::extract(srtm, zion_points, buffer = 1000)`, for example, is a list of vectors, each of which representing the values of cells inside the buffer associated with each point.
In practice this example is a special case of extraction with a polygon selector, described below.

```{r pointextr, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="Locations of points used for raster extraction."}
library(tmap)
library(grid)

rast_poly_point = tm_shape(srtm) +
  tm_raster(palette = terrain_colors, title = "Elevation (m)", 
            legend.show = TRUE, auto.palette.mapping = FALSE, style = "cont") + 
  tm_shape(zion) +
  tm_borders(lwd = 2) + 
  tm_shape(zion_points) + 
  tm_symbols(shape = 1, col = "black") + 
  tm_layout(legend.frame = TRUE, legend.position = c("right", "top"))
rast_poly_point
```

Raster extraction also works with **line** selectors.
To demonstrate this, the code below creates `zion_transect`, a straight line going from northwest to southeast of the Zion National Park, illustrated in Figure \@ref(fig:lineextr):A (see section \@ref(vector-data) for a recap on the vector data model):

```{r}
zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) %>%
  st_linestring() %>% 
  st_sfc() %>% 
  st_sf()
```


```{r, eval=FALSE, echo=FALSE}
# Aim: show how extraction works with non-straight lines by
# using this alternative line object:
zion_transect = cbind(c(-113.2, -112.9, -113.2), c(37.45, 37.2, 37.5)) %>%
  st_linestring() %>% 
  st_sfc() %>% 
  st_sf()
```

The utility of extracting heights from a linear selector is illustrated by imagining that you are planning a hike.
The method demonstrated below provides an 'elevation profile' of the route (the line does not need to be straight), useful for estimating how long it will take due to long climbs:

```{r}
transect = raster::extract(srtm, zion_transect, along = TRUE, cellnumbers = TRUE)
```

Note the use of `along = TRUE` and `cellnumbers = TRUE` arguments to return cell IDs *along* the path. 
The result is a list containing a matrix of cell IDs in the first column and elevation values in the second.
The subsequent code chunk first converts this tricky matrix-in-a-list object into a simple data frame, returns the coordinates associated with each extracted cell and finds the associated distances along the transect (see `?geosphere::distm()` for details):

```{r}
transect_df = map_dfr(transect, as_data_frame, .id = "ID")
transect_coords = xyFromCell(srtm, transect_df$cell)
transect_df$dist = c(0, cumsum(geosphere::distGeo(transect_coords)))    
```

The resulting `transect_df` can be used to create elevation profiles, as illustrated in Figure \@ref(fig:lineextr):B.

```{r lineextr, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="Location of a line used for raster extraction (left) and the elevation along this line (right)."}
library(tmap)
library(grid)

zion_transect_points = st_cast(zion_transect, "POINT")[1:2, ]
zion_transect_points$name = c("start", "end")

rast_poly_line = tm_shape(srtm) +
  tm_raster(palette = terrain_colors, title = "Elevation (m)", 
            legend.show = TRUE, auto.palette.mapping = FALSE, style = "cont") + 
  tm_shape(zion) +
  tm_borders(lwd = 2) + 
  tm_shape(zion_transect) + 
  tm_lines(col = "black", lwd = 4) + 
  tm_shape(zion_transect_points) +                                                 
  tm_text("name", bg.color = "white", bg.alpha = 0.75, auto.placement = TRUE) +
  tm_layout(legend.frame = TRUE, legend.position = c("right", "top"))

plot_transect = ggplot(transect_df, aes(dist, srtm)) + 
  geom_line() +
  labs(x = "Distance (m)", y = "Elevation (m a.s.l.)") + 
  theme_bw() +
  theme(plot.margin = unit(c(5.5, 15.5, 5.5, 5.5), "pt"))

grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2, heights = unit(c(0.25, 5), "null"))))
grid.text("A. Line extraction", vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
grid.text("B. Elevation along the line", vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
print(rast_poly_line, vp = viewport(layout.pos.row = 2, layout.pos.col = 1))
print(plot_transect, vp = viewport(layout.pos.row = 2, layout.pos.col = 2))
```

The final type of geographic vector object for raster extraction is **polygons**.
Like lines and buffers, polygons tend to return many raster values per polygon.
This is demonstrated in the command below, which results in a data frame with column names  `ID` (the row number of the polygon) and `srtm` (associated elevation values):

```{r, eval=FALSE, echo=FALSE}
# aim: create zion_many to test multi-polygon results
n = 3
zion_many = st_sample(x = zion, size = n) %>% 
  st_buffer(dist = 0.01) %>% 
  st_sf(data.frame(v = 1:n), geometry = .)
plot(zion_many)

# for continuous data:
zion_srtm_values = raster::extract(x = srtm, y = zion_many, df = TRUE)
group_by(zion_srtm_values, ID) %>% 
  summarize_at(vars(srtm), funs(min, mean, max))

# for categories
zion_nlcd = raster::extract(nlcd, zion_many, df = TRUE, factors = TRUE)
dplyr::select(zion_nlcd, ID, levels) %>% 
  gather(key, value, -ID) %>%
  group_by(ID, key, value) %>%
  tally() %>% 
  spread(value, n, fill = 0)
```


```{r}
zion_srtm_values = raster::extract(x = srtm, y = zion, df = TRUE)
```

Such results can be used to generate summary statistics for raster values per polygon, for example to  to characterize a single region or to compare many regions.
The generation of summary statistics is demonstrated the code below, which creates the object `zion_srtm_df` containing summary statistics for elevation values in Zion National Park (see \@ref(fig:polyextr):A):

```{r}
group_by(zion_srtm_values, ID) %>% 
  summarize_at(vars(srtm), funs(min, mean, max))
```

The preceding code chunk used the **tidyverse** to provide summary statistics for cell values per polygon ID, as described in Chapter \@ref(attr).
The results provide useful summaries, for example that the maximum height in the park is around 2,661 meters (other summary statistics such as standard deviation can also be calculated in this way).
Because there is only one polygon in the example a data frame with a single row is returned, but the method works when multiple selector polygons are used.

The same approach works for counting occurrences of categorical raster values within polygons.
This is illustrated with a land cover dataset (`nlcd`) from the **spDataLarge** package in \@ref(fig:polyextr):B and demonstrated in the code below:

```{r, warning=FALSE, message=FALSE}
zion_nlcd = raster::extract(nlcd, zion, df = TRUE, factors = TRUE)
dplyr::select(zion_nlcd, ID, levels) %>% 
  gather(key, value, -ID) %>%
  group_by(ID, key, value) %>%
  tally() %>% 
  spread(value, n, fill = 0)
```

```{r polyextr, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="Area used for continuous (left) and categorical (right) raster extraction."}
library(tmap)
library(grid)

rast_poly_srtm = tm_shape(srtm) +
  tm_raster(palette = terrain_colors, title = "Elevation (m)", 
            legend.show = TRUE, auto.palette.mapping = FALSE, style = "cont") + 
  tm_shape(zion) +
  tm_borders(lwd = 2) +
  tm_layout(legend.frame = TRUE, legend.position = c("left", "bottom"))

landcover_cols = c("#476ba0", "#aa0000", "#b2ada3", "#68aa63", "#a58c30", "#c9c977", "#dbd83d", "#bad8ea")

rast_poly_nlcd = tm_shape(nlcd) +
  tm_raster(palette = landcover_cols, style = "cat", title = "Land cover", 
            legend.show = TRUE) + 
  tm_shape(zion) +
  tm_borders(lwd = 2) +
  tm_layout(legend.frame = TRUE, legend.position = c("left", "bottom"))

grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2, heights = unit(c(0.25, 5), "null"))))
grid.text("A. Continuous data extraction", vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
grid.text("B. Categorical data extraction", vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
print(rast_poly_srtm, vp = viewport(layout.pos.row = 2, layout.pos.col = 1))
print(rast_poly_nlcd, vp = viewport(layout.pos.row = 2, layout.pos.col = 2))
```

So far we have seen how `raster::extract()` is a flexible way of extracting raster cell values from a range of input geographic objects.
An issue with the function, however, is that it is slow.
If this is a problem it is useful to know about alternatives and work-arounds, three of which are presented below.

- **Parallelization**: this approach works when using many geographic vector selector objects by splitting them into groups and extracting cell values independently for each group (see `?raster::clusterR()` for details of this approach).
<!-- tabularaster (ref to the vignette - https://cran.r-project.org/web/packages/tabularaster/vignettes/tabularaster-usage.html)-->
- Use the **velox** package [@hunziker_velox:_2017], which provides a fast method for extracting raster data that fits in memory (see the packages [`extract`](https://hunzikp.github.io/velox/extract.html) vignette for details).
- Using **R-GIS bridges** (see Chapter \@ref(gis)): efficient calculation of raster statistics from polygons can be found in the SAGA function `saga:gridstatisticsforpolygons`, for example, which can be accessed via **RQGIS**.
<!-- Methods similar to `raster::extract` can be found in GRASS GIS (e.g. v.rast.stats) -->
<!-- https://grass.osgeo.org/grass74/manuals/v.rast.stats.html - test -->
<!-- https://twitter.com/mdsumner/status/976978499402571776 -->
<!-- https://gist.github.com/mdsumner/d0b26238321a5d2c2c2ba663ff684183 -->

## Rasterization {#rasterization}

Rasterization is the conversion of vector objects into their representation in raster objects.
Usually, the output raster is used for quantitative analysis (e.g. analysis of terrain) or modeling.
As we saw in Chapter \@ref(spatial-class) the raster data model has some characteristics that make it conducive to certain methods.
Furthermore, the process of rasterization can help simplify datasets because the resulting values all have the same spatial resolution: rasterization can be seen as a special type of geographic data aggregation.

The **raster** package contains the function `rasterize()` for doing this work.
Its first two arguments are `x`, vector object to be rasterized and `y`, a 'template raster' object defining the extent, resolution and CRS of the output.
The geographic resolution of the input raster has a major impact on the results: if it is too low (cell size is too large) the result may miss the full geographic variability of the vector data; if it is too high computational times may be excessive.
There are no simple rules to follow when deciding an appropriate geographic resolution, which is heavily dependent on the intended use of the results.

To demonstrate rasterization in action, we will use a template raster that has the same extent and CRS as the input vector data `cycle_hire_osm_projected` (a dataset on cycle hire points in London illustrated in Figure \@ref(fig:vector-rasterization1):A) and spatial resolution of 1000 meters:

```{r}
cycle_hire_osm_projected = st_transform(cycle_hire_osm, 27700)
raster_template = raster(extent(cycle_hire_osm_projected), resolution = 1000,
                         crs = st_crs(cycle_hire_osm_projected)$proj4string)
```

Rasterization is a very flexible operation: the results depend not only on the nature of the template raster, but also on the type of input vector (e.g. points, polygons) and a variety arguments taken by the `rasterize()` function.

To illustrate this flexibility we will try three different approaches rasterization.
First we create a raster representing the presence or absence of cycle hire points (known as presence/absence rasters).
In this case `rasterize()` requires only one argument in addition to `x` and `y` (the aformentioned vector and raster objects): a value to be transferred to all non-empty cells specified by `field` (results illustrated Figure \@ref(fig:vector-rasterization1):B).

```{r}
ch_raster1 = rasterize(cycle_hire_osm_projected, raster_template, field = 1)
```

The `fun` argument specifies summary statistics used to covert multiple observations in close proximity into associate cells in the raster object.
By default `fun = 'last` is used but other options such as `fun = "count"` can be used,  in this case to count the number of cycle hire points in each grid cell (the results of this operation are illustrated in Figure \@ref(fig:vector-rasterization1):C).

```{r}
ch_raster2 = rasterize(cycle_hire_osm_projected, raster_template, 
                       field = 1, fun = "count")
```

The new output, `ch_raster2`, shows the number of cycle hire points in each grid cell.
The cycle hire locations have different numbers of bicycles described by the `capacity` variable, raising the question, what's the capacity in each grid cell?
To calculate that we must `sum` the field (`"capacity"`), resulting in output illustrated in Figure \@ref(fig:vector-rasterization1):D, calculated with the following command (other summary functions such as `mean` could be used):

```{r}
ch_raster3 = rasterize(cycle_hire_osm_projected, raster_template, 
                       field = "capacity", fun = sum)
```

```{r vector-rasterization1, echo=FALSE, fig.cap="Examples of point's rasterization.", warning=FALSE}
source("code/11-vector-rasterization1.R")
```

Another dataset based on California's polygons and borders (created below) illustrates raterization of lines.
After casting the polygon objects into a multilinestring, a template raster is created, with a resolution of a 0.5 degree:

```{r}
california = dplyr::filter(us_states, NAME == "California")
california_borders = st_cast(california, "MULTILINESTRING")
raster_template2 = raster(extent(california), resolution = 0.5,
                         crs = st_crs(california)$proj4string)
```

Line rasterization is demonstrated in the code below.
In the resulting raster, all cells that are touched by a line get a value, as illustrated in Figure \@ref(fig:vector-rasterization2):A.

```{r}
california_raster1 = rasterize(california_borders, raster_template2)
```

Polygon rasterization, by contrast, selects only cells whose centroids are inside the selector polygon, as illustrated in Figure \@ref(fig:vector-rasterization2):B.

```{r}
california_raster2 = rasterize(california, raster_template2)
```

<!-- getCover? -->
<!-- the fraction of each grid cell that is covered by the polygons-->
<!-- ```{r, echo=FALSE, eval=FALSE} -->
<!-- california_raster3 = rasterize(california, raster_template2, getCover = TRUE) -->
<!-- r3po = tm_shape(california_raster3) + -->
<!--   tm_raster(legend.show = TRUE, title = "Values: ", style = "fixed", breaks = c(0, 1, 25, 50, 75, 100)) + -->
<!--   tm_shape(california) + -->
<!--   tm_borders() + -->
<!--   tm_layout(outer.margins = rep(0.01, 4), -->
<!--             inner.margins = rep(0, 4)) -->
<!-- ``` -->

<!-- It is also possible to use the `field` or `fun` arguments for lines and polygons rasterizations. -->

```{r vector-rasterization2, echo=FALSE, fig.cap="Examples of line and polygon rasterizations.", warning=FALSE}
source("code/11-vector-rasterization2.R")
```

As with `raster::extract()`,  `raster::rasterize()` works well for most cases but is not performance optimized. 
Fortunately, there are several alternatives, including the `fasterize::fasterize()` and `gdalUtils::gdal_rasterize()`. 
The former is much (100 times+) faster than `rasterize()` but is currently limited to polygon rasterization.
The latter is part of GDAL and therefore requires a vector file (instead of an `sf` object) and rasterization parameters (instead of a `Raster*` template object) as inputs.^[See more at http://www.gdal.org/gdal_rasterize.html.]

## Spatial vectorization

Spatial vectorization is the counterpart of rasterization \@ref(rasterization), but in the opposite direction.
It involves converting spatially continuous raster data into spatially discrete vector data such as points, lines or polygons.

```{block2 type="rmdnote"}
Be careful with the wording!
In R vectorization refers to the possibility of replacing `for`-loops and alike by doing things like `1:10 / 2` (see also @wickham_advanced_2014).
```

The simplest form of vectorization is to convert the centroids of raster cells into points.
`rasterToPoints()` does exactly this for all non-`NA` raster grid cells (Figure \@ref(fig:raster-vectorization1)).
Setting the `spatial` parameter to `TRUE` ensures the output is a spatial object, not a matrix.

```{r}
elev_point = rasterToPoints(elev, spatial = TRUE) %>% 
  st_as_sf()
```

```{r raster-vectorization1, echo=FALSE, fig.cap="Raster and point representation of `elev`.", warning=FALSE}
source("code/11-raster-vectorization1.R")
```

Another common type of spatial vectorization is the creation of contour lines representing lines of continuous height or temperatures (isotherms) for example.
We will use a real-world digital elevation model (DEM) because the artificial raster `elev` produces parallel lines (task: verify this and explain why this happens).
<!-- because when creating it we made the upper left corner the lowest and the lower right corner the highest value while increasing cell values by one from left to right. -->
Contour lines can be created with the **raster** function `rasterToContour()`, which is itself a wrapper around `contourLines()`, as demonstrated below:

```{r, eval = FALSE}
# not shown
data(dem, package = "RQGIS")
plot(dem, axes = FALSE)
contour(dem, add = TRUE)
```

Contours can be added to existing plots with functions such as `contour()`, `rasterVis::contourplot()` or `tmap::tm_iso()`.
As illustrated in Figure \@ref(fig:contour) (which was created using the **tmap** package described in Chapter \@ref(adv-map)), isolines can be labelled.

```{r contour, echo=FALSE, message=FALSE, fig.cap = "DEM hillshade of the southern flank of Mt. Mongón overlaid with contour lines.", warning=FALSE}
library(tmap)
data("dem", package = "RQGIS")
# create hillshade
hs = hillShade(slope = terrain(dem, "slope"), aspect = terrain(dem, "aspect"))
# create contour
cn = rasterToContour(dem)
rect = tmaptools::bb_poly(hs)
bbx = tmaptools::bb(hs, xlim = c(-.02, 1), ylim = c(-.02, 1), relative = TRUE)

tm_shape(hs, bbox = rect) +
	tm_grid(col = "black", n.x = 2, n.y = 2, labels.inside.frame = FALSE,
	        labels.rot = c(0, 90)) +
	tm_raster(palette = gray(0:100 / 100), n = 100, legend.show = FALSE) +
	tm_shape(dem) +
	tm_raster(alpha = 0.5, palette = terrain.colors(25),
	          auto.palette.mapping = FALSE, legend.show = FALSE) +
	tm_shape(cn) + 
	tm_lines(col = "white") +
	tm_text("level", col = "white") +
	qtm(rect, fill = NULL) +
	tm_layout(outer.margins = c(0.04, 0.04, 0.02, 0.02), frame = FALSE)
```

The final type of vectorisation involves conversion of rasters to polygons.
This can be done with `raster::rasterToPolygons()`, wich converts each raster cell into a polygon consisting of five coordinates, all of which are stored in memory (explaining why rasters are often fast compared with vectors!)

This is illustrated below by converting the `grain` object into polygons and subsequently dissolving borders between polygons with the same attribute values.
Attributes in this case are stored in a collumn called `layer` (see section \@ref(geometry-unions) and Figure \@ref(fig:raster-vectorization2)).
(Note: a convenient alternative for converting rasters into polygons is `spex::polygonize()` which by default returns an `sf` object.)

```{r}
grain_poly = rasterToPolygons(grain) %>% 
  st_as_sf()
grain_poly2 = grain_poly %>% 
  group_by(layer) %>%
  summarize()
```

```{r raster-vectorization2, echo=FALSE, fig.cap="Illustration of vectorization of raster (left) into polygon (center) and polygon aggregation (right).", warning=FALSE}
source("code/11-raster-vectorization2.R")
```

<!-- ## distances? -->

## Exercises

The next two exercises will use a vector (`random_points`) and raster dataset (`ndvi`) from the **RQGIS** package.
It also uses a polygonal 'convex hull' derived from the vector dataset (`ch`) to represent the area of interest:
```{r, message=FALSE}
library(RQGIS)
data(random_points)
data(ndvi)
ch = st_combine(random_points) %>% 
  st_convex_hull()
```
1. Crop the `ndvi` raster using (1) the `random_points` dataset and (2) the `ch` dataset.
Are there any difference in the output maps?
Next, mask `ndvi` using these two datasets.
Can you see any difference now?
How can you explain that?
```{r, echo=FALSE, eval=FALSE}
plot(ndvi)
plot(st_geometry(random_points), add = TRUE)
plot(ch, add = TRUE)

ndvi_crop1 = crop(ndvi, as(random_points, "Spatial")) 
ndvi_crop2 = crop(ndvi, as(ch, "Spatial")) 
plot(ndvi_crop1)
plot(ndvi_crop2)

ndvi_mask1 = mask(ndvi, as(random_points, "Spatial")) 
ndvi_mask2 = mask(ndvi, as(ch, "Spatial")) 
plot(ndvi_mask1)
plot(ndvi_mask2)
```
1. Firstly, extract values from `ndvi` at the points represented in `random_points`.
Next, extract average values of `ndvi` using a 90 buffer around each point from `random_points` and compare these two sets of values. 
When would extracting values by buffers be more suitable than by points alone?
```{r, echo=FALSE, eval=FALSE}
random_points_buf = st_buffer(random_points, dist = 90)
plot(ndvi)
plot(st_geometry(random_points_buf), add = TRUE)
plot(ch, add = TRUE)
random_points$ndvi = extract(ndvi, random_points, buffer = 90, fun = mean)
random_points$ndvi2 = extract(ndvi, random_points)
plot(random_points$ndvi, random_points$ndvi2)
```
1. Subset points higher than 3100 meters in New Zealand (the `nz_height` object) and create a template raster with a resolution of 3km. 
Using these objects:
    - Count numbers of the highest points in each grid cell.
    - Find the maximum elevation in each grid cell.
```{r, echo=FALSE, eval=FALSE}
nz_height3100 = dplyr::filter(nz_height, elevation > 3100)
new_graticule = st_graticule(nz_height3100, datum = 2193)
plot(nz_height3100$geometry, graticule = new_graticule, axes = TRUE)
nz_template = raster(extent(nz_height3100), resolution = 3000,
                         crs = st_crs(nz_height3100)$proj4string)
nz_raster = rasterize(nz_height3100, nz_template, 
                       field = "elevation", fun = "count")
plot(nz_raster)
nz_raster2 = rasterize(nz_height3100, nz_template, 
                       field = "elevation", fun = max)
plot(nz_raster2)
```
1. Polygonize the `grain` dataset and filter all squares representing clay.
    - Name two advantages and disadvantages of vector data over raster data.
    -  At which points would it be useful to convert rasters to vectors in your work?

```{r, echo=FALSE, eval=FALSE}
grain_poly = rasterToPolygons(grain) %>% 
  st_as_sf()
levels(grain)
clay = dplyr::filter(grain_poly, layer == 1)
plot(clay)
# advantages: can be used to subset other vector objects
# can do affine transformations and use sf/dplyr verbs
# disadvantages: better consistency, fast processing on some operations, functions developed for some domains
```