R Markdown

Install required libraries if not already installed

if (!requireNamespace(“sf”, quietly = TRUE)) install.packages(“sf”) if (!requireNamespace(“terra”, quietly = TRUE)) install.packages(“terra”) if (!requireNamespace(“dplyr”, quietly = TRUE)) install.packages(“dplyr”)

Load required libraries

library(sf) library(terra) library(dplyr)

Load the raster and shapefile data

b5 <- rast(“./rs/LC08_044034_20170614_B5.tif”) b4 <- rast(“./rs/LC08_044034_20170614_B4.tif”)

Calculate NDVI

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

Load California places shapefile

ca_places <- st_read(“./ca_places/ca_places.shp”)

#Reading layer ca_places' from data source/cloud/project/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

Filter polygons for Bethel Island and Oakley

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

Calculate median NDVI for Bethel Island

bethel_ndvi <- extract(ndvi, bethel, fun = median) bethel_median_ndvi <- bethel_ndvi[, 2]

Calculate median NDVI for Oakley

oakley_ndvi <- extract(ndvi, oakley, fun = median) oakley_median_ndvi <- oakley_ndvi[, 2]

Optional: Extract all NDVI values and perform t-test

bethel_values <- extract(ndvi, bethel)[, 2] oakley_values <- extract(ndvi, oakley)[, 2]

Perform t-test

t_test_result <- t.test(bethel_values, oakley_values) print(t_test_result)

##Welch Two Sample t-test #data: bethel_values and oakley_values #t = 45.979, df = 26353, p-value < 2.2e-16 #alternative hypothesis: true difference in means is not equal to 0 #95 percent confidence interval: # 0.07298078 0.07948010 #sample estimates: #mean of x mean of y #0.4124472 0.3362167