Module 12 Deliverable - GEOG 5680 Summer 2025

Author

Vivian Strange

Exercises 1 & 2 from Module 12.

GEOG 5680 - 06/17/2025

setwd("~/Desktop/geog5680/Module Deliverables/Module 12")
library(ggplot2)
library(terra)

terra 1.8.54

library(RColorBrewer)
library(sf)

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

library(viridis)

Loading required package: viridisLite

Exercise 1

countries = st_read("countries.shp", quiet = TRUE)

Plot Median GDP: (gdp_md_est)

summary(countries$gdp_md_est)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
     -99    13160    43270   393277   232900 15094000

my_breaks = c(-100, 0, 13160, 43270, 232900, 15094000)
st_crs(4326)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

ggplot() + geom_sf(data=countries, aes(fill = gdp_md_est)) + scale_fill_continuous(name="GDP (millions of USD)",trans = "log", breaks = my_breaks, labels = my_breaks) + labs(x="Longitude (WGS 84)", y="Latitude (WGS 84)", Title="Median GDP By Country")

Warning in transformation$transform(x): NaNs produced

Warning in scale_fill_continuous(name = "GDP (millions of USD)", trans = "log",
: log-2.718282 transformation introduced infinite values.

Warning in transformation$transform(breaks): NaNs produced

Plot Population: (pop_est)

summary(countries$pop_est)

      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
       -99    3441790    9035536   38273988   25946220 1338612970

pop_breaks = c(0, 267722594, 535445188, 803167782, 1070890376, 1338612970)
st_crs(4326)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

ggplot() + geom_sf(data=countries, aes(fill = pop_est)) + scale_fill_viridis(name="Estimated Population", option="viridis", breaks=pop_breaks) + labs(x="Longitude (WGS 84)", y="Latitude (WGS 84)", Title="Estimated Population By Country")

Exercise 2

Calculate the median NDVI for Bethel Island and Oakley.

m12_bethel = st_read("bethel_m12_StrangeV.shp")

Reading layer `bethel_m12_StrangeV' from data source 
  `/Users/vivianstrange/Desktop/geog5680/Module Deliverables/Module 12/bethel_m12_StrangeV.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 16 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 616158.3 ymin: 4207683 xmax: 622256.6 ymax: 4212060
Projected CRS: WGS 84 / UTM zone 10N

st_read("bethel_m12_StrangeV.shp")

Reading layer `bethel_m12_StrangeV' from data source 
  `/Users/vivianstrange/Desktop/geog5680/Module Deliverables/Module 12/bethel_m12_StrangeV.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 16 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 616158.3 ymin: 4207683 xmax: 622256.6 ymax: 4212060
Projected CRS: WGS 84 / UTM zone 10N

m12_oakley = st_read("oakley_m12_StrangeV.shp")

Reading layer `oakley_m12_StrangeV' from data source 
  `/Users/vivianstrange/Desktop/geog5680/Module Deliverables/Module 12/oakley_m12_StrangeV.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 16 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 609198.2 ymin: 4203114 xmax: 620975.2 ymax: 4208891
Projected CRS: WGS 84 / UTM zone 10N

st_read("oakley_m12_StrangeV.shp")

Reading layer `oakley_m12_StrangeV' from data source 
  `/Users/vivianstrange/Desktop/geog5680/Module Deliverables/Module 12/oakley_m12_StrangeV.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 16 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 609198.2 ymin: 4203114 xmax: 620975.2 ymax: 4208891
Projected CRS: WGS 84 / UTM zone 10N

m12_ndvi = rast("ndvi_m12_StrangeV.tif")
rast("ndvi_m12_StrangeV.tif")

class       : SpatRaster 
size        : 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      : ndvi_m12_StrangeV.tif 
name        : LC08_044034_20170614_B5 
min value   :               -0.959918 
max value   :                0.883263

Extract ndvi for each location. Calculate median using fun.

bethel_median_ndvi = extract(m12_ndvi, m12_bethel, fun='median')
extract(m12_ndvi, m12_bethel, fun='median')

  ID LC08_044034_20170614_B5
1  1               0.4180438

The median NDVI for Bethel Island is 0.4180438

oakley_median_ndvi = extract(m12_ndvi, m12_oakley, fun='median')
extract(m12_ndvi, m12_oakley, fun='median')

  ID LC08_044034_20170614_B5
1  1               0.3005387

Exercises 1 & 2 from Module 12.

GEOG 5680 - 06/17/2025

Exercise 1

Plot Median GDP: (gdp_md_est)

Plot Population: (pop_est)

Exercise 2

Calculate the median NDVI for Bethel Island and Oakley.

Extract ndvi for each location. Calculate median using fun.

The median NDVI for Bethel Island is 0.4180438

The median NDVI for Oakley is 0.3005387