library(ggplot2)
library(terra)
## terra 1.7.55
library(RColorBrewer)
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(viridis)
## Loading required package: viridisLite

Exercise 1

The zipfile countries.zip is a global shapefile with various socio-economic indicators for different countries. Load this file into R and make plots of any two of the following variables (the variable names are given in brackets).Try different color scales and transformations of the data to get the most informative maps

Read in file

countries <- st_read("./countries/countries.shp")
## Reading layer `countries' from data source 
##   `/Users/benwilliamson/Documents/geog5680/module12/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
plot(st_geometry(countries))

Set the CRS

st_crs(countries)
## Coordinate Reference System: NA
countries <- st_set_crs(countries, 4326)

st_crs(countries)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
st_crs(countries)$epsg
## [1] 4326
class(countries)
## [1] "sf"         "data.frame"

Select the Median GDP column

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
countries2 <- countries %>%
  select(gdp_md_est)

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"
names(countries2)
## [1] "gdp_md_est" "geometry"
plot(countries2)

ggplot Options for Median GDP

ggplot() + 
  geom_sf(data = countries2, aes(fill = gdp_md_est)) +
  theme_bw()

my_breaks = c(0, 10000, 100000, 1000000)
ggplot() + 
  geom_sf(data = countries2, aes(fill = gdp_md_est)) +
  scale_fill_continuous(trans = "log",
                        breaks = my_breaks, labels = my_breaks) +
  theme_bw()
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in discrete y-axis

ggplot() + 
  geom_sf(data = countries2, aes(fill = gdp_md_est)) +
  scale_fill_viridis(option = "inferno", trans = "log",
                     breaks = my_breaks, labels = my_breaks) +
  theme_bw() +
  ggtitle("Median GDP")
## Warning in self$trans$transform(x): NaNs produced
## Warning in self$trans$transform(x): Transformation introduced infinite values
## in discrete y-axis

Select the Population column

library(dplyr)
countries4 <- countries %>%
  select(pop_est)

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"
names(countries4)
## [1] "pop_est"  "geometry"
plot(countries4)

ggplot Options for Income Group

ggplot() + 
  geom_sf(data = countries4, aes(fill = pop_est)) +
  theme_bw()

my_breaks = c(0, 10000, 100000, 1000000)
ggplot() + 
  geom_sf(data = countries4, aes(fill = pop_est)) +
  scale_fill_continuous(trans = "log",
                        breaks = my_breaks, labels = my_breaks) +
  theme_bw()
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in discrete y-axis

ggplot() + 
  geom_sf(data = countries4, aes(fill = pop_est)) +
  scale_fill_viridis(option = "rocket", trans = "log",
                     breaks = my_breaks, labels = my_breaks) +
  theme_bw() +
  ggtitle("Population")
## Warning in self$trans$transform(x): NaNs produced
## Warning in self$trans$transform(x): Transformation introduced infinite values
## in discrete y-axis

Exercise 2

Using the NDVI values you calculated in the raster section, calculate the median NDVI for Bethel Island and Oakley

Reading in ca_places and rs raster data

ca_places <- st_read("./ca_places/ca_places.shp")
## Reading layer `ca_places' from data source 
##   `/Users/benwilliamson/Documents/geog5680/module12/ca_places/ca_places.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1618 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.2695 ymin: 32.53432 xmax: -114.229 ymax: 41.99317
## Geodetic CRS:  WGS 84
b5 <- rast("./rs/LC08_044034_20170614_B5.tif")
b5
## class       : SpatRaster 
## dimensions  : 1245, 1497, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610) 
## source      : LC08_044034_20170614_B5.tif 
## name        : LC08_044034_20170614_B5 
## min value   :            0.0008457669 
## max value   :            1.0124315023
writeRaster(b5, 
            filename = "./b5.tif",
            overwrite = TRUE)
crs(b5)
## [1] "PROJCRS[\"WGS 84 / UTM zone 10N\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 10N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-123,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],\n        BBOX[0,-126,84,-120]],\n    ID[\"EPSG\",32610]]"
ca_places <- st_transform(ca_places, 32610)
plot(b5, main = "Landsat 8 (B2)")
plot(st_geometry(ca_places), add = TRUE)

Finding the mean NDVI for Bethel Island

Extract the Bethel Island data from ca_places

Extracting b5 for Bethel Island

bethel <- ca_places %>% 
  dplyr::filter(NAME == "Bethel Island")
b5 <- rast("./rs/LC08_044034_20170614_B5.tif")
b5_sub_bethel <- crop(b5, bethel)
plot(b5_sub_bethel)
plot(st_geometry(bethel), add = TRUE)

bethel <- ca_places %>% 
  dplyr::filter(NAME == "Bethel Island")
b5_sub_bethel <- mask(b5_sub_bethel, mask = bethel)
plot(b5_sub_bethel)
plot(st_geometry(bethel), add = TRUE)

Extracting b4 for Bethel Island

b4 <- rast('./rs/LC08_044034_20170614_B4.tif')
b4_sub_bethel <- crop(b4, bethel)
plot(b4_sub_bethel)
plot(st_geometry(bethel), add = TRUE)

bethel <- ca_places %>% 
  dplyr::filter(NAME == "Bethel Island")
b4_sub_bethel <- mask(b4_sub_bethel, mask = bethel)
plot(b4_sub_bethel)
plot(st_geometry(bethel), add = TRUE)

Calculating the median NDVI for Bethel Island

ndvi_bethel <- (b5_sub_bethel - b4_sub_bethel) / (b5_sub_bethel + b4_sub_bethel)
ndvi_bethel
## class       : SpatRaster 
## dimensions  : 146, 203, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 616170, 622260, 4207680, 4212060  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610) 
## source(s)   : memory
## varname     : LC08_044034_20170614_B5 
## name        : LC08_044034_20170614_B5 
## min value   :              -0.6461217 
## max value   :               0.8458711

Plot the median NDVI for Bethel Island

plot(ndvi_bethel, col=rev(terrain.colors(10)), main = "NDVI for Bethel Island")

Make a histogram to check for outliers

hist(ndvi_bethel, main = "NDVI for Bethel Island")

Finding the mean NDVI for Oakley

Extract the Oakley data from ca_places

Extracting b5 for Oakley

oakley <- ca_places %>% 
  dplyr::filter(NAME == "Oakley")
b5 <- rast("./rs/LC08_044034_20170614_B5.tif")
b5_sub_oakley <- crop(b5, oakley)
plot(b5_sub_oakley)
plot(st_geometry(oakley), add = TRUE)

oakley <- ca_places %>% 
  dplyr::filter(NAME == "Oakley")
b5_sub_oakley <- mask(b5_sub_oakley, mask = oakley)
plot(b5_sub_oakley)
plot(st_geometry(oakley), add = TRUE)

Extracting b4 for Oakley

b4 <- rast('./rs/LC08_044034_20170614_B4.tif')
b4_sub_oakley <- crop(b4, oakley)
plot(b4_sub_oakley)
plot(st_geometry(oakley), add = TRUE)

oakley <- ca_places %>% 
  dplyr::filter(NAME == "Oakley")
b4_sub_oakley <- mask(b4_sub_oakley, mask = oakley)
plot(b4_sub_oakley)
plot(st_geometry(oakley), add = TRUE)

Calculating the median NDVI for Oakley

ndvi_oakley <- (b5_sub_oakley - b4_sub_oakley) / (b5_sub_oakley + b4_sub_oakley)
ndvi_oakley
## class       : SpatRaster 
## dimensions  : 192, 392, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 609210, 620970, 4203120, 4208880  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610) 
## source(s)   : memory
## varname     : LC08_044034_20170614_B5 
## name        : LC08_044034_20170614_B5 
## min value   :              -0.4714017 
## max value   :               0.8438498

Plot the median NDVI for Oakley

plot(ndvi_oakley, col=rev(terrain.colors(10)), main = "NDVI for Oakley")

Make a histogram to check for outliers

hist(ndvi_oakley, main = "NDVI for Oakley")