2022_safegraph

library(sf)

## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

library(tmap)

## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)

library(data.table)

setwd("C:\\Users\\zhaiw\\Desktop\\safegraph_anlaysis")

shapefile <- st_read("./bexar/bexar_cbg.shp")

## Reading layer `bexar_cbg' from data source 
##   `C:\Users\zhaiw\Desktop\safegraph_anlaysis\bexar\bexar_cbg.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1084 features and 16 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -98.80655 ymin: 29.11443 xmax: -98.1169 ymax: 29.76071
## Geodetic CRS:  NAD83

df_May <- fread('./2022_data/df_May.txt')
df_June <- fread('./2022_data/df_June.txt')
df_Aug <- fread('./2022_data/df_Aug.txt')

df <-merge(df_May, df_June, by = "GEOID", all = TRUE)
df <-merge(df, df_Aug, by = "GEOID", all = TRUE) 

df$Cooler <- rowMeans(df[, c("May_4", "June_28", "Aug_24")],na.rm = TRUE)
df$Warmer <- rowMeans(df[, c("May_18", "June_14", "Aug_10")],na.rm = TRUE)

shapefile_df <- merge(shapefile, df, by = "GEOID")

cooler_map<- tm_shape(shapefile_df) +tm_fill(col = "Cooler")+ tm_layout(title = "Cooler Trip Map")
warmer_map<- tm_shape(shapefile_df) +tm_fill(col = "Warmer")+ tm_layout(title = "Warmer Trip Map")
council_d <-st_read("./Redistricted_Council_Districts_2022/Redistricted_Council_Districts_2022/RedistrictedCouncilDistricts2022.shp")

## Reading layer `RedistrictedCouncilDistricts2022' from data source 
##   `C:\Users\zhaiw\Desktop\safegraph_anlaysis\Redistricted_Council_Districts_2022\Redistricted_Council_Districts_2022\RedistrictedCouncilDistricts2022.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 10 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -10999870 ymin: 3399373 xmax: -10934130 ymax: 3469007
## Projected CRS: WGS 84 / Pseudo-Mercator

cooler_map <-cooler_map+ tm_shape(council_d) + tm_polygons(border.col = "black", alpha = 0)
warmer_map <-warmer_map+ tm_shape(council_d) + tm_polygons(border.col = "black", alpha = 0)

# Create a layout with two maps in a row
layout <- tmap_arrange(
  cooler_map, warmer_map,
  nrow = 1, ncol = 2,  # Arrange the maps in a single row
  widths = c(0.5, 0.5)  # Adjust the width of each map
)

# Plot the maps side by side
tmap_options(check.and.fix = TRUE) 
map0 <- tmap_arrange(layout = layout)
map0

#tmap_save(map0, filename = "two_maps.png", width = 10, height = 5, dpi = 300)

shapefile_df$Percentage <- 100*(shapefile_df$Warmer -shapefile_df$Cooler)/shapefile_df$Cooler

#df$Percentage <- 100*(df$Warmer -df$Cooler)/df$Cooler
#write.csv(df, file = "merged_df.csv", row.names = FALSE)


council_d <-st_read("./Redistricted_Council_Districts_2022/Redistricted_Council_Districts_2022/RedistrictedCouncilDistricts2022.shp")

## Reading layer `RedistrictedCouncilDistricts2022' from data source 
##   `C:\Users\zhaiw\Desktop\safegraph_anlaysis\Redistricted_Council_Districts_2022\Redistricted_Council_Districts_2022\RedistrictedCouncilDistricts2022.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 10 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -10999870 ymin: 3399373 xmax: -10934130 ymax: 3469007
## Projected CRS: WGS 84 / Pseudo-Mercator

map <- tm_shape(shapefile_df) +tm_fill(col = "Percentage",breaks = c(-100, -50, -25, 0, 25, 50, 100))+ 
  tm_layout(title = "Trip change (%)")
map <- map +
  tm_shape(council_d) +
  tm_polygons(border.col = "black", alpha = 0)

map

## Variable(s) "Percentage" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

#tmap_save(map, filename = "map_output.png", width = 5, height = 5, dpi = 300)

library(ggplot2)
library(sf)
UTFVI <- fread ('utfvi_BGs_all.csv')
UTFVI$GEOID <- as.character(UTFVI$GEOID)
merge_df <- merge(UTFVI, shapefile_df, by= 'GEOID')

correlation <- cor.test(merge_df$utfvi_BGs, merge_df$Percentage)

my_plot <- ggplot(merge_df, aes(x = utfvi_BGs, y = Percentage)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "UTFVI", y = "Mobility change (%)")

my_plot <- my_plot +
  annotate("text", x = Inf, y = Inf, hjust = 1, vjust = 1,
           label = paste("Correlation:", round(correlation$estimate, 2)),
           size = 5) +
  annotate("text", x = Inf, y = Inf, hjust = 1, vjust = 3,
           label = paste("p-value:", format(correlation$p.value, scientific = TRUE)),
           size = 5)
my_plot

## `geom_smooth()` using formula = 'y ~ x'

#ggsave("plot_output.png", plot = my_plot, width = 5, height = 4, dpi = 300)