Module 12 Exercise

Part 1

Load all of the important libraries.

library(sf)

Linking to GEOS 3.13.1, GDAL 3.10.2, PROJ 9.5.1; sf_use_s2() is TRUE

library(terra)

terra 1.8.54

library(RColorBrewer)
library(viridis)

Loading required package: viridisLite

library(ggplot2)

Read in the countries file

countries <- st_read("countries/countries.shp", quiet = TRUE)
names(countries)

 [1] "SP_ID"      "scalerank"  "featurecla" "labelrank"  "sovereignt"
 [6] "sov_a3"     "adm0_dif"   "level"      "type"       "admin"     
[11] "adm0_a3"    "geou_dif"   "geounit"    "gu_a3"      "su_dif"    
[16] "subunit"    "su_a3"      "brk_diff"   "name"       "name_long" 
[21] "brk_a3"     "brk_name"   "brk_group"  "abbrev"     "postal"    
[26] "formal_en"  "formal_fr"  "note_adm0"  "note_brk"   "name_sort" 
[31] "name_alt"   "mapcolor7"  "mapcolor8"  "mapcolor9"  "mapcolor13"
[36] "pop_est"    "gdp_md_est" "pop_year"   "lastcensus" "gdp_year"  
[41] "economy"    "income_grp" "wikipedia"  "fips_10"    "iso_a2"    
[46] "iso_a3"     "iso_n3"     "un_a3"      "wb_a2"      "wb_a3"     
[51] "woe_id"     "adm0_a3_is" "adm0_a3_us" "adm0_a3_un" "adm0_a3_wb"
[56] "continent"  "region_un"  "subregion"  "region_wb"  "name_len"  
[61] "long_len"   "abbrev_len" "tiny"       "homepart"   "geometry"  

Look at the range for GDP and create some breaks for the legend

range(countries$gdp_md_est)

[1]      -99 15094000

my_breaks_gdp = c(0, 10000, 100000, 1000000, 10000000)

Plot GDP by country

ggplot() +
  geom_sf(data = countries, aes(fill = gdp_md_est)) +
  scale_fill_viridis(option = "magma", trans = "log",
                     breaks = my_breaks_gdp, labels = my_breaks_gdp) +
  theme_bw()

Warning in transformation$transform(x): NaNs produced

Warning in scale_fill_gradientn(colours = viridisLite::viridis(256, alpha, :
log-2.718282 transformation introduced infinite values.

Look at the range for population and create some breaks for the legend

range(countries$pop_est)

[1]        -99 1338612970

my_breaks_pop = c(0, 100000, 1000000, 10000000, 100000000, 1000000000)

Plot population by country

ggplot() +
  geom_sf(data = countries, aes(fill = pop_est)) +
  scale_fill_viridis(option = "viridis", trans = "log",
                     breaks = my_breaks_pop, labels = my_breaks_pop)

Warning in transformation$transform(x): NaNs produced

Warning in scale_fill_gradientn(colours = viridisLite::viridis(256, alpha, :
log-2.718282 transformation introduced infinite values.

Part 2

Read in the files.

ca_places <- st_read("ca_places/ca_places.shp", quiet = TRUE)
b5 <- rast("rs/LC08_044034_20170614_B5.tif")
b4 <- rast("rs/LC08_044034_20170614_B4.tif")

Reprojecting the CA raster.

ca_places <- st_transform(ca_places, 32610)

Calculating the NDVI

ndvi <- (b5 - b4) / (b5 + b4)

Extracting data for Bethel Island

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:terra':

    intersect, union

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

bethel <- ca_places |>
  filter(NAME == "Bethel Island")
b5_bethel_med <- extract(b5, bethel, fun = 'median')
b4_bethel_med <- extract(b4, bethel, fun = 'median')

Calculate the median NDVI for Bethel Island

medNDVIbethel <- (b5_bethel_med - b4_bethel_med) / (b5_bethel_med + b4_bethel_med)
medNDVIbethel

  ID LC08_044034_20170614_B5
1  0                 0.46711

Extracting the data for Oakley

oakley <- ca_places |>
  filter(NAME == "Oakley")
b5_oak_med <- extract(b5, oakley, fun = 'median')
b4_oak_med <- extract(b4, oakley, fun = 'median')

Calculate the median NDVI for Oakley

medNDVIoak <- (b5_oak_med - b4_oak_med) / (b5_oak_med + b4_oak_med)
medNDVIoak

  ID LC08_044034_20170614_B5
1  0               0.3259331