library(ggplot2)
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
library(terra)
## terra 1.7.78
countries <- st_read("countries.shp")
## Reading layer `countries' from data source 
##   `/Users/rebeccawulff/Desktop/GEOG5680/module12/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

Exercise 1

Make plots of any two of the variables below. Try different color scales and transformations of the data to get the most informative maps.

Graph 1

my_breaks = c(0, 10000, 100000, 1000000)
ggplot() + 
  geom_sf(data = countries, aes(fill = gdp_md_est)) +
  scale_fill_viridis(option = "plasma", trans = "log",
                     breaks = my_breaks, labels = my_breaks) +
  coord_sf(xlim = c(-90,-30), ylim = c(-55,12)) +
  ggtitle("gdp_md_est in South America")
## Warning in transformation$transform(x): NaNs produced
## Warning in scale_fill_gradientn(colours = viridisLite::viridis(256, alpha, :
## log-2.718282 transformation introduced infinite values.

Graph 2

my_breaks <- c(0, 1000000, 10000000, 100000000)
ggplot() + 
  geom_sf(data = countries, aes(fill = pop_est)) +
  scale_fill_viridis(option = "viridis", trans = "log",
                     breaks = my_breaks, labels = my_breaks) +
  coord_sf(xlim = c(-15, 60), ylim = c(-40, 45)) +
  ggtitle("pop_est in Africa") +
  theme_minimal()
## Warning in transformation$transform(x): NaNs produced
## Warning in scale_fill_gradientn(colours = viridisLite::viridis(256, alpha, :
## log-2.718282 transformation introduced infinite values.

Exercise 2

Using the NDVI values you calculated in the raster section, calculate the median NDVI for Bethel Island and Oakley. (Optional - for a slightly more challenging exercise, extract all NDVI values for these two places and compare using a t -test.)

ca_places <- st_read("./ca_places/ca_places.shp")
## Reading layer `ca_places' from data source 
##   `/Users/rebeccawulff/Desktop/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
b4 <- rast('./rs/LC08_044034_20170614_B4.tif')
b5 <- rast('./rs/LC08_044034_20170614_B5.tif')

ndvi = (b5-b4) / (b5+b4)
bethelisland = ca_places %>% 
  dplyr::filter(NAME == "Bethel Island")
bethelisland
## Simple feature collection with 1 feature and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -121.6763 ymin: 38.00883 xmax: -121.6072 ymax: 38.04844
## Geodetic CRS:  WGS 84
##   STATEFP PLACEFP  PLACENS   GEOID          GEOIDFQ          NAME
## 1      06   06210 02407834 0606210 1600000US0606210 Bethel Island
##            NAMELSAD LSAD CLASSFP PCICBSA MTFCC FUNCSTAT    ALAND  AWATER
## 1 Bethel Island CDP   57      U1       N G4210        S 13599384 1468307
##      INTPTLAT     INTPTLON                       geometry
## 1 +38.0281297 -121.6298509 MULTIPOLYGON (((-121.6763 3...
oakley = ca_places %>% 
  dplyr::filter(NAME == "Oakley")
oakley
## Simple feature collection with 1 feature and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -121.7564 ymin: 37.96878 xmax: -121.6226 ymax: 38.02108
## Geodetic CRS:  WGS 84
##   STATEFP PLACEFP  PLACENS   GEOID          GEOIDFQ   NAME    NAMELSAD LSAD
## 1      06   53070 02411294 0653070 1600000US0653070 Oakley Oakley city   25
##   CLASSFP PCICBSA MTFCC FUNCSTAT    ALAND AWATER    INTPTLAT     INTPTLON
## 1      C1       N G4110        A 41091677 789638 +37.9961485 -121.6936235
##                         geometry
## 1 MULTIPOLYGON (((-121.7562 3...
bethelisland_ndvi = extract(ndvi, bethelisland, fun = median, na.rm = TRUE) 
## Warning: [extract] transforming vector data to the CRS of the raster
bethelisland_median_ndvi = bethelisland_ndvi[, 2]
bethelisland_median_ndvi
## [1] 0.4180438
oakley_ndvi = extract(ndvi, oakley, fun = median, na.rm = TRUE)
## Warning: [extract] transforming vector data to the CRS of the raster
oakley_median_ndvi = oakley_ndvi[, 2]
oakley_median_ndvi
## [1] 0.3005387