exercise12

Exercise 01

0. Setting

library(sf)

Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(ggplot2)
library(viridis)

Loading required package: viridisLite

1. Load the dataset

countries <- st_read("./data/countries/countries.shp", quiet = TRUE)

2. Verify

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"  

class(countries)

[1] "sf"         "data.frame"

3. GDP

ggplot() +
  geom_sf(data = countries, aes(fill = gdp_md_est)) +
  scale_fill_viridis(option = "magma", trans = "log") +
  theme_bw()

Warning in transformation$transform(x): NaNs produced

Warning in scale_fill_gradientn(colours = viridisLite::viridis(256, alpha, :
log-2.718282 transformation introduced infinite values.

4. Population

ggplot() +
  geom_sf(data = countries, aes(fill = pop_est)) +
  scale_fill_viridis(option = "viridis", trans = "log") +
  theme_bw()

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 02

0. Setting

library(sf)
library(terra)

terra 1.9.27

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:terra':

    intersect, union

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

1. Load the dataset

b5 <- rast("./data/rs/LC08_044034_20170614_B5.tif")
b4 <- rast("./data/rs/LC08_044034_20170614_B4.tif")
ca_places <- st_read("./data/ca_places/ca_places.shp", quiet = TRUE)

2. Match CRS

ca_places <- st_transform(ca_places, 32610)

3. NDVI

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

4. Extract Bethel Island and Oakley

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

5. Median NDVI for each place

extract(ndvi, bethel, fun = "median")

  ID LC08_044034_20170614_B5
1  1               0.4180438

extract(ndvi, oakley, fun = "median")

  ID LC08_044034_20170614_B5
1  1               0.3005387

bethel_vals <- extract(ndvi, bethel)
oakley_vals <- extract(ndvi, oakley)

6. Extract only number

bethel_ndvi <- bethel_vals[, 2]
oakley_ndvi <- oakley_vals[, 2]

7. t-test

t.test(bethel_ndvi, oakley_ndvi)

    Welch Two Sample t-test

data:  bethel_ndvi and oakley_ndvi
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