Module 12 - Spatial Data in R

Exercises 1 & 2 – Socio-Economic Mapping and NDVI Analysis

# Load libraries
library(sf)
library(ggplot2)
library(RColorBrewer)
library(viridis)
library(scales)
library(terra)
library(dplyr)

# Exercise 1

# Set working directory for shapefile
setwd("~/Downloads/Intro to R/Module 12/countries")

# Load country shapefile
countries <- st_read("countries.shp", quiet = TRUE)

# GDP map
GDP_plot <- ggplot(countries) +
  geom_sf(aes(fill = gdp_md_est)) +
  scale_fill_viridis_c(option = "plasma", trans = "log10", labels = comma_format()) +
  labs(title = "Median GDP by Country", fill = "GDP (USD)") +
  theme_minimal()

# Population map
POP_plot <- ggplot(countries) +
  geom_sf(aes(fill = pop_est)) +
  scale_fill_viridis_c(option = "cividis", trans = "log10", labels = comma_format()) +
  labs(title = "Population by Country", fill = "Population") +
  theme_minimal()

# Income group map
INCOME_plot <- ggplot(countries) +
  geom_sf(aes(fill = income_grp)) +
  scale_fill_brewer(palette = "Set3", na.value = "gray80") +
  labs(title = "Income Group by Country", fill = "Income Group") +
  theme_minimal()

# Display maps
print(GDP_plot)

print(POP_plot)

print(INCOME_plot)

# Exercise 2

# Set working directory for raster & places
setwd("~/Downloads/Intro to R/Module 12")

# Load and rescale Landsat bands
b5 <- rast("rs/LC08_044034_20170614_B5.tif") / 10000
b4 <- rast("rs/LC08_044034_20170614_B4.tif") / 10000

# Compute NDVI
ndvi <- (b5 - b4) / (b5 + b4)

# Plot NDVI raster
plot(ndvi, col = rev(terrain.colors(10)), main = "NDVI Map")

# Load city shapefile and select Bethel Island & Oakley
places <- st_read("ca_places/ca_places.shp", quiet = TRUE)
bethel <- places %>% filter(NAME == "Bethel Island")
oakley <- places %>% filter(NAME == "Oakley")

# Reproject to match raster
bethel <- st_transform(bethel, crs(ndvi))
oakley <- st_transform(oakley, crs(ndvi))

# Create 5000m buffers
bethel_buf <- st_buffer(bethel, dist = 5000)
oakley_buf <- st_buffer(oakley, dist = 5000)
bethel_buf <- st_transform(bethel_buf, crs(ndvi))
oakley_buf <- st_transform(oakley_buf, crs(ndvi))

# Extract NDVI values from buffers
ndvi_bethel <- extract(ndvi, vect(bethel_buf))[[2]]
ndvi_oakley <- extract(ndvi, vect(oakley_buf))[[2]]

# Remove NA values
ndvi_bethel <- na.omit(ndvi_bethel)
ndvi_oakley <- na.omit(ndvi_oakley)

# Summary stats
message("Median NDVI in 5000m buffer:")
message("  Bethel Island: ", round(median(ndvi_bethel), 4))
message("  Oakley: ", round(median(ndvi_oakley), 4))

# Histograms
par(mfrow = c(1,2))
hist(ndvi_bethel, breaks = 40, main = "Bethel Island NDVI", xlab = "NDVI", col = "darkgreen")
hist(ndvi_oakley, breaks = 40, main = "Oakley NDVI", xlab = "NDVI", col = "steelblue")

par(mfrow = c(1,1))

# T-test comparison
if (sd(ndvi_bethel) > 0 && sd(ndvi_oakley) > 0) {
  message("\nT-test Result (NDVI, 5000m buffer):")
  print(t.test(ndvi_bethel, ndvi_oakley))
  message("The t-test compares whether the average NDVI differs significantly between the buffered areas of Bethel Island and Oakley.")
} else {
  message("\nT-test skipped: One or both areas have no NDVI variation.")
}

## 
##  Welch Two Sample t-test
## 
## data:  ndvi_bethel and ndvi_oakley
## t = 20.173, df = 344180, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.01501996 0.01825272
## sample estimates:
## mean of x mean of y 
## 0.3418600 0.3252237