Module 12
AJ Riopelle
2024-06-18
Exercise 1 Make multiple simple maps with new datasets.
- Load in correct packages
library(sf)## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUElibrary(viridis)## Loading required package: viridisLitelibrary(ggplot2)
library(RColorBrewer)
library(terra)## terra 1.7.78- Create graph based on Median GDP
countries = st_read("C:/GEOG 5680/module12/countries/countries.shp")## Reading layer `countries' from data source
## `C:\GEOG 5680\module12\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: NAggplot() + geom_sf(data=countries, aes(fill = gdp_md_est)) +theme_bw()
- Create a graph based on Population
ggplot() + geom_sf(data=countries, aes(fill = pop_est))+ scale_fill_viridis(option = "virdis") +theme_bw()## Warning in viridisLite::viridis(256, alpha, begin, end, direction, option):
## Option 'virdis' does not exist. Defaulting to 'viridis'.
- Create graph based on Income Group
ggplot() + geom_sf(data =countries, aes(fill = income_grp))+theme_bw()
Exercise 2 Using the NDVI values you calculated in the raster section, calculate the median NDVI for Bethel Island and Oakley. - Reading in ca_places and rs raster data
ca_places = st_read("ca_places/ca_places.shp")## Reading layer `ca_places' from data source
## `C:\GEOG 5680\module12\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 84b5 = rast("./rs/LC08_044034_20170614_B5.tif")
b5## class : SpatRaster
## dimensions : 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 : LC08_044034_20170614_B5.tif
## name : LC08_044034_20170614_B5
## min value : 0.0008457669
## max value : 1.0124315023writeRaster(b5, filename="./b5.tif", overwrite = TRUE)
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]]"ca_places = st_transform(ca_places, 32610)
plot(b5, main = "Landsat 8 (B2)")
plot(st_geometry(ca_places), add = TRUE)
-Finding mean NDVI for Bethel Island Extract Bethel Island data from ca_places,b5(crop)
bethel <- ca_places %>%
dplyr::filter(NAME == "Bethel Island")
b5 <- rast("./rs/LC08_044034_20170614_B5.tif")
b5_sub_bethel <- crop(b5, bethel)
plot(b5_sub_bethel)
plot(st_geometry(bethel), add = TRUE)
Mask data from last step
bethel <- ca_places %>%
dplyr::filter(NAME == "Bethel Island")
b5_sub_bethel <- mask(b5_sub_bethel, mask = bethel)
plot(b5_sub_bethel)
plot(st_geometry(bethel), add = TRUE)
Extract Bethel Island from ca_places, b4(crop)
b4 <- rast('./rs/LC08_044034_20170614_B4.tif')
b4_sub_bethel <- crop(b4, bethel)
plot(b4_sub_bethel)
plot(st_geometry(bethel), add = TRUE)
Mask data from last step
bethel <- ca_places %>%
dplyr::filter(NAME == "Bethel Island")
b4_sub_bethel <- mask(b4_sub_bethel, mask = bethel)
plot(b4_sub_bethel)
plot(st_geometry(bethel), add = TRUE)
Summary stats for Bethel Island NDVI
ndvi_bethel <- (b5_sub_bethel - b4_sub_bethel) / (b5_sub_bethel + b4_sub_bethel)
ndvi_bethel## class : SpatRaster
## dimensions : 146, 203, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 616170, 622260, 4207680, 4212060 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## source(s) : memory
## varname : LC08_044034_20170614_B5
## name : LC08_044034_20170614_B5
## min value : -0.6461217
## max value : 0.8458711Median Value for Bethel Island NDVI
global(ndvi_bethel,median)## global
## LC08_044034_20170614_B5 NAPlot the median NDVI for Bethel Island
plot(ndvi_bethel, col=rev(terrain.colors(10)), main = "NDVI for Bethel Island")
-Finding the mean NDVI for Oakley Extract the Oakley data from ca_places, b5(crop)
oakley <- ca_places %>%
dplyr::filter(NAME == "Oakley")
b5 <- rast("./rs/LC08_044034_20170614_B5.tif")
b5_sub_oakley <- crop(b5, oakley)
plot(b5_sub_oakley)
plot(st_geometry(oakley), add = TRUE)
Mask data from last step
oakley <- ca_places %>%
dplyr::filter(NAME == "Oakley")
b5_sub_oakley <- mask(b5_sub_oakley, mask = oakley)
plot(b5_sub_oakley)
plot(st_geometry(oakley), add = TRUE)
Extract the Oakley data from ca_places, b4(crop)
b4 <- rast('./rs/LC08_044034_20170614_B4.tif')
b4_sub_oakley <- crop(b4, oakley)
plot(b4_sub_oakley)
plot(st_geometry(oakley), add = TRUE)
Mask data from last step
oakley <- ca_places %>%
dplyr::filter(NAME == "Oakley")
b4_sub_oakley <- mask(b4_sub_oakley, mask = oakley)
plot(b4_sub_oakley)
plot(st_geometry(oakley), add = TRUE)
Summary stats for Oakley NDVI
ndvi_oakley = (b5_sub_oakley - b4_sub_oakley)/(b5_sub_oakley+b4_sub_oakley)
ndvi_oakley## class : SpatRaster
## dimensions : 192, 392, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 609210, 620970, 4203120, 4208880 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## source(s) : memory
## varname : LC08_044034_20170614_B5
## name : LC08_044034_20170614_B5
## min value : -0.4714017
## max value : 0.8438498Median Value for Oakley NDVI
global(ndvi_oakley,median)## global
## LC08_044034_20170614_B5 NAPlot the Median NDVI for Oakley
plot(ndvi_oakley, col=rev(terrain.colors(10)), main = "NDVI for Oakley")