Raster Data Models for Geographic Information

name: 1
class: center middle main-title section-title-4

# Raster Data Models for Geographic Information

.class-info[

**Session 7**

.light[HES597: Introduction to Spatial Data in R<br>
Boise State University Human-Environment Systems<br>
Fall 2021]

]

---

name: outline
class: title title-inv-5

# Plan for today

What do we mean by raster data?

Types of raster data

Creating raster data in `R`

Operations on raster data

Functional programming and benchmarking with raster data

---
class: center middle section-title section-title-2 animated fadeIn
# Defining the raster data model
---
<figure>
  <img src="img/04/landscape_geometry.jpeg" alt="ZZZ" title="ZZZ" width="100%">
  <figcaption>Image Source: QGIS User's manual </figcaption>
</figure>

---
name: defining
class: title title-inv-3

# Defining the raster data model
.pull-left[
*   __Vector data__ describe the "exact" locations of features on a landscape (including a Cartesian landscape)

* __Raster data__ represent spatially continuous phenomena (`NA` is possible)

* Depict the alignment of data on a regular lattice (often a square)
  * Operations mimic those for `matrix` objects in `R`

* Geometry is implicit; the spatial extent and number of rows and columns define the cell size
]
.pull-right[
<img src="04-slides_files/figure-html/elevrast-1.png" width="432" style="display: block; margin: auto;" />
]
---
# Working with rasters in `R`

* `raster` - the original workhorse package; built on `sp`, `rgdal`, and `rgeos`
  * `RasterLayer`, `RasterStack`, and `RasterBrick` classes

* `terra` - relatively new; developed by the `raster` folks, but designed to be much faster
   * `SpatRaster` and `SpatVector` classes

* `stars` - developed by `sf` package developers; `tidyverse` compatible; designed for spatio-temporal data
  * `stars` class
  *  Crosswalk between `raster` and `stars` is available [here](https://github.com/r-spatial/stars/wiki/How-%60raster%60-functions-map-to-%60stars%60-functions)
   
---
# Rasters in `R`
.pull-left[

```r
matrix <- matrix(1:16, nrow=4)
matrix
```

```
##      [,1] [,2] [,3] [,4]
## [1,]    1    5    9   13
## [2,]    2    6   10   14
## [3,]    3    7   11   15
## [4,]    4    8   12   16
```

```r
rstr <- raster::raster(matrix)
rstr
```

```
## class      : RasterLayer 
## dimensions : 4, 4, 16  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : 1, 16  (min, max)
```
]
.pull-right[
<img src="04-slides_files/figure-html/mtxrst-1.png" width="360" style="display: block; margin: auto;" />

.small[Note: you must have `raster` or `terra` loaded for `plot()` to work on `Rast*` objects; otherwise you get `Error in as.double(y) : cannot coerce type 'S4' to vector of type 'double'`]
]
---
# Types of raster data
.pull-left[
* __Regular__: constant cell size; axes aligned with Easting and Northing

* __Rotated__: constant cell size; axes not aligned with Easting and Northing

* __Sheared__: constant cell size; axes not parallel

* __Rectilinear__: cell size varies along a dimension

* __Curvilinear__: cell size and orientation dependent on the other dimension

__Note__: `stars` provides suport for rectilinear and curvilinear rasters (not `raster` or `terra`)
]
.pull-right[
<img src="04-slides_files/figure-html/rastype-1.png" width="432" style="display: block; margin: auto;" />
]
---
# Types of raster data

* __Continuous__: numeric data representing a measurement (e.g., elevation, precipitation)

* __Categorical__: integer data representing factors (e.g., land use, land cover)
---
# Continuous raster data
.pull-left[

```r
mintemp <- rast("ftp://ftp.hafro.is/pub/data/rasters/Iceland_minbtemp.tif")
mintemp
```

```
## class       : SpatRaster 
## dimensions  : 340, 375, 1  (nrow, ncol, nlyr)
## resolution  : 2000, 2000  (x, y)
## extent      : -1050226, -300225.9, -699984.3, -19984.32  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=laea +lat_0=69 +lon_0=-4 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## source      : Iceland_minbtemp.tif 
## name        : Iceland_minbtemp 
## min value   :       -0.9982879 
## max value   :         8.603114
```
]
.pull-right[
<img src="04-slides_files/figure-html/plotmin-1.png" width="504" style="display: block; margin: auto;" />
]
---
# Categorical raster data

.pull-left[
* Requires a classification matrix and coercion to factor

* `levels` allows you to define the categories

```r
# Create classification matrix
cm <- matrix(c(
  -2, 2, 0,
  2, 4, 1,
  4, 10, 2), ncol = 3, byrow = TRUE)

# Create a raster with integers
temp_reclass <- classify(mintemp, cm)
tempcats <- c("cold", "mild", "warm")
levels(temp_reclass) <- tempcats
cats(temp_reclass)
```

```
## [[1]]
##   ID category
## 1  0     cold
## 2  1     mild
## 3  2     warm
```
]
.pull-right[

<img src="04-slides_files/figure-html/catplot-1.png" width="504" style="display: block; margin: auto;" />
]
---
# Stacks, bricks, cubes

* When data are aligned in space and/or time, more efficient to process as 'cubes'
---
name: operations
class: center middle section-title section-title-2 animated fadeIn
# Operations on raster data
---
class: title title-inv-3

# Manipulating the geography of raster data

* Projections

* Changing resolutions

* Cropping

* Combining datasets
---
# Projecting raster data
.pull-left[

```r
r <- rast(xmin=-110, xmax=-90, ymin=40, ymax=60, ncols=40, nrows=40)
values(r) <- 1:ncell(r)
plot(r)
```

<img src="04-slides_files/figure-html/rastinit-1.png" width="288" style="display: block; margin: auto;" />
]
.pull-right[
* Transformation from lat/long to planar CRS involves some loss of precision
* New cell values estimated using overlap with original cells
* Interpolation for continuous data, nearest neighbor for categorical data
* Equal-area projections are preferred; especially for large areas
]
---
# Projecting raster data
.pull-left[
* simple method; alignment not guarenteed

```r
newcrs <- "+proj=robin +datum=WGS84"
pr1 <- terra::project(r, newcrs)
plot(pr1)
```

<img src="04-slides_files/figure-html/newproj1-1.png" width="288" style="display: block; margin: auto;" />
]
.pull-right[
* providing a template to ensure alignment

```r
x <- rast(pr1)
# Set the cell size
res(x) <- 200000
pr3 <- terra::project(r, x)
plot(pr3)
```

<img src="04-slides_files/figure-html/newproj2-1.png" width="288" style="display: block; margin: auto;" />
]
---
# Changing resolutions
* `aggregate`, `disaggregate`, `resample` allow changes in cell size
* `aggregate` requires a function (e.g., `mean()` or `min()`) to determine what to do with the grouped values

* `resample` allows changes in cell size __and__ shifting of cell centers (slower)
---
# Changing resolutions
.pull-left[

```r
r <- rast()
values(r) <- 1:ncell(r)
plot(r)
```

<img src="04-slides_files/figure-html/agg1-1.png" width="504" style="display: block; margin: auto;" />
]
.pull-right[

```r
ra <- aggregate(r, 20)
plot(ra)
```

<img src="04-slides_files/figure-html/agg2-1.png" width="504" style="display: block; margin: auto;" />
]
---
# Changing resolutions
.pull-left[

```r
ra <- aggregate(r, 20)
plot(ra)
```

<img src="04-slides_files/figure-html/agg3-1.png" width="504" style="display: block; margin: auto;" />
]
.pull-right[

```r
rd <- disaggregate(r, 20)
plot(rd)
```

<img src="04-slides_files/figure-html/agg4-1.png" width="504" style="display: block; margin: auto;" />
]
---
# Changing Resolutions

```r
r <- rast(nrow=3, ncol=3, xmin=0, xmax=10, ymin=0, ymax=10)
values(r) <- 1:ncell(r)
s <- rast(nrow=25, ncol=30, xmin=1, xmax=11, ymin=-1, ymax=11)
x <- resample(r, s, method="bilinear")
```

<img src="04-slides_files/figure-html/resampplot-1.png" width="504" style="display: block; margin: auto;" />
---
# Cropping

* Useful for subsetting analyses and cartography
* requires an extent object (or an object with one)
* With vectors, need to add `mask()` to crop precisely to the polygon
---
# Cropping
.pull-left[

```r
ext(r)
```

```
## SpatExtent : 0, 10, 0, 10 (xmin, xmax, ymin, ymax)
```

```r
r1 <- crop(r, ext(-50,0,0,30))
```
<img src="04-slides_files/figure-html/plotcrop-1.png" width="504" style="display: block; margin: auto;" />
]

.pull-right[

```r
r2 <- crop(r, ext(-10,50,-20, 10))
opar <- par(no.readonly =TRUE)
```

<img src="04-slides_files/figure-html/plotcrop2-1.png" width="504" style="display: block; margin: auto;" />
]
---
# Combining datasets

* Useful when mosaicing satellite imagery

* Also helpful when using functional programming to utilize computing resources more efficiently (we'll talk about this on Thurs)

```r
m <- merge(r1, r2)
plot(m)
```

<img src="04-slides_files/figure-html/mergex-1.png" width="360" style="display: block; margin: auto;" />
---

# Operations on raster values

* Raster math

* Summarizing and focal functions

* Reclassification

---
# Raster math

* Generally works the same as matrix operations (all layers must be aligned)

```r
r <- rast(ncol=5, nrow=5)
values(r) <- 1:ncell(r)
r2 <- r*2
r3 <- r + r2
```

<img src="04-slides_files/figure-html/rastmathplot-1.png" width="504" style="display: block; margin: auto;" />
---
# Focal Functions and Summmaries

* Allow 'moving window' summaries of raster objects
.pull-left[

```r
r <- rast(ncols=10, nrows=10, ext(0, 10, 0, 10))
values(r) <- 1:ncell(r)
f <- focal(r, w=3, fun="mean") 
f2 <- focal(r, w=3, fun="sum") 
```
]
.pull-right[
<img src="04-slides_files/figure-html/plotfocal-1.png" width="504" style="display: block; margin: auto;" />
]

---
# Reclassifying

.pull-left[

```r
r <- rast(ncol=3, nrow=2)
values(r) <- 1:ncell(r)
values(r)
```

```
##      lyr.1
## [1,]     1
## [2,]     2
## [3,]     3
## [4,]     4
## [5,]     5
## [6,]     6
```
]
.pull-right[

```r
s <- app(r, fun=function(x){ x[x < 4] <- NA; return(x)} )
as.matrix(s)
```

```
##      lyr.1
## [1,]    NA
## [2,]    NA
## [3,]    NA
## [4,]     4
## [5,]     5
## [6,]     6
```
]
---
# Reclassifying

* `app` function works like `apply` 
* beginning of functional programming
* functional programming is useful when the same operation is applied repeatedly (and independently)

---

# Converting vectors to rasters

* Sometimes necessary to convert between data models

* `raster::rasterize`, `terra::rasterize`, `stars::st_rasterize`, and `fasterize::fasterize` all will convert polygons to raster data

* `stars::st_polygonize` will work in the opposite direction
---
# A note about support

* We talked briefly about the `agr` option in the `st_sf()` function

* `agr` refers to the attribute-geometry-relationship which can be:
  * __constant__ = applies to every point in the geometry (lines and polygons are just lots of points)
  * __identity__ = a value unique to a geometry
  * __aggregate__ = a single value that integrates data across the geometry

* __Support__ is the area to which an attribute applies.

* Rasters can have __point__ (attribute refers to the cell center) or __cell__ (attribute refers to an area similar to the pixel) support
---