library(sf)
library(terra)
library(dplyr)
library(ggplot2)
library(viridis)
library(RColorBrewer)
module12 <- "C:/Users/yolom/OneDrive/Desktop/GEOG5680/Module 12"GEOG 5680 Module 12 Exercises
Exercise 1
countries_zip <- normalizePath(
file.path(module12, "countries.zip"),
winslash = "/",
mustWork = TRUE
)
countries <- st_read(
paste0("/vsizip/", countries_zip, "/countries/countries.shp"),
quiet = TRUE
)
countries <- countries %>%
mutate(
pop_plot = ifelse(pop_est > 0, pop_est, NA),
gdp_plot = ifelse(gdp_md_est > 0, gdp_md_est, NA)
)
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"
[66] "pop_plot" "gdp_plot" Population
ggplot() +
geom_sf(
data = countries,
aes(fill = pop_plot),
color = "gray40",
linewidth = 0.1
) +
scale_fill_viridis(
option = "mako",
trans = "log10",
na.value = "gray85",
name = "Population"
) +
theme_bw() +
labs(title = "World Population by Country")
Median GDP
ggplot() +
geom_sf(
data = countries,
aes(fill = gdp_plot),
color = "gray40",
linewidth = 0.1
) +
scale_fill_viridis(
option = "plasma",
trans = "log10",
na.value = "gray85",
name = "GDP"
) +
theme_bw() +
labs(title = "Median GDP by Country")
Optional: Income Group
ggplot() +
geom_sf(
data = countries,
aes(fill = income_grp),
color = "gray40",
linewidth = 0.1
) +
scale_fill_brewer(
palette = "Set2",
na.value = "gray85",
name = "Income group"
) +
theme_bw() +
labs(title = "Income Group by Country")
Exercise 2
rs_zip <- normalizePath(
file.path(module12, "rs.zip"),
winslash = "/",
mustWork = TRUE
)
ca_places_zip <- normalizePath(
file.path(module12, "ca_places.zip"),
winslash = "/",
mustWork = TRUE
)
ca_places <- st_read(
paste0("/vsizip/", ca_places_zip, "/ca_places/ca_places.shp"),
quiet = TRUE
)
b4 <- rast(paste0("/vsizip/", rs_zip, "/rs/LC08_044034_20170614_B4.tif"))
b5 <- rast(paste0("/vsizip/", rs_zip, "/rs/LC08_044034_20170614_B5.tif"))Calculate NDVI
ndvi <- (b5 - b4) / (b5 + b4)
names(ndvi) <- "NDVI"
ndviclass : SpatRaster
size : 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(s) : memory
varname : LC08_044034_20170614_B5
name : NDVI
min value : -0.959918
max value : 0.883263plot(ndvi, col = rev(terrain.colors(10)), main = "NDVI")
Extract Bethel Island and Oakley
ca_places_utm <- st_transform(ca_places, crs(b5))
bethel <- ca_places_utm %>%
filter(NAME == "Bethel Island")
oakley <- ca_places_utm %>%
filter(NAME == "Oakley")
places <- rbind(bethel, oakley)
placesSimple feature collection with 2 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 609198.2 ymin: 4203114 xmax: 622256.6 ymax: 4212060
Projected CRS: WGS 84 / UTM zone 10N
STATEFP PLACEFP PLACENS GEOID GEOIDFQ NAME
1 06 06210 02407834 0606210 1600000US0606210 Bethel Island
2 06 53070 02411294 0653070 1600000US0653070 Oakley
NAMELSAD LSAD CLASSFP PCICBSA MTFCC FUNCSTAT ALAND AWATER
1 Bethel Island CDP 57 U1 N G4210 S 13599384 1468307
2 Oakley city 25 C1 N G4110 A 41091677 789638
INTPTLAT INTPTLON geometry
1 +38.0281297 -121.6298509 MULTIPOLYGON (((616158.3 42...
2 +37.9961485 -121.6936235 MULTIPOLYGON (((609214.1 42...Median NDVI
median_ndvi <- extract(
ndvi,
vect(places),
fun = median,
na.rm = TRUE
)
median_ndvi$Place <- places$NAME
median_ndvi ID NDVI Place
1 1 0.4180438 Bethel Island
2 2 0.3005387 OakleyNDVI Values by Place
ndvi_values <- extract(ndvi, vect(places))
names(ndvi_values) <- c("ID", "NDVI")
ndvi_values$Place <- places$NAME[ndvi_values$ID]
ndvi_values <- ndvi_values[!is.na(ndvi_values$NDVI), ]
head(ndvi_values) ID NDVI Place
1 1 0.101607996 Bethel Island
2 1 -0.004937939 Bethel Island
3 1 -0.176106054 Bethel Island
4 1 0.228816837 Bethel Island
5 1 0.417747920 Bethel Island
6 1 0.413580546 Bethel Islandggplot(ndvi_values, aes(x = NDVI, fill = Place)) +
geom_histogram(alpha = 0.6, position = "identity", bins = 30) +
theme_bw() +
labs(title = "NDVI Values for Bethel Island and Oakley")
Optional T-Test
t.test(NDVI ~ Place, data = ndvi_values)
Welch Two Sample t-test
data: NDVI by Place
t = 45.979, df = 26353, p-value < 2.2e-16
alternative hypothesis: true difference in means between group Bethel Island and group Oakley is not equal to 0
95 percent confidence interval:
0.07298078 0.07948010
sample estimates:
mean in group Bethel Island mean in group Oakley
0.4124472 0.3362167 Summary
The maps in Exercise 1 show global differences in population, GDP, and income group. The NDVI analysis in Exercise 2 compares vegetation greenness between Bethel Island and Oakley using Landsat 8 imagery.