module12

library(sf)

Warning: package 'sf' was built under R version 4.3.3

Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(ggplot2)

##unzipping file of countries
unzip("countries.zip", exdir = "data/countries")
list.files("data/countries/countries", recursive = TRUE)

[1] "countries.dbf"         "countries.geojson.txt" "countries.shp"        
[4] "countries.shx"        

countries <- st_read(dsn = "./data/countries/countries/countries.shp",
                layer = "countries",
                drivers = "ESRI Shapefile")

options:        ESRI Shapefile 
Reading layer `countries' from data source 
  `C:\Users\bergh\OneDrive\Documents\GEOG5680\module12\data\countries\countries\countries.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 177 features and 64 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
CRS:           NA

##plotting map of median gdp, changing log to make data more readable
ggplot(countries) +
  geom_sf(aes(fill = gdp_md_est)) +
  scale_fill_viridis_c(trans = "log", na.value = "grey90") +
  labs(title = "Median GDP (log-transformed)", fill = "GDP ($M)")

Warning in transformation$transform(x): NaNs produced

Warning in scale_fill_viridis_c(trans = "log", na.value = "grey90"):
log-2.718282 transformation introduced infinite values.

##plotting map of estimated population, changing log to make data more readable
ggplot(countries) +
  geom_sf(aes(fill = pop_est)) +
  scale_fill_viridis_c(trans = "log", na.value = "grey90") +
  labs(title = "Estimated Population (log-transformed)", fill = "Population")

Warning in transformation$transform(x): NaNs produced

Warning in scale_fill_viridis_c(trans = "log", na.value = "grey90"):
log-2.718282 transformation introduced infinite values.

##running code from module to get ndvi data up
library(terra)

Warning: package 'terra' was built under R version 4.3.3

terra 1.8.29

library(sf)
library(dplyr)

Warning: package 'dplyr' was built under R version 4.3.3

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

unzip(zipfile = "./data/rs/rs.zip", exdir = "./data/rs/rs")

Warning in unzip(zipfile = "./data/rs/rs.zip", exdir = "./data/rs/rs"): error 1
in extracting from zip file

list.files("./data/rs/rs", full.names = TRUE)

 [1] "./data/rs/rs/LC08_044034_20170614_B1.tif" 
 [2] "./data/rs/rs/LC08_044034_20170614_B10.tif"
 [3] "./data/rs/rs/LC08_044034_20170614_B11.tif"
 [4] "./data/rs/rs/LC08_044034_20170614_B2.tif" 
 [5] "./data/rs/rs/LC08_044034_20170614_B3.tif" 
 [6] "./data/rs/rs/LC08_044034_20170614_B4.tif" 
 [7] "./data/rs/rs/LC08_044034_20170614_B5.tif" 
 [8] "./data/rs/rs/LC08_044034_20170614_B6.tif" 
 [9] "./data/rs/rs/LC08_044034_20170614_B7.tif" 
[10] "./data/rs/rs/LC08_044034_20170614_B8.tif" 
[11] "./data/rs/rs/LC08_044034_20170614_B9.tif" 

b5 <- rast("./data/rs/rs/LC08_044034_20170614_B5.tif")

unzip("ca_places.zip", exdir = "data/ca_places")
list.files("data/ca_places/ca_places", recursive = TRUE)

[1] "ca_places.dbf" "ca_places.prj" "ca_places.shp" "ca_places.shx"

countries <- st_read(dsn = "./data/ca_places/ca_places/ca_places.shp",
                layer = "ca_places",
                drivers = "ESRI Shapefile")

options:        ESRI Shapefile 
Reading layer `ca_places' from data source 
  `C:\Users\bergh\OneDrive\Documents\GEOG5680\module12\data\ca_places\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

b4 <- rast("./data/rs/rs/LC08_044034_20170614_B4.tif")  

b5 <- rast("./data/rs/rs/LC08_044034_20170614_B5.tif") 

crs(b4); crs(b5)

[1] "PROJCRS[\"WGS 84 / UTM zone 10N\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 10N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-123,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n        AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],\n        BBOX[0,-126,84,-120]],\n    ID[\"EPSG\",32610]]"

[1] "PROJCRS[\"WGS 84 / UTM zone 10N\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 10N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-123,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n        AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],\n        BBOX[0,-126,84,-120]],\n    ID[\"EPSG\",32610]]"

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

##calculating median values for bethel and oakley
ca_places <- st_read("data/ca_places/ca_places/ca_places.shp")

Reading layer `ca_places' from data source 
  `C:\Users\bergh\OneDrive\Documents\GEOG5680\module12\data\ca_places\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

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

bethel <- st_transform(bethel, crs(ndvi))
oakley <- st_transform(oakley, crs(ndvi))

bethel_v <- vect(bethel)
oakley_v <- vect(oakley)

bethel_ndvi <- extract(ndvi, bethel_v)[[1]]
oakley_ndvi <- extract(ndvi, oakley_v)[[1]]

median_bethel <- median(bethel_ndvi, na.rm = TRUE)
median_oakley <- median(oakley_ndvi, na.rm = TRUE)

print(median_bethel)

[1] 1

print(median_oakley)

[1] 1

##calculating median values for bethel and oakley

---