Geospatial Analysis with R

class: center, middle, inverse, title-slide

.title[
# Geospatial Analysis with R
]
.subtitle[
## Class 8
]

---

```r
library(sf)
library(dplyr)
library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)
data(world.cities, package = "maps")

world <- ne_countries(scale = "medium", returnclass = "sf")
afr_capitals <- world.cities %>% filter(capital == 1) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>% 
  st_intersection(., world %>% filter(continent == "Africa"))
p <- world %>% filter(continent == "Africa") %>% 
  ggplot() + geom_sf(aes(fill = name), lwd = 0.2) + 
  geom_sf(data = afr_capitals, col = "blue", size = 0.5) + 
  scale_fill_grey(guide = FALSE) + theme_minimal()
ggsave(here::here("external/slides/figures/africa_capitals.png"), 
       width = 5, height = 4, dpi = 300, bg = "transparent")
```

---

```r
library(stars)
tif <- system.file("tif/L7_ETMs.tif", package = "stars")
x <- read_stars(tif)
plot(x, axes = TRUE)
plot(x[, , , 4:2], axes = TRUE)

s3url <- glue::glue("/vsis3/africa-fieldboundary-labels/",
                    "tiles_prod/cat2/992155_2023-06/",
                    "image.tif")  # not accessible
p <- read_stars(s3url, proxy = TRUE)
plot(p, rgb = 4:2)
plot(p, rgb = 3:1)
plot(p[, , , 1:3])
```

See [here](https://r-spatial.github.io/stars/articles/stars1.html) for more details.

---
<img src="class8_files/figure-html/unnamed-chunk-3-1.png" style="display: block; margin: auto;" />
---
# Programming tip

- Whenever possible, let the computer do your counting for you:

```r
a <- 1:1001
a[length(a)]
```

```
## [1] 1001
```
- Rather than

```r
a[1001]
```

```
## [1] 1001
```
- Imagine this case:

```r
set.seed(1)
b <- 1:sample(1:10000, 1)
b[length(b)]
```

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

---

## Useful functions
- `print`
- You can print just about anything, but as a generic function, `print` works differently for different types of objects.

```r
a <- "test"
df <- data.frame(numbers = 1:3, abc = letters[1:3])
print(a) # a data frame CAN hold multiple data types. 
```

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

```r
print(df)
```

```
##   numbers abc
## 1       1   a
## 2       2   b
## 3       3   c
```
---

```r
methods(print) ## all the different ways to print...
```

```
##   [1] print.acf*                                          
##   [2] print.AES*                                          
##   [3] print.anova*                                        
##   [4] print.aov*                                          
##   [5] print.aovlist*                                      
##   [6] print.ar*                                           
##   [7] print.Arima*                                        
##   [8] print.arima0*                                       
##   [9] print.AsIs                                          
##  [10] print.aspell*                                       
##  [11] print.aspell_inspect_context*                       
##  [12] print.bibentry*                                     
##  [13] print.Bibtex*                                       
##  [14] print.browseVignettes*                              
##  [15] print.bslib_fragment*                               
##  [16] print.bslib_page*                                   
##  [17] print.by                                            
##  [18] print.cachem*                                       
##  [19] print.changedFiles*                                 
##  [20] print.check_bogus_return*                           
##  [21] print.check_code_usage_in_package*                  
##  [22] print.check_compiled_code*                          
##  [23] print.check_demo_index*                             
##  [24] print.check_depdef*                                 
##  [25] print.check_details*                                
##  [26] print.check_details_changes*                        
##  [27] print.check_doi_db*                                 
##  [28] print.check_dotInternal*                            
##  [29] print.check_make_vars*                              
##  [30] print.check_nonAPI_calls*                           
##  [31] print.check_package_code_assign_to_globalenv*       
##  [32] print.check_package_code_attach*                    
##  [33] print.check_package_code_data_into_globalenv*       
##  [34] print.check_package_code_startup_functions*         
##  [35] print.check_package_code_syntax*                    
##  [36] print.check_package_code_unload_functions*          
##  [37] print.check_package_compact_datasets*               
##  [38] print.check_package_CRAN_incoming*                  
##  [39] print.check_package_datalist*                       
##  [40] print.check_package_datasets*                       
##  [41] print.check_package_depends*                        
##  [42] print.check_package_description*                    
##  [43] print.check_package_description_encoding*           
##  [44] print.check_package_license*                        
##  [45] print.check_packages_in_dir*                        
##  [46] print.check_packages_used*                          
##  [47] print.check_po_files*                               
##  [48] print.check_pragmas*                                
##  [49] print.check_Rd_line_widths*                         
##  [50] print.check_Rd_metadata*                            
##  [51] print.check_Rd_xrefs*                               
##  [52] print.check_RegSym_calls*                           
##  [53] print.check_S3_methods_needing_delayed_registration*
##  [54] print.check_so_symbols*                             
##  [55] print.check_T_and_F*                                
##  [56] print.check_url_db*                                 
##  [57] print.check_vignette_index*                         
##  [58] print.checkDocFiles*                                
##  [59] print.checkDocStyle*                                
##  [60] print.checkFF*                                      
##  [61] print.checkRd*                                      
##  [62] print.checkRdContents*                              
##  [63] print.checkReplaceFuns*                             
##  [64] print.checkS3methods*                               
##  [65] print.checkTnF*                                     
##  [66] print.checkVignettes*                               
##  [67] print.citation*                                     
##  [68] print.cli_ansi_html_style*                          
##  [69] print.cli_ansi_string*                              
##  [70] print.cli_ansi_style*                               
##  [71] print.cli_boxx*                                     
##  [72] print.cli_diff_chr*                                 
##  [73] print.cli_doc*                                      
##  [74] print.cli_progress_demo*                            
##  [75] print.cli_rule*                                     
##  [76] print.cli_sitrep*                                   
##  [77] print.cli_spark*                                    
##  [78] print.cli_spinner*                                  
##  [79] print.cli_tree*                                     
##  [80] print.codoc*                                        
##  [81] print.codocClasses*                                 
##  [82] print.codocData*                                    
##  [83] print.colorConverter*                               
##  [84] print.compactPDF*                                   
##  [85] print.condition                                     
##  [86] print.connection                                    
##  [87] print.CRAN_package_reverse_dependencies_and_views*  
##  [88] print.css*                                          
##  [89] print.data.frame                                    
##  [90] print.Date                                          
##  [91] print.default                                       
##  [92] print.dendrogram*                                   
##  [93] print.density*                                      
##  [94] print.difftime                                      
##  [95] print.dist*                                         
##  [96] print.Dlist                                         
##  [97] print.DLLInfo                                       
##  [98] print.DLLInfoList                                   
##  [99] print.DLLRegisteredRoutines                         
## [100] print.document_context*                             
## [101] print.document_position*                            
## [102] print.document_range*                               
## [103] print.document_selection*                           
## [104] print.dummy_coef*                                   
## [105] print.dummy_coef_list*                              
## [106] print.ecdf*                                         
## [107] print.eigen                                         
## [108] print.factanal*                                     
## [109] print.factor                                        
## [110] print.family*                                       
## [111] print.fileSnapshot*                                 
## [112] print.findLineNumResult*                            
## [113] print.formula*                                      
## [114] print.fseq*                                         
## [115] print.ftable*                                       
## [116] print.function                                      
## [117] print.getAnywhere*                                  
## [118] print.glm*                                          
## [119] print.glue*                                         
## [120] print.hashtab*                                      
## [121] print.hclust*                                       
## [122] print.help_files_with_topic*                        
## [123] print.hexmode                                       
## [124] print.HoltWinters*                                  
## [125] print.hsearch*                                      
## [126] print.hsearch_db*                                   
## [127] print.htest*                                        
## [128] print.html*                                         
## [129] print.html_dependency*                              
## [130] print.htmltools.selector*                           
## [131] print.htmltools.selector.list*                      
## [132] print.infl*                                         
## [133] print.integrate*                                    
## [134] print.isoreg*                                       
## [135] print.json*                                         
## [136] print.key_missing*                                  
## [137] print.kmeans*                                       
## [138] print.knitr_kable*                                  
## [139] print.Latex*                                        
## [140] print.LaTeX*                                        
## [141] print.libraryIQR                                    
## [142] print.lifecycle_warnings*                           
## [143] print.listof                                        
## [144] print.lm*                                           
## [145] print.loadings*                                     
## [146] print.loess*                                        
## [147] print.logLik*                                       
## [148] print.ls_str*                                       
## [149] print.medpolish*                                    
## [150] print.MethodsFunction*                              
## [151] print.mtable*                                       
## [152] print.NativeRoutineList                             
## [153] print.news_db*                                      
## [154] print.nls*                                          
## [155] print.noquote                                       
## [156] print.numeric_version                               
## [157] print.object_size*                                  
## [158] print.octmode                                       
## [159] print.packageDescription*                           
## [160] print.packageInfo                                   
## [161] print.packageIQR*                                   
## [162] print.packageStatus*                                
## [163] print.pairwise.htest*                               
## [164] print.person*                                       
## [165] print.POSIXct                                       
## [166] print.POSIXlt                                       
## [167] print.power.htest*                                  
## [168] print.ppr*                                          
## [169] print.prcomp*                                       
## [170] print.princomp*                                     
## [171] print.proc_time                                     
## [172] print.quosure*                                      
## [173] print.quosures*                                     
## [174] print.R6*                                           
## [175] print.R6ClassGenerator*                             
## [176] print.raster*                                       
## [177] print.Rd*                                           
## [178] print.recordedplot*                                 
## [179] print.restart                                       
## [180] print.RGBcolorConverter*                            
## [181] print.rlang_box_done*                               
## [182] print.rlang_box_splice*                             
## [183] print.rlang_data_pronoun*                           
## [184] print.rlang_dict*                                   
## [185] print.rlang_dyn_array*                              
## [186] print.rlang_envs*                                   
## [187] print.rlang_error*                                  
## [188] print.rlang_fake_data_pronoun*                      
## [189] print.rlang_lambda_function*                        
## [190] print.rlang_message*                                
## [191] print.rlang_trace*                                  
## [192] print.rlang_warning*                                
## [193] print.rlang_zap*                                    
## [194] print.rlang:::list_of_conditions*                   
## [195] print.rle                                           
## [196] print.rlib_bytes*                                   
## [197] print.rlib_error_3_0*                               
## [198] print.rlib_trace_3_0*                               
## [199] print.roman*                                        
## [200] print.sass*                                         
## [201] print.sass_bundle*                                  
## [202] print.sass_layer*                                   
## [203] print.scalar*                                       
## [204] print.sessionInfo*                                  
## [205] print.shiny.tag*                                    
## [206] print.shiny.tag.env*                                
## [207] print.shiny.tag.list*                               
## [208] print.shiny.tag.query*                              
## [209] print.simple.list                                   
## [210] print.smooth.spline*                                
## [211] print.socket*                                       
## [212] print.srcfile                                       
## [213] print.srcref                                        
## [214] print.stepfun*                                      
## [215] print.stl*                                          
## [216] print.stringr_view*                                 
## [217] print.StructTS*                                     
## [218] print.subdir_tests*                                 
## [219] print.summarize_CRAN_check_status*                  
## [220] print.summary.aov*                                  
## [221] print.summary.aovlist*                              
## [222] print.summary.ecdf*                                 
## [223] print.summary.glm*                                  
## [224] print.summary.lm*                                   
## [225] print.summary.loess*                                
## [226] print.summary.manova*                               
## [227] print.summary.nls*                                  
## [228] print.summary.packageStatus*                        
## [229] print.summary.ppr*                                  
## [230] print.summary.prcomp*                               
## [231] print.summary.princomp*                             
## [232] print.summary.table                                 
## [233] print.summary.warnings                              
## [234] print.summaryDefault                                
## [235] print.table                                         
## [236] print.tables_aov*                                   
## [237] print.terms*                                        
## [238] print.ts*                                           
## [239] print.tskernel*                                     
## [240] print.TukeyHSD*                                     
## [241] print.tukeyline*                                    
## [242] print.tukeysmooth*                                  
## [243] print.undoc*                                        
## [244] print.vctrs_bytes*                                  
## [245] print.vctrs_sclr*                                   
## [246] print.vctrs_unspecified*                            
## [247] print.vctrs_vctr*                                   
## [248] print.vignette*                                     
## [249] print.warnings                                      
## [250] print.xfun_raw_string*                              
## [251] print.xfun_rename_seq*                              
## [252] print.xfun_strict_list*                             
## [253] print.xgettext*                                     
## [254] print.xngettext*                                    
## [255] print.xtabs*                                        
## see '?methods' for accessing help and source code
```
---
## Examining objects
- `ls`, `class`, `str`

```r
a <- "test"
df <- data.frame(numbers = 1:3, abc = letters[1:3])
ls() ## list objects in global environemtn
```

```
## [1] "a"  "b"  "d"  "df" "i"
```

```r
class(df) ## class
```

```
## [1] "data.frame"
```

```r
str(df) ## structure
```

```
## 'data.frame':	3 obs. of  2 variables:
##  $ numbers: int  1 2 3
##  $ abc    : chr  "a" "b" "c"
```

---
## Sampling
- Sample without replacement

```r
set.seed(1234) ## set a random seed
v <- 1:100
s <- sample(v, 50)
print(s)
```

```
##  [1] 28 80 22  9  5 38 16  4 86 90 70 79 78 14 56 62 93 84 21 40 92 67 96 66 47
## [26] 81 48  3 41 32 42 43  2 54 49 99 51  6 77 29 71 85 57  8 26 17 58 91 60 76
```

```r
print(sort(s))
```

```
##  [1]  2  3  4  5  6  8  9 14 16 17 21 22 26 28 29 32 38 40 41 42 43 47 48 49 51
## [26] 54 56 57 58 60 62 66 67 70 71 76 77 78 79 80 81 84 85 86 90 91 92 93 96 99
```

- Sample with replacement

```r
s2 <- sample(v, 50, replace = T)
print(sort(s2))
```

```
##  [1]  2  3  3  6  6  9 10 17 17 19 20 22 22 25 27 28 30 32 35 36 41 41 48 50 51
## [26] 55 57 58 58 60 61 63 63 65 66 70 70 71 72 76 77 80 83 85 85 86 86 87 90 96
```

---
## Does S2 have duplicates?

```r
s2_unique <- unique(s2)
print(length(s2))
```

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

```r
print(length(s2_unique))
```

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

```r
print(length(s2) == length(s2_unique))
```

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

---
## In-class exercise

- Create the following:
  - `a`: a random vector of integers with 10 elements drawn from 1-20:
    - Use the `sample` function with `set.seed(10)`
    - Name the elements of `a` with a vector of names starting with "V1" and ending with "V10". 
      - Use the `paste0` function to create those names.
      - Create the identical vector of names using the `paste` function.
  - `b`: Using `a` as an index to select from `letters`
  - `d`: Use `rnorm` with a mean = 100 and an sd of 20
    - Use `?rnorm` to see the help guide.
  - Why did I skip `c`? 
  - Create a list `l` from `a`, `b`, `d`. 
    - Assign the names of the vectors in `l` to the `l`'s elements

---

## 2-d structures

- Create the following:
  - `m`: a matrix with three integer columns named "V1", "V2", "V3"
    - Create each column first as its own vector, then combine
    - `V1` = 1:10
    - `V2` is a random sample between 1:100
    - `V3` is drawn from a random uniform distribution between 0 and 50     - Use the same `set.seed(50)` as before
    - Inspect the `str` and `class` of `m`
  - `dat`, a data.frame built from `V1`, `V2`, `V3`, and `V4`
    - `V4` is a random selection of the letters A-E

---

## Functions
### Components

```r
function_name <- function(arg1, arg2 = 1:10, 
                          arg3 = ifelse(arg2 == 2, TRUE, FALSE)) {
  body
}
```

Three components of a function:
- `formals()`: arguments
- `body()`, the code, which returns the last object generated, unless specified with `return(x)`.
- `environment()`, function finds the values

Unnamed functions are **anonymous** functions. (Used in `*apply`)

---
Using `x` in a function does not change its global value.

```r
x <- 1:10
myfun <- function() {
  x * 10
}
myfun()
```

```
##  [1]  10  20  30  40  50  60  70  80  90 100
```

```r
myfun <- function(x) {
  x <- x * 10
  return(x)
}
x <- 10
myfun(x = 20)
```

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

```r
x
```

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

---
Each time you run `myfun`, a new function environment is created.

```r
myfun <- function(x) {
  x <- x * 10
  print(environment())
  return(x)
}
myfun(x)
```

```
## <environment: 0x7fc664158668>
```

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

```r
myfun(x)
```

```
## <environment: 0x7fc6641d0c48>
```

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

---
## Global assignment.
Use `<<-` to change value of global variable within a function.

```r
a <- 10
myfun <- function(x) {
  a <<- x * 10   ## note <<- instead of <- 
  return(a)
}
myfun(5)
```

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

```r
print(a)
```

```
## [1] 50
```
---
## Useful functions

- `which` finds indices where a condition is true.

```r
v <- 10:15
print(v)
```

```
## [1] 10 11 12 13 14 15
```

```r
a <- which(v %% 3 == 0) ## subset to elements divisible by 3
print(a) ## shows indices where condition is true.
```

```
## [1] 3 6
```

---
## Useful functions
- `which.min` finds index of min value

```r
v <- sample(1:20, 10)
print(v)
```

```
##  [1] 16 11  4  3 10 18  5 20  6 12
```

```r
print(which.min(v)) # index of min value
```

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

```r
print(which.max(v)) # index of max value
```

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

---
## data.frame vs data.table vs. tibble
- all 2D structures.
- data.frame = Base R
- tibble = `tidyverse`
- data.table = fast.

For now, we'll stick to data.frame
---

## data.frame indexing

- `data.frame` uses the following to subset: `[*row conditions, *column conditions]`

```r
df <- data.frame(v1 = 1:5, v2 = 6:10)
rownames(df) <- LETTERS[1:5]
print(df)
```

```
##   v1 v2
## A  1  6
## B  2  7
## C  3  8
## D  4  9
## E  5 10
```

---

## data.frame indexing

- Index using names. 
- Empty index `[  ,  'v2]` means "keep all rows"

```r
df[,'v2'] ## column indexing
```

```
## [1]  6  7  8  9 10
```

```r
df[c("A", "B", "D"), ] ## row indexing
```

```
##   v1 v2
## A  1  6
## B  2  7
## D  4  9
```
---
## data.frame subset

- Logical subset

```r
df[df$v1 > 3, ] ## get observations (rows) where first column is larger than 3
```

```
##   v1 v2
## D  4  9
## E  5 10
```

---
## Control structures
### Branching

- Pay attention to `{ }` placement

```r
a <- 5
if(a > 10) {
  print("Greater than 10!")
} else {
  print("Less than or equal to 10")
}
```

```
## [1] "Less than or equal to 10"
```

---
### Looping

```r
b <- 1:3
for(i in b) print(i)
```

```
## [1] 1
## [1] 2
## [1] 3
```

```r
b <- 1:5
a <- 2
for(i in b){
  a <- 2 * a
  print(a)
}
```

```
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
```
---

### *apply
- A special form of looping
- Intended for *applying* a function to data. Uses *anonymous* function.
- 3 main kinds: `sapply`, `lapply`, `apply`

---
### `sapply`
- `sapply` iterates over input and returns a vector.

```r
v <- 1:10
sapply(v, function(x) x + 10) ## adds 10 to each element in v.
```

```
##  [1] 11 12 13 14 15 16 17 18 19 20
```
Use `{ }` for more complicated functions. BUT be careful with order of `{ }`, `( )`

```r
v1 <- 1:10
v2 <- sapply(v1, function(x){
  y <- x^2 
  return(y)
}) #
print(v2)
```

```
##  [1]   1   4   9  16  25  36  49  64  81 100
```
---
### `sapply`
If you don't specify `return`, the last object created will be returned.

```r
v1 <- 1:10
v2 <- sapply(v1, function(x){
  y <- x^2  ## y will be returned
}) #
print(v2)
```

```
##  [1]   1   4   9  16  25  36  49  64  81 100
```
---

### `lapply`
- Similar to `sapply`, except final object is returned as `list`.
- Useful if you need to store more complex objects (data.frame, plot, raster etc.)

```r
v1 <- 1:10
v2 <- lapply(v1, function(x){
  y <- x^2  ## y will be returned
}) #
print(v2)
```

```
## [[1]]
## [1] 1
## 
## [[2]]
## [1] 4
## 
## [[3]]
## [1] 9
## 
## [[4]]
## [1] 16
## 
## [[5]]
## [1] 25
## 
## [[6]]
## [1] 36
## 
## [[7]]
## [1] 49
## 
## [[8]]
## [1] 64
## 
## [[9]]
## [1] 81
## 
## [[10]]
## [1] 100
```
---

### `apply`
- `apply` works well for 2D data, when you want to apply function over a row or column.

```r
v1 <- sample(1:100, 10)
v2 <- sample(1:100, 10)
DF <- data.frame(v1, v2) ## data frame columns will take names of vectors
DF
```

```
##    v1 v2
## 1   9 38
## 2  43 88
## 3  19 30
## 4  22 70
## 5  46 94
## 6  40 19
## 7  29 79
## 8  57 86
## 9  16 14
## 10 66 23
```
---

Use `apply` to get column max value. The index 2 means "apply function to columns".

```r
colMax <- apply(DF, 2, FUN = max)
colMax
```

```
## v1 v2 
## 66 94
```

---

Use `apply` to get row max value. The index 1 means "apply function to rows".

```r
rowMax <- apply(DF, 1, FUN = max)
rowMax
```

```
##  [1] 38 88 30 70 94 40 79 86 16 66
```

We can use `apply` or `sapply` to create a new column in a data frame.

```r
DF$rowMax <- apply(DF, 1, FUN = max)
DF
```

```
##    v1 v2 rowMax
## 1   9 38     38
## 2  43 88     88
## 3  19 30     30
## 4  22 70     70
## 5  46 94     94
## 6  40 19     40
## 7  29 79     79
## 8  57 86     86
## 9  16 14     16
## 10 66 23     66
```