rasterR_slides.knit

class: right, bottom

## Thinking in Rasters

##### Nikhil Kaza
##### Department of City & Regional Planning <br /> University of North Carolina at Chapel Hill
###### updated: 2026-09-22

---

# Fields vs. Objects

This is an age old debate as to which is a better representation of reality.

Objects: Discrete, sharply defined boundaries, have distinct attributes (e.g. buildings, roads, parcels, census tracts)
  - Need to keep track of topology and spatial relations among them

Fields: Something that varies continuously over space. The discretisation is an artefact of the data storage and representation (e.g. temperature, ground level, urbanicity)
  - Topological relationships are embedded

Usually fields are represented by rasters.

---

# What is a raster, really?

A raster is a grid of cells (pixels), each holding a value. Rasters show up well outside geography:

<div style="display:flex; justify-content:center; gap:3%; margin-top:20px;">
<div style="text-align:center; width:22%;">
<img src="./includes/what_is_a_raster/photograph.jpg" style="width:100%; border-radius:4px;">
<div style="font-size:0.8em; margin-top:6px;"><b>Photograph</b><br/>light intensity</div>
</div>
<div style="text-align:center; width:22%;">
<img src="./includes/what_is_a_raster/xray.jpg" style="width:100%; border-radius:4px;">
<div style="font-size:0.8em; margin-top:6px;"><b>X-ray</b><br/>tissue density</div>
</div>
<div style="text-align:center; width:22%;">
<img src="./includes/what_is_a_raster/telescope.gif" style="width:100%; border-radius:4px;">
<div style="font-size:0.8em; margin-top:6px;"><b>Telescope image</b><br/>photon counts</div>
</div>
<div style="text-align:center; width:22%;">
<img src="./includes/what_is_a_raster/drawing.jpg" style="width:100%; border-radius:4px;">
<div style="font-size:0.8em; margin-top:6px;"><b>Scanned drawing</b><br/>ink density</div>
</div>
</div>
<div style="font-size:0.55em; color:#999; text-align:right; margin-top:2px;">Stephan's Quintet, Hubble vs. James Webb — NASA, ESA, CSA, STScI</div>

None of these are geospatial, but they have a coordinate system.

---

# Key features of any raster

- **Grid structure** — rows and columns, forming a rectangular array

- **Pixel values** — each cell holds a number: brightness, elevation, intensity, a class label

- **Resolution** — the size of each cell; smaller cells mean finer detail

- **Bands/layers** — a raster can carry several layers at once (e.g. R, G, B channels, or timestamps), and layers can also stack observations of the same place over time

- **Spatial context** *(optional)* — if a raster is georeferenced, the coordinates of its corners plus its resolution locate every pixel in the world. If not, it's just a grid.

---

# Raster data in urban analytics

.pull-left[
- Satellite images

- Ground level pictures

- Airborne sensors

- Radar
]

.pull-right[
![](./figs/examplerasters.png)
[Image Credit: ESRI](https://developers.arcgis.com/net/latest/uwp/guide/add-raster-data.htm)
]

---
class: inverse, middle, center

# Part I (continued)
## Four ways to think about a raster

---

# 1. A raster as a table

Each pixel is a row. Each band is a column.  
- **A raster is a tibble that happens to have a grid layout.**

---

# 2. A raster as a matrix

Equivalently: each band is its own matrix, and a multi-band raster is a stack of matrices of identical dimensions.

---

# 3. A raster as a graph

Pixels are nodes. Edges connect neighboring pixels. The definition of "neighbour" (rook, queen, k-th order) is a modelling choice.

This framing matters whenever an operation depends on *connectivity* rather than *position* — moving-window filters, cost-distance and regionalisation

---

# 4. A raster as sampled continuity

A raster can also be seen as a discrete sample of an underlying continuous surface — elevation, temperature, reflectance — independent of how finely we chose to grid it.

Resolution is then a sampling-rate decision, and operations like resampling/aggregation are really questions about how to re-sample that surface.

---

# Which framing, when?

- Reaching for `dplyr`-style verbs, joining pixel values to other data, tidy plotting → **table**

- Linear algebra, image filters expressed as matrix operations → **matrix**

- Connectivity, adjacency, shortest paths, region-growing → **graph**

- Interpolation, resampling, "what's the true value between these cells?" → **continuous surface**

None of these is the "correct" one. Pick the framing that matches the question you're asking.

---
class: inverse, middle, center

# Part II
## The raster ecosystem in R
### (it's bigger than geospatial)

---

# Remote sensing & raster visualization

<table class="table" style="font-size: 22px; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> Package </th>
   <th style="text-align:left;"> Description </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> terra </td>
   <td style="text-align:left;"> Modern package for handling spatial raster data efficiently. Supports cropping, resampling, projection. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> rasterVis </td>
   <td style="text-align:left;"> Lattice-style raster visualization and thematic mapping. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> stars </td>
   <td style="text-align:left;"> Spatiotemporal arrays; integrates with sf and the tidyverse. Good for multi-band satellite data. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> satellite </td>
   <td style="text-align:left;"> Satellite imagery processing, including atmospheric correction and vegetation indices. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> RStoolbox </td>
   <td style="text-align:left;"> PCA, supervised classification, spectral indices for remote sensing. </td>
  </tr>
</tbody>
</table>

---

# Environmental & climate raster data

<table class="table" style="font-size: 22px; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> Package </th>
   <th style="text-align:left;"> Description </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> ncdf4 / RNetCDF </td>
   <td style="text-align:left;"> Reading and writing NetCDF files; custom extraction of structured arrays. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> stars </td>
   <td style="text-align:left;"> Also handles NetCDF and spatiotemporal rasters within a tidy workflow. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> geodata </td>
   <td style="text-align:left;"> Global environmental datasets: WorldClim, elevation, land cover. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> elevatr </td>
   <td style="text-align:left;"> Global elevation data from SRTM, AWS, and other sources. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> climate / rnoaa </td>
   <td style="text-align:left;"> Climate and weather time series, NOAA datasets. </td>
  </tr>
</tbody>
</table>

---

# Biological & medical imaging

The pixel grid is the same object as your satellite scene. The domain-specific vocabulary is what changes.

---

# General image processing & computer vision

<table class="table" style="font-size: 22px; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> Package </th>
   <th style="text-align:left;"> Description </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> imager </td>
   <td style="text-align:left;"> Pixel-level operations, filters, edge detection, convolution. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> magick </td>
   <td style="text-align:left;"> Wraps ImageMagick: cropping, compositing, format conversion, annotation. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> OpenImageR </td>
   <td style="text-align:left;"> Preprocessing, feature extraction, similarity metrics for ML pipelines. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> landscapemetrics </td>
   <td style="text-align:left;"> Landscape-ecology metrics (patch, class, landscape level) from raster data. </td>
  </tr>
  <tr>
   <td style="text-align:left;"> imagefx </td>
   <td style="text-align:left;"> Lightweight blur/sharpen/color-adjustment effects. </td>
  </tr>
</tbody>
</table>

---
class: inverse, middle, center

# Part III
## `terra` basics
### thought about the tidyverse way

---

# Why `terra`?

- Successor to Robert Hijmans' `raster` package (`raster` is now in maintenance mode)

- One class, `SpatRaster`, replaces `RasterLayer` / `RasterStack` / `RasterBrick`

- Fast, out-of-memory processing — handles rasters larger than RAM

- Plays well with `sf` (vectors) and the tidyverse (`as.data.frame()`, pipes)

We'll use `terra` throughout. Some older packages (e.g. `gdistance`) still expect the legacy `raster` classes.

---

# Creating a `SpatRaster`

``` r
(rast_empty <- rast())
```

```
## class       : SpatRaster
## size        : 180, 360, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
```

---

``` r
res(rast_empty) <- 20
rast_empty
```

```
## class       : SpatRaster
## size        : 9, 18, 1  (nrow, ncol, nlyr)
## resolution  : 20, 20  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
```

``` r
ncol(rast_empty) <- 30
rast_empty
```

```
## class       : SpatRaster
## size        : 9, 30, 1  (nrow, ncol, nlyr)
## resolution  : 12, 20  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
```

---

# External files

``` r
r1 <- rast(system.file("external/test.grd", package="raster"))
plot(r1, plg=list(cex=1.5))
```

![](figs/unnamed-chunk-11-1.png)

---

# The pixel-table equivalence

Way 1 said a raster *is* a table.

``` r
as.data.frame(r1, xy = TRUE) %>%
  glimpse()
```

```
## Rows: 3,178
## Columns: 3
## $ x    <dbl> 181140, 181180, 181100, 181140, 181180, 181220, 181060, 181100, 1…
## $ y    <dbl> 333620, 333620, 333580, 333580, 333580, 333580, 333540, 333540, 3…
## $ test <dbl> 655.5997, 612.2552, 753.8189, 672.5580, 616.4962, 604.4344, 907.7…
```

Every pixel is a row, `x`/`y` are its coordinates, `test` is its value. From here on, ordinary `dplyr` verbs apply.

---

# Multiple layers = multiple columns

``` r
noise <- init(r1, fun = runif)          # a same-shaped random layer
rastS <- c(r1, noise)                    # bind layers, like bind_cols()
names(rastS) <- c("Original", "Noise")
rsum(rastS)
```

```
## class      : SpatRaster
## size       : 115, 80, 2  (nrow, ncol, nlyr)
## resolution : 40, 40  (x, y)
## names      : Original, Noise
## min values : 138.7071, 2e-04
## max values : 1736.0579, 1
```

---

``` r
plot(rastS)
```

![](figs/unnamed-chunk-14-1.png)

---

# Cell-by-cell algebra

Standard operators (`+`, `-`, `>=`, `max`, ...) work cell-by-cell across layers or against a number — exactly like `mutate()` combining columns of a tibble.

``` r
rastB <- c(r1 * 10, r1 * noise)
names(rastB) <- c('Original', 'Modified')
rsum(rastB)
```

```
## class      : SpatRaster
## size       : 115, 80, 2  (nrow, ncol, nlyr)
## resolution : 40, 40  (x, y)
## names      : Original, Modified
## min values : 1387.0708, 0.0643
## max values : 17360.5798, 1332.2969
```

---

``` r
(rastB[["Original"]] / 100 >= rastB[["Modified"]]) %>%
  plot(main = "Original/100 >= Modified")
```

---

# Moving windows

`focal()` slides a window over the 2-D grid — the spatial cousin of a rolling window over a 1-D time series.

.pull-left[

``` r
# 3x3 mean filter, in space
smoothed <- focal(r1, w = matrix(1/9, 3, 3))
round(values(smoothed)[1:8], 1)
```

```
## [1] NaN NaN NaN NaN NaN NaN NaN NaN
```
]

.pull-right[

``` r
# 3-point mean filter, in time
x <- rnorm(10)
round(zoo::rollmean(x, k = 3, fill = NA), 1)
```

```
##  [1]   NA -0.1  0.1 -0.4 -0.7 -0.6 -0.1  0.3  0.1   NA
```
]

Same idea — local averaging over a neighborhood — different notion of "neighbor" (grid cells vs. adjacent time points).

---

# The same mechanism, side by side

A window of fixed size slides across the data, one step at a time, and summarizes whatever falls inside it. Whether "next" means the next time point or the next grid cell is the only real difference.

![](./includes/moving_window_parallel.gif)

---

# Moving window: local maximum

``` r
rast3 <- focal(r1, w=matrix(1, nrow=5,ncol=5), fun=max, na.rm=TRUE)
res(rast3)
```

```
## [1] 40 40
```

``` r
plot(rast3)
```

---

# Reduce resolution

`aggregate()` pools blocks of cells into one, using a summary function — the raster analogue of grouping rows and summarizing.

``` r
rast3c <- r1 %>% aggregate(fact=5, fun=mean)
res(rast3c)
```

```
## [1] 200 200
```

``` r
plot(rast3c)
```

---

# Categorical rasters

``` r
rastCat <- classify(
  r1,
  rcl = matrix(c(-Inf,300,1,  300,800,2,  800,1100,3,  1100,Inf,4),
               ncol = 3, byrow = TRUE)
)
levels(rastCat) <- data.frame(id = 1:4, class = c('Low','Mid','High','V. High'))
```
<img src="figs/unnamed-chunk-23-1.png" alt="" width="60%" />

The reclassification matrix (breakpoint, breakpoint, new value) is doing exactly what a `case_when()` ladder of conditions would do, one row at a time.

---

# Vector/Raster operations

`terra` objects (`SpatRaster`, `SpatVector`) interoperate with `sf`. Build two polygons and bring them into `terra`:

``` r
cds1 <- structure(c(179169.418312112, 179776.022915626, 179344.40040928,
179064.429053812, 178737.795805766, 179076.094526956, 178726.130332621,
179169.418312112, 330830.748571234, 330527.446269478, 330049.161870553,
329955.838085397, 330037.496397409, 330247.47491401, 330387.460591744,
330830.748571234), .Dim = c(8L, 2L))

cds2 <- structure(c(180440.954884862, 180522.613196873, 179857.681227637,
179951.005012793, 180440.954884862, 332417.252918885, 331717.324530216,
331554.007906193, 331962.29946625, 332417.252918885), .Dim = c(5L, 2L))

polys_sf <- st_sf(ID = c(1,2),
                   geometry = st_sfc(st_polygon(list(cds1)), st_polygon(list(cds2)),
                                      crs = crs(rastCat)))
polys <- vect(polys_sf)
```

---

# Vector/Raster operations

``` r
cat("class   : SpatVector\n",
    "geoms   : ", length(polys), " polygons\n",
    "extent  : ", paste(round(as.vector(ext(polys)), 1), collapse = ", "), "\n",
    "fields  : ", paste(names(polys), collapse = ", "), "\n",
    sep = "")
```

```
## class   : SpatVector
## geoms   : 2 polygons
## extent  : 178726.1, 180522.6, 329955.8, 332417.3
## fields  : ID
```

---

# Vector/Raster operations

``` r
plot(rastCat, col = brewer.pal(n = 4, name = 'Dark2'))
plot(polys, lwd=4, border='blue', add=TRUE)
```

---

# Count percentage of cells

`extract()` returns one row per cell per polygon — a long table you then `group_by()` and `summarize()`, just as you would after a spatial join.

``` r
terra::extract(rastCat, polys) %>%
  group_by(ID, class) %>%
  summarize(count = n(), .groups = "drop_last") %>%
  mutate(pct = count / sum(count))
```

```
## # A tibble: 6 × 4
## # Groups:   ID [2]
##      ID class   count     pct
##   <dbl> <fct>   <int>   <dbl>
## 1     1 Low        92 0.297  
## 2     1 Mid       218 0.703  
## 3     2 Low        31 0.157  
## 4     2 Mid       154 0.782  
## 5     2 High       11 0.0558 
## 6     2 V. High     1 0.00508
```

---

# Least-cost path: raster as graph, again

Recall Way 3 — a raster can be a graph, with pixels as nodes and neighbors as edges. `gdistance` builds exactly that graph to find a cheapest path.

``` r
cost <- classify(r1, rcl = matrix(c(-Inf,200,20, 200,400,30,
                                        400,900,NA, 900,Inf,40),
                                       ncol=3, byrow=TRUE)) # toy cost surface
A <- c(179706.2, 330570.0)
B <- c(180100.9, 331074.3)
plot(cost)
```

---

``` r
library(raster)     # gdistance predates terra; it wants the legacy class
library(gdistance)
conductance <- transition(raster::raster(cost), function(x) 1/mean(x), 8)
AtoB <- shortestPath(conductance, A, B, output = "SpatialLines")
plot(cost)
lines(AtoB, col="red", lwd=2)
```

Not every package has migrated to `terra` — `raster::raster()` converts back to the legacy class when a downstream package needs it.

---

# Some key things to think about

A few properties of rasters surprise people the first time they work with real data

1. `NA` values are used to define shape

2. Data type determines size

3. Reprojecting a raster is costly

4. Units live in metadata, not in the raster itself

---

# 1. `NA` defines the shape

A raster's *extent* is always a rectangle. Almost nothing real is. `NA` carves an arbitrary shape out of that rectangle — "no data" and "outside the thing I care about" are the same value.

``` r
r1 <- rast(system.file("external/test.grd", package="raster"))
par(mfrow = c(1, 2))
plot(r1, main = "Values (blank = NA)")
plot(!is.na(r1), col = c("grey85", "steelblue"), legend = FALSE,
     main = "Same rectangle, valid cells only")
```

---

# 1. `NA` defines the shape

``` r
ncell(r1)                # the full rectangle
```

```
## [1] 9200
```

``` r
sum(!is.na(values(r1)))  # cells actually "inside" the shape
```

```
## [1] 3178
```

Two-thirds of this rectangle is `NA` — not missing data we wish we had, just bookkeeping. (Recall: `as.data.frame()` earlier silently returned only 3,178 rows, for exactly this reason.)

---

# 2. Data type determines size

A cell is just a number, but *how* that number is stored is a choice you make, not a fact about the data.

.pull-left[

``` r
tmp <- tempfile(fileext = ".tif")
file_size <- function(r, dt) {
  writeRaster(r, tmp, datatype = dt,
              overwrite = TRUE)
  file.size(tmp)
}
tibble::tibble(
  datatype = c("FLT8S", "FLT4S",
               "INT2S", "INT1U (mask)"),
  bytes = c(
    file_size(r1, "FLT8S"),
    file_size(r1, "FLT4S"),
    file_size(round(r1), "INT2S"),
    file_size(!is.na(r1), "INT1U"))
)
```
]

.pull-right[

|datatype     | bytes|
|:------------|-----:|
|FLT8S        | 20740|
|FLT4S        | 17200|
|INT2S        |  7701|
|INT1U (mask) |  1955|

Same values, four sizes. `INT2S` rounds to the nearest meter; the #1 mask only needs 0/1, so `INT1U` fits it alone.
]

---

# 3. Reprojecting a raster is costly

Reprojecting a vector just recomputes vertex coordinates. Reprojecting a raster has to build a **new grid** and resample every cell onto it — there's no shortcut.

``` r
box <- as.polygons(ext(r1), crs = crs(r1))
t_vec  <- system.time(project(box, "EPSG:4326"))[["elapsed"]]
t_rast <- system.time(r1_ll <- project(r1, "EPSG:4326"))[["elapsed"]]

tibble::tibble(what = c("vector: recompute corners", "raster: resample ~9,200 cells"),
               seconds = c(t_vec, t_rast))
```

|what                          | seconds|
|:-----------------------------|-------:|
|vector: recompute corners     |   0.002|
|raster: resample ~9,200 cells |   0.008|

``` r
cat("cell count:", ncell(r1), "->", ncell(r1_ll), "(a *different* grid)")
```

```
## cell count: 9200 -> 9975 (a *different* grid)
```

This toy raster still does measurably more work than its vector counterpart, and cost scales with cell count — for a drone- or satellite-sized raster that gap runs to minutes.

---

# 4. Units live in metadata

A cell holds a bare number. Whether that number means meters, feet, or reflectance is a label someone attached — `terra` stores it, but never checks it.

``` r
units(r1)                 # this file ships with no units at all
```

```
## [1] ""
```

``` r
feet <- r1 * 3.28084       # actually converted to feet...
units(feet) <- "m"             # ...but nothing stops a wrong label
units(feet)
```

```
## [1] "m"
```

``` r
linearUnits(r1)   # meters per coordinate unit — THIS is enforced by the CRS
```

```
## [1] 1
```

The **spatial** unit (what an X/Y coordinate means) is guaranteed by the CRS. The **value** unit (what a cell's number means) is not — it's on you to check and document it.

---

# `terra` &harr; tidyverse cheat sheet

<table class="table" style="font-size: 20px; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> Task </th>
   <th style="text-align:left;"> terra </th>
   <th style="text-align:left;"> tidyverse </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Combine columns/layers </td>
   <td style="text-align:left;"> c(r1, r2) </td>
   <td style="text-align:left;"> bind_cols() </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Select a column/layer </td>
   <td style="text-align:left;"> r[["band"]], r$band </td>
   <td style="text-align:left;"> select(), df$col </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Row-wise arithmetic across columns </td>
   <td style="text-align:left;"> r1 - r2, app(), lapp() </td>
   <td style="text-align:left;"> mutate() </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Recode into categories </td>
   <td style="text-align:left;"> classify() </td>
   <td style="text-align:left;"> case_when(), cut() </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Neighborhood/rolling summary </td>
   <td style="text-align:left;"> focal() </td>
   <td style="text-align:left;"> slider/zoo rolling functions </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Coarsen &amp; summarize groups </td>
   <td style="text-align:left;"> aggregate() </td>
   <td style="text-align:left;"> group_by() %&gt;% summarize() </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Join keys &amp; summarize </td>
   <td style="text-align:left;"> extract() + group_by/summarize </td>
   <td style="text-align:left;"> left_join() + group_by/summarize </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Convert to a plain table </td>
   <td style="text-align:left;"> as.data.frame(r, xy=TRUE) </td>
   <td style="text-align:left;"> (it already is one) </td>
  </tr>
</tbody>
</table>