R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot. # Module 12 # Ellie Spencer # Intro to R

setwd(“C:/Users/u1317304/OneDrive - University of Utah/GEOG 5680/countries”)

to begin, here is some pre-exercise practice:

path_to_data <- system.file(“shape/nc.shp”, package=“sf”)

north_carolina <- st_read(path_to_data, quiet = TRUE)

north_carolina <- north_carolina[ , c(“CNTY_ID”, “NAME”, “AREA”, “PERIMETER”)]

north_carolina

point_one <- st_point(c(0, 3))

point_two <- st_point(c(5, 7))

a_line <- st_linestring(c(point_one, point_two))

point_one

a_line

st_geometry_type(a_line)

st_geometry_type(north_carolina[1, ])

st_crs(4326)

st_crs(4326)$proj4string

st_bbox(north_carolina)

utah <- st_read(dsn = “./utahcounty/utahcounty.shp”, layer = “utahcounty”, drivers = “ESRI Shapefile”)

utah <- st_read(“./utahcounty/utahcounty.shp”) plot(st_geometry(utah)) wna_climate <- read.csv(“./WNAclimate.csv”) head(wna_climate) wna_climate <- st_as_sf(wna_climate, coords = c(“LONDD”, “LATDD”), crs = 4326)

wna_climate plot(st_geometry(wna_climate), pch = 19, col = alpha(“darkgreen”, 0.5)) st_write(obj = wna_climate, dsn = “./wnaclim.shp”, layer = “wnaclim”, drivers = “ESRI Shapefile”) st_write(wna_climate, dsn = “./wnaclim.shp”) st_crs(utah) utah <- st_set_crs(utah, 4326)

st_crs(utah) <- 4326 utah <- st_transform(utah, crs = 32612)

st_crs(utah) st_crs(utah)\(epsg st_crs(wna_climate)\)epsg format(st_crs(utah)) class(utah) utah

utah2 <- utah[ , c(“NAME”, “FIPS”, “POP_CURRES”)]

method 2 (dplyr)

library(dplyr) utah2 <- utah %>% select(NAME, FIPS, POP_CURRES)

names(utah) names(utah2) extract(b5, cbind(615000,4199000))

EXERCISE ONE

library(sf) library(ggplot2) library(terra) library(RColorBrewer) library(viridis)

shapefile_path <- (“C:317304- University of Utah.shp”) countries_sf <-st_read(shapefile_path) countries_terra <- vect(“C:317304- University of Utah.shp”) countries <- st_as_sf(countries_terra)

head(countries)

Median GDP

ggplot(data = countries) + geom_sf(aes(fill = gdp_md_est)) + scale_fill_viridis(option = “plasma”, trans = “log10”) + theme_minimal() + labs(title = “Median GDP of Countries”, fill = “GDP (log scale)”)

Income Group

ggplot(data = countries) + geom_sf(aes(fill = income_grp)) + scale_fill_brewer(palette = “Set3”) + theme_minimal() + labs(title = “Income Group of Countries”, fill = “Income Group”)

EXERCISE TWO

setwd(“C:/Users/u1317304/OneDrive - University of Utah/GEOG 5680/countries”)

library(raster) library(terra) library(sf) library(ggplot2) library(RColorBrewer) library(dplyr)

Load the raster data (replace with actual file paths)

b4_path <- “path_to_band4.tif” b5_path <- “path_to_band5.tif” b4 <- raster(b4_path) b5 <- raster(b5_path)

Calculate NDVI

ndvi <- (b5 - b4) / (b5 + b4)

Plot NDVI

plot(ndvi, col=rev(terrain.colors(10)), main = “NDVI”) hist(ndvi, main = “NDVI”)

Threshold to identify water areas (NDVI < 0)

water <- ndvi < 0 plot(water)

Filter Bethel Island and Oakley from ca_places

bethel <- ca_places %>% dplyr::filter(NAME == “Bethel Island”) oakley <- ca_places %>% dplyr::filter(NAME == “Oakley”)

Convert sf objects to Spatial objects for raster package compatibility

bethel_sp <- as(bethel, “Spatial”) oakley_sp <- as(oakley, “Spatial”)

Mask and crop NDVI raster with Bethel Island shapefile

ndvi_bethel <- mask(ndvi, bethel_sp) ndvi_bethel <- crop(ndvi_bethel, bethel_sp)

Mask and crop NDVI raster with Oakley shapefile

ndvi_oakley <- mask(ndvi, oakley_sp) ndvi_oakley <- crop(ndvi_oakley, oakley_sp)

Calculate median NDVI for Bethel Island

median_ndvi_bethel <- cellStats(ndvi_bethel, stat = ‘median’, na.rm = TRUE) cat(“Median NDVI for Bethel Island:”, median_ndvi_bethel, “”)

Calculate median NDVI for Oakley

median_ndvi_oakley <- cellStats(ndvi_oakley, stat = ‘median’, na.rm = TRUE) cat(“Median NDVI for Oakley:”, median_ndvi_oakley, “”)

Plot NDVI for Bethel Island

plot(ndvi_bethel, col=rev(terrain.colors(10)), main = “NDVI - Bethel Island”) plot(st_geometry(bethel), add = TRUE)

Plot NDVI for Oakley

plot(ndvi_oakley, col=rev(terrain.colors(10)), main = “NDVI - Oakley”) plot(st_geometry(oakley), add = TRUE)

Extract NDVI values for random points in Bethel Island

random_pnts_bethel <- st_sample(bethel, size = 20) ndvi_values_bethel <- extract(ndvi, st_coordinates(random_pnts_bethel))

Extract NDVI values for random points in Oakley

random_pnts_oakley <- st_sample(oakley, size = 20) ndvi_values_oakley <- extract(ndvi, st_coordinates(random_pnts_oakley))

Combine and plot histogram of NDVI values for Bethel Island and Oakley

ndvi_combined <- data.frame( ID = rep(c(“Bethel Island”, “Oakley”), each = 20), NDVI = c(ndvi_values_bethel, ndvi_values_oakley) )

ggplot(ndvi_combined, aes(x = NDVI, fill = ID)) + geom_histogram(alpha = 0.7, position = ‘identity’, bins = 30) + labs(title = “Histogram of NDVI Values”, x = “NDVI”, fill = “Location”) + theme_minimal()