Geography 13

class: center, middle, inverse, title-slide

# Geography 13
## Lecture 18: Raster Manipulation
### Mike Johnson

---

# Yesterday

- Introduced the raster data structure

- The idea of color scales through RGB channels and color ramps

---

# 1. Bounding Box / Extents

**Geometries** have extents that define the maximum and minimum coverage of the shape in a coordinate reference system

<div class="figure" style="text-align: center">
<img src="lec-img/18-geom-extent.png" alt="Image Source: National Ecological Observatory Network (NEON)" width="496" />
<p class="caption">Image Source: National Ecological Observatory Network (NEON)</p>
</div>

---

# 2. Extent

- When dealing with objects, the extent (or bbox) is derived from the coordinate set

- When dealing with raster data, the extent is a fondational component of the raster data structure
  
--

- That is, we need to know the area the raster is covering!

<div class="figure" style="text-align: center">
<img src="lec-img/17-raster-extent.png" alt="Image Source: National Ecological Observatory Network (NEON)" width="474" />
<p class="caption">Image Source: National Ecological Observatory Network (NEON)</p>
</div>
---

# 3. Discretization

Once we know the **extent**, we need to know _how_ that space is split up

Two complimentary bit of information can tell us this:

- Resolution (res)
   - Number of row and number of columns (nrow/ncol)
  
<div class="figure" style="text-align: center">
<img src="lec-img/17-raster-res.png" alt="Image Source: National Ecological Observatory Network (NEON)" width="500" />
<p class="caption">Image Source: National Ecological Observatory Network (NEON)</p>
</div>
   
---

# So,

A raster is made of an **extent**, and a **resolution** / row-column structure

- A vector of values fill that structure (same way a vector in R can have diminisons)

- These values are often scaled to integers to reduce file size
  
--

- Values are referenced in cartisian space, based on cell index

- A CRS along with the extent, can provide spatial reference / coordinates

---

# ENVI HDR Files

<div class="figure" style="text-align: center">
<img src="lec-img/18-raster.png" alt="Image Source: National Ecological Observatory Network (NEON)" width="474" />
<p class="caption">Image Source: National Ecological Observatory Network (NEON)</p>
</div>
---

# ENVI Data Ingest

---

---

Almost all remote sensing / image analysis begins with the same basic steps:

1. Identifying an area of interest (AOI)
  
  2. Identifying and downloading the relevant images or products
  
  3. Analyzing the raster products

The definition of a AOI is critical because raster data in continuous, therefore we need to define the bounds of the study rather then the bounds of the objects

- **But**, objects often (even typically) define our bounds

---
## Yesterdays Assignment

Find elevation data for Goleta:

- Define the AOI

```r
bb = read_csv("data/uscities.csv") %>%
  st_as_sf(coords = c("lng", "lat"), crs = 4326) %>% 
  filter(city == "Goleta") %>% 
  st_transform(5070) %>% 
  st_buffer(5000) %>% 
  st_bbox() %>% 
  st_as_sfc() %>% 
  st_as_sf()
```

- Read data from elevation map tiles, for a specific zoom, and crop to the AOI

```r
elev = elevatr::get_elev_raster(bb, z = 11) %>% crop(bb)
writeRaster(elev, filename = "data/goleta-elev.tif", overwrite = TRUE)
```

---

- The resulting raster ...

```r
(elev = raster("data/goleta-elev.tif"))
```

```
class      : RasterLayer 
dimensions : 318, 318, 101124  (nrow, ncol, ncell)
resolution : 31.45065, 31.45065  (x, y)
extent     : -2157752, -2147751, 1530395, 1540396  (xmin, xmax, ymin, ymax)
crs        : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
source     : goleta-elev.tif 
names      : goleta.elev 
values     : -70, 501  (min, max)
```

---

# Raster Values

```r
v = values(elev)
class(v)
```

```
[1] "numeric"
```

```r
length(v)
```

```
[1] 101124
```

```r
elev
```

---

# Raster Values

```r
# The length of the vector is equal to the rows * columns
length(v) == nrow(elev) * ncol(elev)
```

```
[1] TRUE
```

```r
# The span of the x extent divided by the resolution equals the raster rows
(xmax(elev) - xmin(elev) ) / res(elev)[1] == nrow(elev) 
```

```
[1] TRUE
```

```r
# The span of the x extent divided by the number of rows equals the raster resolution
(xmax(elev) - xmin(elev) ) / nrow(elev) == res(elev)[1] 
```

```
[1] TRUE
```

---
## All image files are the same!

```r
download.file(url = "https://a.tile.openstreetmap.org/18/43803/104352.png", destfile = "data/104352.png")
img = readPNG("data/104352.png")
class(img)
```

```
[1] "array"
```

```r
typeof(img)
```

```
[1] "double"
```

```r
dim(img)
```

```
[1] 256 256   3
```

---

```r
img[1,1,1:3]
```

```
[1] 1.0000000 0.9607843 0.8980392
```

```r
rgb(1.0000000, 0.9607843, 0.8980392)
```

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

<div class="figure" style="text-align: center">
<img src="lec-img/18-color-picker.png" alt="Google Color Picker" width="50%" />
<p class="caption">Google Color Picker</p>
</div>

---

.pull-left[
![](data/104352.png)
]

.pull-right[
![](lec-img/18-color-picker.png)
]
---

# Raster Algebra

So our raster **data** is stored as a large numeric array/vector

Many generic functions allow for simple algebra on Raster objects,

These include:

- normal algebraic operators such as `+`, `-`, `*`, `/`

- logical operators such as `>`, `>=`, `<`,` ==`, `!`

- functions like `abs`, `round`, `ceiling`, `floor`, `trunc`, `sqrt`, `log`, `log10`, `exp`, `cos`, `sin`, `atan`, `tan`, `max`, `min`, `range`, `prod`, `sum`, `any`, `all`

--
*** 
In these functions you can mix `raster` objects with numbers, as long as the *first* argument is a raster object.

That means you can add 100 to a raster object but not a raster object to 100

```r
# GOOD
raster + 100

# BAD
100 + raster
```

---

### For example:

```r
elev + 100
```

```
class      : RasterLayer 
dimensions : 318, 318, 101124  (nrow, ncol, ncell)
resolution : 31.45065, 31.45065  (x, y)
extent     : -2157752, -2147751, 1530395, 1540396  (xmin, xmax, ymin, ymax)
crs        : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
source     : memory
names      : goleta.elev 
values     : 30, 601  (min, max)
```

```r
log10(elev)
```

```
class      : RasterLayer 
dimensions : 318, 318, 101124  (nrow, ncol, ncell)
resolution : 31.45065, 31.45065  (x, y)
extent     : -2157752, -2147751, 1530395, 1540396  (xmin, xmax, ymin, ymax)
crs        : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
source     : memory
names      : layer 
values     : -Inf, 2.699838  (min, max)
```

---

# Replacement

- Raster values can be replaced on a conditional statements
- Doing this changes the underlying data!
- If you want to retain the original data, you must make a copy of the base layer

.pull-left[

```r
plot(elev)
```

<img src="lecture-18_files/figure-html/unnamed-chunk-21-1.png" width="432" />
]

.pull-right[

```r
*elev2 = elev
*elev2[elev2 <= 0] = NA
plot(elev2)
```

<img src="lecture-18_files/figure-html/unnamed-chunk-22-1.png" width="432" />
]

---

# Modifying a raster

When we want to modify the **extent** of a raster we can _clip_ it to a new bounds

`crop`: lets you reduce the extent of a raster to the extent of another, overlapping object:

.pull-left[

```r
#remotes::install_github("mikejohnson51/AOI")
iv = AOI::aoi_get("Isla Vista") %>% 
  st_transform(crs(elev2))
plot(elev2)
plot(iv, add = TRUE, col = NA)
```

<img src="lecture-18_files/figure-html/unnamed-chunk-23-1.png" width="432" />
]

.pull-right[

```r
*iv_elev = crop(elev2, iv)
plot(iv_elev)
```

<img src="lecture-18_files/figure-html/unnamed-chunk-24-1.png" width="432" />
]

---

# Modifying the underlying data:

`mask`: mask takes an input object (sf, sp, or raster) and set anything not undelying the input to a new value (default = NA)

```r
library(osmdata)
osm = osmdata::opq(st_transform(iv,4326)) %>% 
  add_osm_feature("water", "lagoon") %>% 
  osmdata_sf()

(poly = osm$osm_polygons %>% 
  st_transform(crs(elev)))
```

```
Simple feature collection with 4 features and 8 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -2154933 ymin: 1531819 xmax: -2149682 ymax: 1533700
Projected CRS: NAD83 / Conus Albers
             osm_id     description intermittent natural
23144668   23144668            <NA>         <NA>    <NA>
446941168 446941168            <NA>         <NA>   water
827758986 827758986 Devereux Puddle          yes   water
897309035 897309035            <NA>         <NA>    <NA>
          salt
23144668  <NA>
446941168 <NA>
827758986  yes
897309035 <NA>
                                                                                                       source
23144668                                                                       nhd_import_v0.1_20080127120618
446941168                                                                                                <NA>
827758986 https://www.reddit.com/r/UCSantaBarbara/comments/hpf7s2/shelter_in_place_day_113_we_present_to_you/
897309035                                                                                                <NA>
          tidal  water                       geometry
23144668   <NA>   <NA> POLYGON ((-2152197 1531843,...
446941168  <NA> lagoon POLYGON ((-2149682 1533126,...
827758986   yes lagoon POLYGON ((-2154753 1532808,...
897309035  <NA>   <NA> POLYGON ((-2152663 1532489,...
```

---

.pull-left[

```r
plot(iv_elev)
plot(poly, col = "blue", add  = TRUE)
```

]
.pull-right[

```r
ma = raster::mask(iv_elev, poly)
ma2 = raster::mask(iv_elev, poly, inverse = TRUE)
plot(stack(iv_elev, ma, ma2))
```

<img src="lecture-18_files/figure-html/unnamed-chunk-27-1.png" width="432" />
]

---
# What is `mask` doing?

```r
NA * 7
```

```
[1] NA
```

```r
mask_r = fasterize::fasterize(poly, iv_elev, background = NA)
plot(mask_r)
```

---

```r
base_mask = mask_r * iv_elev
plot(base_mask)
```

---
# Base raster::rasterize is slow...

This means on large rasters, mask is slow..

```r
system.time({
  fasterize::fasterize(poly, elev2, background = NA)
})
```

```
   user  system elapsed 
  0.001   0.001   0.001 
```

```r
system.time({
  raster::rasterize(poly, elev2, background = NA)
})
```

```
   user  system elapsed 
  0.434   0.049   0.485 
```

---

# Crop or/and mask

- Crop is more efficient then mask
- Often you will want to mask and crop a raster
- The correct way to do this is crop _then_ mask

```r
cm = crop(iv_elev, poly) %>%  mask(poly)
plot(cm)
```

<img src="lecture-18_files/figure-html/unnamed-chunk-32-1.png" width="432" />
---

# Aggregate and disaggregate

- `aggregate` and `disaggregate` allow for changing the _resolution_  of a Raster object.

- This is similar to the zoom scaling on a web map except the scale factor is not set to 2

- For aggregate, you need to specify a function determining what to do with the grouped cell values (default = mean).

.pull-left[

```r
plot(iv_elev)
plot(rpoly, add = T)
```

<img src="lecture-18_files/figure-html/unnamed-chunk-34-1.png" width="432" />
]

.pull-right[

```r
*agg = aggregate(iv_elev, 10, fun = max)
plot(agg)
plot(rpoly, add = T)
```

<img src="lecture-18_files/figure-html/unnamed-chunk-35-1.png" width="432" />
]

---

### `calc`

Just like a vector, we can apply functions over a raster with `calc`

These types of formulas are very useful for thresholding analysis

*Question: separate Goleta into the higher and lower elevations*

```r
FUN = function(x){ ifelse(x < mean(x), 1, 2) }
```

```r
*elev3 = calc(elev, FUN)
plot(elev3, col = c("red", "blue"))
```

---

# Summary Values
`cellStats`: computes statistics for the values of each layer in a Raster* object.

```r
cellStats(elev, mean)
```

```
[1] 44.509
```

```r
mean(values(elev), na.rm = TRUE)
```

```
[1] 44.509
```
---

# Why not just `mean()`

In the raster package, functions like `max`, `min`, and `mean`, return a new Raster* object (with a value computed for each cell).

In contrast, `cellStats` returns a single value, computed from the all the values of a layer.

.pull-left[

```r
s = stack(elev, elev^2, elev*.5)
mean(s) %>% plot()
```

<img src="lecture-18_files/figure-html/unnamed-chunk-39-1.png" width="432" />
]

.pull-right[

```r
cellStats(s, mean)
```

```
goleta.elev.1 goleta.elev.2 goleta.elev.3 
      44.5090     8186.4259       22.2545 
```
]

---

### reclassify

A common taks in raster analysis in reclassifying existing values to new values or bins:

```r
(quarts = cellStats(elev2, fivenum))
```

```
[1]   1  13  37 102 501
```

```r
(rcl = data.frame(quarts[1:4], quarts[2:5], 1:4))
```

```
  quarts.1.4. quarts.2.5. X1.4
1           1          13    1
2          13          37    2
3          37         102    3
4         102         501    4
```
---

```r
reclassify(elev2, rcl, include.lowest=TRUE ) %>% 
  plot(col = viridis::viridis(4))
```

---

### Real Example: Rainfall Regions of California

```r
#remotes::install_github("mikejohnson51/climateR")
library(climateR)
AOI = USAboundaries::us_states() %>% 
  dplyr::filter(name == "California")

system.time({ prcp = climateR::getTerraClim(AOI, "prcp", startDate = "2000-01-01", endDate = '2005-12-31') })
```

```
Spherical geometry (s2) switched off
Spherical geometry (s2) switched on
```

```
   user  system elapsed 
  1.829   0.460   6.069 
```
--

```r
prcp$terraclim_prcp
```

```
class      : RasterStack 
dimensions : 229, 247, 56563, 72  (nrow, ncol, ncell, nlayers)
resolution : 0.04166667, 0.04166667  (x, y)
extent     : -124.4167, -114.125, 32.5, 42.04167  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : X2000.01, X2000.02, X2000.03, X2000.04, X2000.05, X2000.06, X2000.07, X2000.08, X2000.09, X2000.10, X2000.11, X2000.12, X2001.01, X2001.02, X2001.03, ... 
min values :      0.0,      0.2,      1.6,      3.5,      0.0,      0.0,      0.0,      0.0,      0.0,      0.0,      2.6,      0.0,      0.0,      3.3,      5.0, ... 
max values :    154.8,    605.2,    597.6,    165.5,    291.1,    124.0,     59.2,     65.7,     99.0,     59.6,    235.2,    151.4,    157.7,    301.4,    319.2, ... 
```

---

.pull-left[
### cellStats

```r
plot(cellStats(prcp$terraclim_prcp, max), type = "l", 
     ylab = "rainfall", xlab = "month since 2000-01")
lines(cellStats(prcp$terraclim_prcp, min), type = "l")
lines(cellStats(prcp$terraclim_prcp, mean), type = "l", col = "darkred", lwd = 2)
```

<img src="lecture-18_files/figure-html/unnamed-chunk-45-1.png" width="432" />
]

.pull-right[
### mean()

```r
plot(mean(prcp$terraclim_prcp), col = blues9)
plot(AOI, add =TRUE, col = NA, lwd = 2)
```

<img src="lecture-18_files/figure-html/unnamed-chunk-46-1.png" width="432" />
]

---

```r
quarts = cellStats(prcp$terraclim_prcp, fivenum)
(quarts = rowMeans(quarts))
```

```
[1]   1.194444   9.252778  19.231944  43.847222 222.505556
```

```r
(rcl = data.frame(quarts[1:4], quarts[2:5], 1:4))
```

```
  quarts.1.4. quarts.2.5. X1.4
1    1.194444    9.252778    1
2    9.252778   19.231944    2
3   19.231944   43.847222    3
4   43.847222  222.505556    4
```
---

```r
reclassify(mean(prcp$terraclim_prcp), 
           rcl, include.lowest=TRUE ) %>% 
   plot(col = blues9)
```

---

# Daily Assignment

In the same Rmd file as yesterday:

- Start with your Goleta elevation raster

- Use `calc` to threshold your raster so that values greater then 0 are equal to 1 and values less (or equal) to zero are equal to NA

- `Multiply` that raster by the original elevation data to isolate the land cells

- `Reclassify` the raster into 6 classes with breaks at even 100m intervals (think of these as topo lines)

- stack your three rasters

- use `setNames()` to change the names of the layers

- plot the stack using the `viridis::viridis` color palette

---

Your final result should look similar to this:

---

# Submission

Submit your knit HTML file to Gauchospace