Research Questions
1.) Is there a correlation between CO2 output from facilities with agricultural production throughout the US.
Map of the US, set buffers around them and find agricultural producers within that buffer. Look for correlations in their production (see if they drop/rise at the same times, spatial autocorrelation, etc.)
2.) What is causing the correlation in the data? (CASE STUDIES)
Isolate the top 5 cells with the greatest correlation between c02 and agricultural production.
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
## Warning in CPL_gdal_init(): GDAL Error 1: libcfitsio.so.10: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libcfitsio.so.10: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libpodofo.so.2: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libpodofo.so.2: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libarrow.so.1500: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libarrow.so.1500: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libmariadb.so.3: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libmariadb.so.3: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libarrow_dataset.so.1500: cannot open
## shared object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libarrow_dataset.so.1500: cannot open
## shared object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libcfitsio.so.10: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libcfitsio.so.10: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libpodofo.so.2: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libpodofo.so.2: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libarrow.so.1500: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libarrow.so.1500: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libmariadb.so.3: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libmariadb.so.3: cannot open shared
## object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libarrow_dataset.so.1500: cannot open
## shared object file: No such file or directory
## Warning in CPL_gdal_init(): GDAL Error 1: libarrow_dataset.so.1500: cannot open
## shared object file: No such file or directory
## Linking to GEOS 3.12.0, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE
library(dplyr)
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(spData)
library(ggplot2)
library(ggthemes)
library(spdep)
library(spatialreg)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
##
## Attaching package: 'spatialreg'
##
## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption
library(GWmodel)
## Loading required package: robustbase
## Loading required package: sp
## Loading required package: Rcpp
## Welcome to GWmodel version 2.3-2.
library(tidyr)
library(terra)
## Warning: libcfitsio.so.10: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libcfitsio.so.10: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libpodofo.so.2: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libpodofo.so.2: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libarrow.so.1500: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libarrow.so.1500: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libmariadb.so.3: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libmariadb.so.3: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libarrow_dataset.so.1500: cannot open shared object file: No such file
## or directory (GDAL error 1)
## Warning: libarrow_dataset.so.1500: cannot open shared object file: No such file
## or directory (GDAL error 1)
## Warning: libcfitsio.so.10: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libcfitsio.so.10: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libpodofo.so.2: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libpodofo.so.2: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libarrow.so.1500: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libarrow.so.1500: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libmariadb.so.3: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libmariadb.so.3: cannot open shared object file: No such file or
## directory (GDAL error 1)
## Warning: libarrow_dataset.so.1500: cannot open shared object file: No such file
## or directory (GDAL error 1)
## Warning: libarrow_dataset.so.1500: cannot open shared object file: No such file
## or directory (GDAL error 1)
## terra 1.7.71
##
## Attaching package: 'terra'
##
## The following object is masked from 'package:tidyr':
##
## extract
library(leaflet)
library(urbnmapr)
library(cowplot)
##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:ggthemes':
##
## theme_map
##
## The following object is masked from 'package:lubridate':
##
## stamp
library(usdata) #https://rdrr.io/cran/usdata/man/abbr2state.html
library(magrittr)
##
## Attaching package: 'magrittr'
##
## The following objects are masked from 'package:terra':
##
## extract, inset
##
## The following object is masked from 'package:purrr':
##
## set_names
##
## The following object is masked from 'package:tidyr':
##
## extract
setwd("/home/luc/Documents/R Class/Final Project")
raw_facilities_data <- read.csv("co2_data.csv") %>%
filter(!(state %in% c("AK", "HI", "PR", "VI", "GU", "MP", "AS")))
facilities_data <- read.csv("co2_data.csv") %>%
filter(!(state %in% c("AK", "HI", "PR", "VI", "GU", "MP", "AS"))) %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
st_transform("EPSG:4326")
counties <- get_urbn_map("counties")
counties.sf <- get_urbn_map("counties", sf=TRUE) %>%
st_transform("EPSG:4326")
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## Warning in CPL_crs_from_input(x): GDAL Message 1: CRS EPSG:2163 is deprecated.
## Its non-deprecated replacement EPSG:9311 will be used instead. To use the
## original CRS, set the OSR_USE_NON_DEPRECATED configuration option to NO.
states.sf <- get_urbn_map("states", sf=TRUE) %>%
st_transform("EPSG:4326")
## old-style crs object detected; please recreate object with a recent sf::st_crs()
state_emissions <- st_join(states.sf, facilities_data) %>%
group_by(state) %>%
summarize(total_emissions = sum(co2_emissions, na.rm = TRUE)) %>%
na.omit()
county_emissions <- st_join(counties.sf, facilities_data) %>%
group_by(county_name) %>%
summarize(total_emissions = sum(co2_emissions, na.rm = TRUE),
total_facs = n_distinct(facility_id)) %>%
na.omit()
crop_sales_data <- read.csv("crop_sales_yearly.csv") %>%
group_by(State) %>%
summarize(total_sales = sum(as.numeric(Value), na.rm = TRUE)) %>%
na.omit() %>%
mutate(State = state2abbr(State))
## Warning: There were 27 warnings in `summarize()`.
## The first warning was:
## ℹ In argument: `total_sales = sum(as.numeric(Value), na.rm = TRUE)`.
## ℹ In group 3: `State = "ARIZONA"`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 26 remaining warnings.
state_facilities_co2_emissions_corr <- st_join(state_emissions, facilities_data) %>%
group_by(state.y) %>%
mutate(total_facilities = n_distinct(facility_id)) %>%
na.omit()
soil_health <- rast("soil_health.tif") %>%
project("EPSG:4326") #https://datacarpentry.org/r-raster-vector-geospatial/03-raster-reproject-in-r.html
county_sac <- merge(counties, county_emissions, by.x = "county_name", by.y = "county_name") %>%
group_by(county_name) %>%
merge(crop_sales_data, by.x = "state_abbv", by.y = "State") %>%
na.omit()
county_sac.sf <- st_as_sf(county_sac, coords = c("long", "lat"), crs = 4326) %>%
st_transform("EPSG:4326")
unique_breaks <- unique(county_emissions$total_emissions)
pal_fun <- colorQuantile("YlOrRd", unique_breaks, n = 9, na.color = NA)
pal_fun2 <- colorQuantile("YlOrRd", NULL, n = 9, na.color = NA)
pop_up <- paste("Facility Name: ", facilities_data$facility_name, "<br>",
"CO2 Emissions: ", facilities_data$co2_emissions, "<br>",
"State: ", facilities_data$state, "<br>",
"State Emmissions: ", counties.sf$total_emissions, "<br>",
"County: ", counties.sf$county_name, "<br>",
"City: ", facilities_data$city, "<br>")
leaflet() %>%
leaflet::addPolygons(data = counties.sf,
stroke = FALSE,
#color = "black",
#weight = .5,
fillColor = ~pal_fun(county_emissions$total_emissions),
fillOpacity = 1,
smoothFactor = 0.5,
group = "CO2 Emissions") %>%
leaflet::addRasterImage(soil_health,
group = "Soil Health",
colors = pal_fun2) %>%
leaflet::addCircleMarkers(data = facilities_data,
radius = 1,
stroke = FALSE,
fillOpacity = 1.5,
popup = pop_up,
group = "Facilities") %>%
addLegend("bottomleft",
pal=pal_fun,
values=county_emissions$total_emissions,
title = 'CO2 Emissions',
group = "CO2 Emissions") %>%
addLegend("bottomleft",
pal=pal_fun2,
values=values(soil_health),
title = "Soil Health",
group = "Soil Health") %>%
addTiles(group = "OpenStreetMap") %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addLayersControl(baseGroups = c("OpenStreetMap", "CO2 Emissions", "Soil Health"),
overlayGroups = c("Facilities"),
options = layersControlOptions(collapsed = FALSE))
#Bar graph of top 10 states' co2 levels
top_state_emissions <- state_emissions %>%
arrange(desc(total_emissions)) %>%
head(10) %>%
ggplot(aes(x = reorder(state, total_emissions), y = total_emissions)) +
geom_bar(stat = "identity", fill = "blue") +
coord_flip() +
labs(title = "Top 10 States by CO2 Emissions",
x = "State",
y = "CO2 Emissions (tons)") +
theme_minimal()
top_crop_sales <- crop_sales_data %>%
arrange(desc(total_sales)) %>%
head(10) %>%
ggplot(aes(x = reorder(State, total_sales), y = total_sales)) +
geom_bar(stat = "identity", fill = "green") +
coord_flip() +
labs(title = "Top 10 States by Crop Sales",
x = "State",
y = "Total Sales (USD)") +
theme_minimal()
state_emissions_crop_sales_corr <- merge(state_emissions, crop_sales_data, by.x = "state", by.y = "State") %>%
na.omit()
state_facilities_co2_emissions_corr_plot <- ggplot(state_facilities_co2_emissions_corr, aes(x = total_facilities, y = total_emissions)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Number of Facilities vs. CO2 Emissions",
x = "Number of Facilities",
y = "CO2 Emissions (tons)") +
theme_minimal()
state_emissions_crop_sales_corr_plot <- ggplot(state_emissions_crop_sales_corr, aes(x = total_emissions, y = total_sales)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "CO2 Emissions vs. Crop Sales",
x = "CO2 Emissions (tons)",
y = "Total Sales (USD)") +
theme_minimal()
top_plots <- plot_grid(top_state_emissions, top_crop_sales, state_facilities_co2_emissions_corr_plot, state_emissions_crop_sales_corr_plot, ncol = 2) %>% print()
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
# lm_lag <- errorsarlm(total_emissions ~ total_sales + total_facs,
# data = county_sac,
# listw = county_sac_list2,
# zero.policy = TRUE,
# na.action = na.omit)
# summary(lm_lag)
summary(lm1 <- lm(total_emissions ~ total_sales + total_facs, data = county_sac))
##
## Call:
## lm(formula = total_emissions ~ total_sales + total_facs, data = county_sac)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15655364 -904106 844666 892526 39479693
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.481e+06 1.564e+04 -94.681 < 2e-16 ***
## total_sales -5.681e-06 1.139e-06 -4.989 6.06e-07 ***
## total_facs 6.362e+05 8.710e+02 730.379 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5207000 on 208802 degrees of freedom
## Multiple R-squared: 0.7188, Adjusted R-squared: 0.7188
## F-statistic: 2.668e+05 on 2 and 208802 DF, p-value: < 2.2e-16
neighbors <- dnearneigh(st_coordinates(county_sac.sf), 0, .05)
#run the Moran's I test for regression residuals
county_sac_list <- nb2listw(neighbors, style="W", zero.policy = TRUE)
lm.morantest(lm1, county_sac_list)
##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = total_emissions ~ total_sales + total_facs, data =
## county_sac)
## weights: county_sac_list
##
## Moran I statistic standard deviate = 388.5, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 4.306665e-01 -1.170494e-05 1.228922e-06
outlier_colors <- localmoran(county_sac$total_emissions,
county_sac_list,
TRUE,
na.action = na.omit)
meanVal <- mean(county_sac$total_emissions)
outlier_colors %<>% as_tibble() %>%
set_colnames(c("Ii","E.Ii","Var.Ii","Z.Ii","Pr(z > 0)")) %>%
#Create a new categorical value
#This variable will call insignificant units "insignificant"
#Outliers are determined by LI <0 - units higher than the mean are HL and units lower than the mean are LH
#Clusters are units with LI >= 0, and classified as HH or LL based on whether they are above or below the mean
mutate(coType = case_when(
`Pr(z > 0)` > 0.05 ~ "NA",
`Pr(z > 0)` <= 0.05 & Ii >= 0 & county_sac$total_emissions >= meanVal ~ "HH",
`Pr(z > 0)` <= 0.05 & Ii >= 0 & county_sac$total_emissions < meanVal ~ "LL",
`Pr(z > 0)` <= 0.05 & Ii < 0 & county_sac$total_emissions >= meanVal ~ "HL",
`Pr(z > 0)` <= 0.05 & Ii < 0 & county_sac$total_emissions < meanVal ~ "LH"))
county_sac$coType <- outlier_colors$coType %>% replace_na("NA") %>% as.factor()
counties_auto.sf <- merge(counties.sf, county_sac[c("county_name", "coType")], by = "county_name") %>%
st_transform("EPSG:4326") %>%
na.omit() %>%
select(-county_fips) %>%
select(-state_fips) %>%
select(-state_name) %>%
select(-state_abbv) %>%
select(-fips_class) %>%
group_by(county_name)
ggplot(counties_auto.sf) +
geom_sf(aes(fill=coType)) +
scale_fill_manual(values = c('red','brown','blue','cyan','white'), name='Clusters & \nOutliers') +
labs(title = "Spatial AC Map for County CO2 Emissions") +
theme_minimal()
