Module 12

Author

MJ Kemp

Exercise 1

Load in countries.zip file and format properly

library(RColorBrewer)
library(terra)

terra 1.9.27

library(viridis)

Loading required package: viridisLite

library(sf)

Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

library(ggplot2)
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

ca_places = st_read("ca_places/ca_places.shp", quiet = TRUE)
path_to_data = system.file("shape/nc.shp", package = "sf")
countries = st_read("countries/countries.shp")

Reading layer `countries' from data source 
  `C:\Users\u1221910\OneDrive\Desktop\GEOG 5680\Data\Module 12\countries\countries.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 177 features and 64 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
CRS:           NA

countries = st_set_crs(countries, 4326)

Isolate and plot specific socio-economic columns - Population (pop_est) & Income Group (income_gdp)

countries.econ = countries |> select(name, pop_est, income_grp)
plot(countries.econ)

Adjusting color scales & transformations to better interpret data

ggplot() + geom_sf(data=countries.econ, aes(fill = pop_est)) + scale_fill_viridis(option = "viridis") + theme_bw()

ggplot() + geom_sf(data=countries.econ, aes(fill = income_grp)) + theme_bw()

Exercise 2

Read in the NDVI values calculated in the raster section

b2 = rast('rs/LC08_044034_20170614_B2.tif')
b3 = rast('rs/LC08_044034_20170614_B3.tif')
b4 = rast('rs/LC08_044034_20170614_B4.tif')
b5 = rast('rs/LC08_044034_20170614_B5.tif')
stack = c(b5, b4, b3, b2)
ndvi = (b5 - b4) / (b5 + b4)
plot(ndvi, col=rev(terrain.colors(10)), main = "NDVI")

Bethel = ca_places |> dplyr::filter(NAME == "Bethel Island")
Oakley = ca_places |> dplyr::filter(NAME == "Oakley")

Calculate the median NDVI for Bethel Island and Oakley

extract(b5, Bethel, fun='median')

Warning: [extract] transforming vector data to the CRS of the raster

  ID LC08_044034_20170614_B5
1  1               0.3779947

extract(b5, Oakley, fun='median')

Warning: [extract] transforming vector data to the CRS of the raster

  ID LC08_044034_20170614_B5
1  1               0.2823575

b5_sub.BethelOakley = extract(b5, rbind(Bethel, Oakley))

Warning: [extract] transforming vector data to the CRS of the raster

names(b5_sub.BethelOakley) = c("ID", "B5")
ggplot(b5_sub.BethelOakley, aes(x=B5, fill=as.factor(ID))) + geom_histogram(alpha = 0.7, position = 'identity')

`stat_bin()` using `bins = 30`. Pick better value `binwidth`.