1. Coding

Q1

Bike Demand [3 points]

Hint: For [Bike Demand] part (a), do not prematurely summarize the data. We are expecting a line associated with each day.

[Bike Demand] This problem asks you to visualize a dataset(https://uwmadison.box.com/shared/static/f16jmkkskylfl1hnd5rpslzduja929g2.csv) of hourly bikeshare demand. Provide your code and make sure it is readable.

Q1a

Make a line plot of bike demand (count) by hour, faceted out across the 7 days of the week (weekday).

A1a.

library(ggplot2)
library(readr)

# Load the dataset
bike_data <- read_csv("https://uwmadison.box.com/shared/static/f16jmkkskylfl1hnd5rpslzduja929g2.csv")

# Adjust the date to be of the proper Date type
bike_data$dteday <- as.Date(bike_data$dteday)

# Create a line plot of bike demand by hour for each day, faceted by the day of the week
ggplot(bike_data, aes(x = hr, y = count, group = dteday)) +
  #each day in the dataset will have its unique line across the hours of that day
  geom_line(aes(color = dteday), alpha = 0.5, size = 0.7) +
  #facet_wrap(~ weekday, ncol = 1, scales = "free_y") +
  facet_wrap(~ weekday, ncol = 1) +
  labs(title = "Hourly Bike Demand for Each Individual Day by Weekday",
       x = "Hour of the Day",
       y = "Count of Bike Rentals") +
  theme_minimal() +
  theme(legend.position = "none")

Q1b

Create a new summary data.frame giving the 25 and 75 percent quantiles of demand (count) for each hour (hr) by day of the week (weekday) combination, separately within each year (yr) that the data was collected. A1b.

library(dplyr)

# Calculate the 25th and 75th percentiles for each hour and day of the week, within each year
summary_quantiles <- bike_data %>%
  group_by(yr, weekday, hr) %>%
  summarise(
    Q25 = quantile(count, 0.25),
    Q75 = quantile(count, 0.75),
    .groups = 'drop'  # Drop grouping structure after summarising
  )

# View the summary data frame
print(summary_quantiles)

## # A tibble: 336 × 5
##       yr weekday    hr   Q25   Q75
##    <dbl>   <dbl> <dbl> <dbl> <dbl>
##  1     0       0     0  39   108. 
##  2     0       0     1  33.5  79  
##  3     0       0     2  29    65.8
##  4     0       0     3  13.5  35  
##  5     0       0     4   3    11  
##  6     0       0     5   4     9  
##  7     0       0     6   3    16.8
##  8     0       0     7  10    37.2
##  9     0       0     8  33.5  92  
## 10     0       0     9  58   172. 
## # ℹ 326 more rows

Q1c

Using a ribbon plot, overlay the quantities from (b) onto your plot from part (a). Use color to distinguish between the ribbons for the first and second year that the data were collected.

A1c.

library(ggplot2)
library(readr)
library(dplyr)

# First, join the summary quantiles back to the original data for plotting
bike_data <- bike_data %>%
  left_join(summary_quantiles, by = c("yr", "weekday", "hr"))

# Plotting the line plot with ribbon overlay
ggplot(bike_data, aes(x = hr, y = count, group = dteday)) +
  geom_line(aes(color = as.factor(yr)), alpha = 0.3, size = 0.5) +  # Plot lines for each day
  geom_ribbon(aes(ymin = Q25, ymax = Q75, fill = as.factor(yr)), alpha = 0.2) +  # Add ribbons for quantiles
  #facet_wrap(~ weekday, ncol = 1, scales = "free_y") +
  facet_wrap(~ weekday, ncol = 1) +
  labs(title = "Hourly Bike Demand with Quantile Ribbons by Year and Weekday",
       x = "Hour of the Day",
       y = "Count of Bike Rentals",
       color = "Year",
       fill = "Year") +
  theme_minimal() +
  theme(legend.position = "bottom")  # Adjust legend position

Q1d

Provide a brief description of some takeaways from the final visualization.

A1d.

Takeways from the final visualization:
In general the final visualization overlays quantile ribbons onto the hourly bike demand data, faceted by weekday and distinguished by year. This visualization could help inform decisions on bike fleet management and infrastructure planning, and highlight the importance of considering temporal patterns in urban transportation planning.

(1) Demand Patterns: There is a clear cyclical pattern of bike demand across the hours of the day, which likely correlates with typical daily routines such as commuting in the morning and evening hours.

(2) Weekday vs. Weekend: 
i. The demand on weekdays shows pronounced peaks, especially during commuting hours (likely early morning hour7 and late afternoon hour17). 
ii. On weekends, the demand curve is smoother and the peaks are less pronounced, suggesting a more even distribution of bike usage throughout the day. We could also see the demand increases more gradually, peaks around midday, and then tapers off more slowly than on weekdays.

(3) Yearly Comparison: 
i. By examining the colors of the ribbons (with the two different years represented by different colors), we could see there is a general increase or decrease in demand from one year(year 0) to the next(year 1). 
ii.And the overlay shows consistency year over year, with some slight variations in the width of the ribbons, which indicates changes in variability.

(4) Demand Variability: The ribbons' width shows variability in demand. The certain hours having wider ribbons indicates more variability in bike demand during those hours. Conversely, narrower ribbons indicate more consistency in demand.

(5) Outliers: The individual lines, "year0 weekday2 hour14-16 are above the year1", that it s stray far from the bounds of the ribbon, which represent those days with unusually high demand in year0. 
These are points that significantly deviate from the quantile range, potentially indicating special events or weather impacts.

(6) Evening vs. Morning Demand:
The evening peak is wider or higher than the morning peak in weekdays, indicating that bike usage patterns are not symmetric throughout the day.

(7) Visual Clarity: The visualization seems clear, but with a lot of overlapping lines that could potentially obscure some details. It's important to note that when a lot of data is plotted simultaneously, the clarity of specific trends on individual days can be compromised.

Q2

Q2f

Geospatial Datasets [3 points].

For (f), do not use (a) or (c), since we’ve did those earlier.

[Geospatial Datasets] For each of the datasets below, specify whether it is in a vector or raster data format. If it is in a vector data format, explain which types of geometries it contains (e.g., a point or linestring). Explain your reasoning.

NOTE: I MOVE THE PART (a) to PART (e) to the Discussion part at about the bottom of the assignment!
Here is only the code part, which is PART (f).

Visualize one of these datasets.(choose from b.Africa Population 2020(https://github.com/krisrs1128/stat436_s24/raw/main/data/afripop2020.tif) OR d. Wisconsin EV Charging(https://github.com/krisrs1128/stat436_s24/raw/main/data/ev.geojson) OR e. Zion Elevation(https://github.com/krisrs1128/stat436_s24/raw/main/data/landsat.tif))

A2f.

# Visualize Wisconsin EV Charging Stations (https://github.com/krisrs1128/stat436_s24/raw/main/data/ev.geojson)

# Load the necessary libraries
library(sf)
library(ggplot2)
library(plotly)

# Read the GeoJSON data from the URL
ev_charging_url <- "https://github.com/krisrs1128/stat436_s24/raw/main/data/ev.geojson"
ev_data <- st_read(ev_charging_url, quiet = TRUE)

# Create a ggplot object
gg <- ggplot(data = ev_data) +
  geom_sf(aes(color = FUEL_TYPE_), size = 2) + 
  ggtitle("Wisconsin EV Charging Stations") +
  theme_minimal()

# Convert ggplot object to a plotly interactive plot
ggplotly(gg, tooltip = c("STATION_NA", "CITY", "FUEL_TYPE_"))

Q3

Political Book Recommendations [4 points]

Network and Hierarchical Data

[Political Book Recommendations] In this problem, we’ll study a network dataset of Amazon bestselling US Politics books. Books are linked by an edge if they appeared together in the recommendations (“customers who bought this book also bought these other books”).

Q3a

The code below reads in the edges and nodes associated with the network. The edges dataset(https://github.com/krisrs1128/stat992_f23/blob/6c4130bddbdfc9ef90537c794cdca47773643752/activities/week10/political-books-edges.csv) only contains IDs of co-recommended books, while the nodes data(https://github.com/krisrs1128/stat992_f23/blob/6c4130bddbdfc9ef90537c794cdca47773643752/activities/week10/political-books-nodes.csv) includes attributes associated with each book. Build a tbl_graph object to store the graph.

A3a.

library(tidygraph)
library(readr)
library(dplyr)

#Read the Data
edge_data_path <- "https://raw.githubusercontent.com/krisrs1128/stat992_f23/6c4130bddbdfc9ef90537c794cdca47773643752/activities/week10/political-books-edges.csv"
node_data_path <- "https://github.com/krisrs1128/stat992_f23/raw/6c4130bddbdfc9ef90537c794cdca47773643752/activities/week10/political-books-nodes.csv"


# edges <- read_csv(edge_data_path, col_types = cols(
#   source = col_character(),
#   target = col_character(),
#   weight = col_double()
# ))
# 
# nodes <- read_csv(node_data_path, col_types = cols(
#   id = col_character(),
#   label = col_character(),
#   political_ideology = col_character()
# ))

edges <- read_csv(edge_data_path, col_types = "cci")
nodes <- read_csv(node_data_path, col_types = "ccc")

#Create a tbl_graph
library(tidygraph)

# Convert nodes and edges to factors to ensure they align correctly
edges <- edges %>%
  mutate(
    source = as.factor(source),
    target = as.factor(target)
  )

nodes <- nodes %>%
  mutate(id = as.factor(id))

# Create the graph
graph <- tbl_graph(nodes = nodes, edges = edges, directed = FALSE)

#Examine the Graph
graph

## # A tbl_graph: 105 nodes and 441 edges
## #
## # An undirected simple graph with 1 component
## #
## # Node Data: 105 × 3 (active)
##    id    label                      political_ideology
##    <fct> <chr>                      <chr>             
##  1 0     1000 Years for Revenge     neutral           
##  2 1     Bush vs. the Beltway       conservative      
##  3 2     Charlie Wilson's War       conservative      
##  4 3     Losing Bin Laden           conservative      
##  5 4     Sleeping With the Devil    neutral           
##  6 5     The Man Who Warned America conservative      
##  7 6     Why America Slept          neutral           
##  8 7     Ghost Wars                 neutral           
##  9 8     A National Party No More   conservative      
## 10 9     Bush Country               conservative      
## # ℹ 95 more rows
## #
## # Edge Data: 441 × 3
##    from    to weight
##   <int> <int>  <int>
## 1     1     2      1
## 2     1     3      1
## 3     1     4      1
## # ℹ 438 more rows

Q3b

Use the result from part (a) to visualize the network as a node-link diagram. Include the book’s title in the node label, and shade in the node according to political ideology.

A3b.

#Color nodes according to the political_ideology attribute and label them with the book title.
library(ggraph)

# Convert nodes and edges to factors to ensure they align correctly
edges <- edges %>%
  mutate(
    source = as.factor(source),
    target = as.factor(target)
  )

nodes <- nodes %>%
  mutate(id = as.factor(id))

# Create the graph
graph <- tbl_graph(nodes = nodes, edges = edges, directed = FALSE)

# Visualize the network
ggraph(graph, layout = 'fr') + 
  geom_edge_link(color = 'grey50', alpha = 0.5) +
  #geom_node_point(aes(fill = political_ideology), size = 5) +
   geom_node_point(aes(color = political_ideology), size = 2) +
  geom_node_text(aes(label = label), repel = TRUE, 
                 color = 'black', size = 2, 
                 fontface = "bold", 
                 check_overlap = TRUE) +
  #scale_fill_manual(values = c('conservative' = 'red', 'liberal' = 'blue', 'neutral' = 'grey')) +
  theme_void() +
  ggtitle("Network of Political Books") +
  theme(legend.position = "bottom") +
  guides(fill = guide_legend(title = "Political Ideology"))

Q3c

Create the analogous adjacency matrix visualization. Provide examples of visual queries that are easy to answer using one encoding but not the other (i.e., what is easy to see in the node-link view vs. what is easy to see in the adjacency matrix).

A3c.

#An adjacency matrix from a graph object using heatmaps where cells are filled based on the presence or absence of edges between the nodes (books).
library(igraph)
library(ggplot2)

# Convert the graph to an adjacency matrix
adj_matrix <- as_adjacency_matrix(graph, attr = "weight", sparse = FALSE)

# Create a data frame suitable for ggplot2
adj_df <- reshape2::melt(adj_matrix)

# Give names to the variables
names(adj_df) <- c("source", "target", "weight")

# Plot the adjacency matrix
ggplot(adj_df, aes(x = source, y = target, fill = as.factor(weight))) +
  geom_tile(color = "white") +
  scale_fill_manual(values = c('1' = 'darkblue', '0' = 'white'), guide = FALSE) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1), 
        axis.title = element_blank()) +
  labs(fill = "Connection") +
  ggtitle("Adjacency Matrix of Political Books Network")

Visual Queries Comparison

Visual Queries Comparison:
(1) Node-Link Diagram (Part b):

i. Easy to see: The overall structure of the network, clusters of tightly connected nodes (books), and possible outliers.

ii.Difficult to see: Precise relationships between individual nodes when the network is large and dense, due to overlapping edges and nodes.

(2) Adjacency Matrix (Part c):

i. Easy to see: Direct comparisons of connections for specific nodes since rows and columns align with individual nodes. It's also easier to spot if a specific node (book) has many or few connections overall due to the row or column being predominantly filled or empty.And possibly the exact count of connections when the matrix isn't too large.

ii.Difficult to see: The overall structure and clustering of the network can be less intuitive in an adjacency matrix, especially for those not accustomed to this form of representation.

Q4

Colony Collapse, parts a – c [4 points]

[Colony Collapse] In this problem, we will study a dataset describing factors that might be leading to colony collapse disorder among bees. Since there are multiple stressors for each state × timepoint combination, we will use clustering to create a summarized “ecological stress” profile.

Q4a

The code below reads in data from its original Tidy Tuesday source and replaces NAs with 0’s.

We will summarize each state × timepoint combination by its profile of stress_pct across each stressor.

Reshape the data so that stressor appears along columns, with the associated percentage contained within the table.

Each row is therefore a profile of the amount of different stressors at a particular location and time.

A4a.

library(tidyverse)

# Load the data and replace NAs with 0s
stressor <- read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2022/2022-01-11/stressor.csv') %>%
  mutate(stress_pct = replace_na(stress_pct, 0))

# Reshape the data to wide format
stress_profile <- stressor %>%
  pivot_wider(
    names_from = stressor, 
    values_from = stress_pct,
    values_fill = list(stress_pct = 0) # Fill missing stress_pct with 0
  ) %>%
  # Group by year, months, and state to ensure each combination has a single row
  group_by(year, months, state) %>%
  # Summarize all stress_pct, if there are duplicates
  summarize(across(starts_with("Varroa mites"):starts_with("Unknown"), sum, na.rm = TRUE)) %>%
  ungroup() # Ungroup for further analysis

# Display the resulting data frame
print(stress_profile)

## # A tibble: 1,222 × 9
##     year months state `Varroa mites` Other pests/parasite…¹ Disesases Pesticides
##    <dbl> <chr>  <chr>          <dbl>                  <dbl>     <dbl>      <dbl>
##  1  2015 April… Alab…           16.7                   42.5       0          2.3
##  2  2015 April… Ariz…            8.4                   32.1       0.5       20.1
##  3  2015 April… Arka…           65.9                   63.4       7         22.6
##  4  2015 April… Cali…           50.9                   21.5       7.8       23  
##  5  2015 April… Colo…           38.6                    1.8       4.6        2.8
##  6  2015 April… Conn…           10.2                    0.4       1.9        0.8
##  7  2015 April… Flor…           44.2                   26.4       2.7       17.7
##  8  2015 April… Geor…           52.8                   17.5       2.9       28.7
##  9  2015 April… Hawa…           60.1                   61         0.3        0.1
## 10  2015 April… Idaho           25.5                    6.7       6.5        6  
## # ℹ 1,212 more rows
## # ℹ abbreviated name: ¹`Other pests/parasites`
## # ℹ 2 more variables: Other <dbl>, Unknown <dbl>

Q4b

Apply K-means to the profiles constructed in part (a). You may choose K. Visualize and briefly interpret the resulting centroids.

A4b.

library(tidyverse)
library(cluster)

# Assuming stress_profile is the dataframe from part a with the stressor percentages
set.seed(123)

# Before applying kmeans, ensure that the dataset has only numeric data for the clustering
# We'll select only the columns with stress percentages
numeric_data <- stress_profile %>% 
  select(`Varroa mites`, `Other pests/parasites`, `Disesases`, `Pesticides`, `Other`, `Unknown`) %>%
  na.omit() # Remove rows with NA values, if any

# Applying K-means clustering
k_means_result <- kmeans(numeric_data, centers = 3)

# Add the cluster assignments to your stress_profile data
stress_profile$cluster <- k_means_result$cluster

# Calculate the centroids of the clusters for interpretation
centroids <- as.data.frame(k_means_result$centers)
centroids$cluster <- seq_len(nrow(centroids))

# Visualize the centroids
centroids_long <- centroids %>%
  pivot_longer(-cluster, names_to = "stressor", values_to = "mean_stress_pct")

ggplot(centroids_long, aes(x = stressor, y = mean_stress_pct, group = cluster)) +
  geom_line(aes(color = as.factor(cluster))) +
  geom_point(aes(color = as.factor(cluster))) +
  theme_minimal() +
  labs(title = "Cluster Centroids of Ecological Stress Profiles",
       x = "Stressor",
       y = "Mean Stress Percentage",
       color = "Cluster")

# Print the resulting centroids for reference
print(centroids)

##   Varroa mites Other pests/parasites Disesases Pesticides     Other  Unknown
## 1     13.96603              4.261587  1.449524   2.373175  3.939524 2.983651
## 2     41.06180             12.510438  4.603967   8.523591  7.508977 4.593946
## 3     64.39027             38.223009  9.564602  17.947788 10.610619 7.077876
##   cluster
## 1       1
## 2       2
## 3       3

Briefly interpret the resulting centroids

Based on the centroid visualization and the data table, each cluster represents a different profile of ecological stress factors affecting bee colonies.

(1) Cluster 1 (Red): This group represents states with lower average stress percentages across all factors. The primary stressor is Varroa mites, followed by unknown factors. Other pests/parasites, diseases, pesticides, and other factors are relatively low. This cluster could represent states where bee colonies are relatively healthier or the management practices are effectively mitigating most stressors.

(2) Cluster 2 (Green): This cluster shows a moderate level of stress from Varroa mites, which is significantly higher than in Cluster 1 but lower than in Cluster 3. Other stress factors are present at moderate levels, with pesticides and unknown factors also contributing noticeably. This could indicate regions where multiple stressors are present, possibly interplaying to affect the health of bee colonies, but none are as dominant as Varroa mites in Cluster 3.

(3) Cluster 3 (Blue): States in this cluster are characterized by very high stress from Varroa mites, which is substantially higher than the other two clusters. The other significant stressor is other pests/parasites. Pesticides, diseases, and unknown factors also contribute to stress, but to a lesser degree. This profile is indicative of areas where Varroa mites and pests/parasites are a significant problem and likely the main contributors to colony collapse disorder.

Conclusion: 
These centroids highlight the different stressor profiles that could be useful for targeted interventions. For instance, efforts in Cluster 3 regions might focus on controlling Varroa mite infestations, while Cluster 2 might require a more balanced approach addressing several factors. The centroids provide valuable insights into which stressors are most prevalent in different clusters, guiding research and policy to better support bee colony health.

Q4c

Design and implement a visualization to show the change in cluster memberships over time. Briefly interpret your results in context.

A4c.

library(tidyverse)
library(lubridate)

# Load the data and replace NAs with 0s
stressor <- read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2022/2022-01-11/stressor.csv') %>%
  mutate(stress_pct = replace_na(stress_pct, 0))

# Reshape the data to wide format
stress_profile <- stressor %>%
  pivot_wider(
    names_from = stressor, 
    values_from = stress_pct,
    values_fill = list(stress_pct = 0) # Fill missing stress_pct with 0
  ) %>%
  group_by(year, months, state) %>%
  summarize(across(starts_with("Varroa mites"):starts_with("Unknown"), sum, na.rm = TRUE), .groups = 'drop')

# Set the seed and perform k-means clustering
set.seed(123)
k_means_result <- kmeans(stress_profile %>% select(-year, -months, -state), centers = 3)
stress_profile$cluster <- k_means_result$cluster

# Combine the year and months into a single 'time_period' column
stress_profile <- stress_profile %>%
  mutate(months = case_when(
    months == "January-March" ~ "01",
    months == "April-June" ~ "04",
    months == "July-September" ~ "07",
    months == "October-December" ~ "10"
  )) %>%
  unite(time_period, c("year", "months"), sep = "-") %>%
  mutate(time_period = as.Date(paste0(time_period, "-01")))

# Count the number of states in each cluster for each time period
cluster_membership <- stress_profile %>%
  group_by(time_period, cluster) %>%
  summarise(count = n(), .groups = 'drop')

# Plot the cluster memberships over time
ggplot(cluster_membership, aes(x = time_period, y = count, group = cluster, color = as.factor(cluster))) +
  geom_line() +
  labs(
    title = "Change in Cluster Memberships Over Time",
    x = "Time Period",
    y = "Count of States in each Cluster",
    color = "Cluster"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotate the x-axis labels for readability

# Print the graph
ggsave("cluster_membership_over_time.png", width = 12, height = 6, dpi = 300)

Briefly interpret the results in context

The above visualization shows the count of states in each of the three clusters over a series of time periods. Here's a brief interpretation of the results:

(1) Cluster 1 (Red): 

i. This cluster exhibits the most dynamic changes over time, with counts of states sometimes spiking dramatically. Such behavior could indicate acute episodes where the stressors corresponding to this cluster intensify significantly. For example, if Cluster 1 is characterized by higher levels of pesticide stress, the peaks could coincide with periods of increased pesticide application, such as during particular farming seasons. The volatility could also reflect changes in beekeeping practices, new legislation regarding pesticide use, or even the spread of diseases affecting bee colonies.

ii. (A spike in Cluster 1 might signal a widespread challenge that bee colonies face during specific periods, necessitating targeted interventions during those times.)

(2) Cluster 2 (Green): 

i. The variation in this cluster is moderate, suggesting that the stressors it represents affect a consistent number of states, but with some fluctuations. If Cluster 2 has stressors such as Varroa mites at moderate levels, the fluctuations could correspond to the effectiveness of mite control measures. These fluctuations might also reflect the success of interventions or natural cycles of stressor prevalence (like lifecycle patterns of parasites).

ii.The moderate fluctuations in Cluster 2 might imply that although certain stressors are consistently present, their impact is controlled to some extent, though not as effectively as in Cluster 3.

(3) Cluster 3 (Blue): The relatively stable and low count of states in this cluster suggests a specific combination of stressors that are either well-managed across these regions or inherently less prevalent. If, for instance, Cluster 3 represents areas less affected by human agricultural activity, this could explain why fewer states fall into this category and why their numbers don't change as dramatically. The stable low count may also be due to these states employing successful strategies to combat colony collapse disorder.

Q4d[Bonus]

[Bonus] The code below reads in spatial data associated with each state. Design and implement a visualization that encodes spatial variation across clusters, either at a single timepoint or over all time. An example result is shown in Figure 10.

A4d.

library(tidyverse)
library(rnaturalearth)
library(sf)
library(rmapshaper)

# Load the natural earth data for the states
states <- ne_states(country = "United States of America", returnclass = "sf") %>%
  st_as_sf() %>%
  ms_simplify() %>% # reduce size for easy plotting
  filter(!(name %in% c("Alaska", "Hawaii"))) %>% # excluding non-continental states
  select(name, geometry) %>%
  rename(state = name)

# stress_profile has a 'cluster' and 'state' column
# visualizing for a specific time period

# For this example, let's aggregate by the most recent cluster assignment per state
latest_cluster <- stress_profile %>%
  arrange(desc(time_period)) %>%
  group_by(state) %>%
  slice(1) %>% # Select the most recent time period per state
  ungroup()

# Join this with the spatial data
spatial_clusters <- states %>%
  left_join(latest_cluster, by = "state")

# Plotting the spatial variation across clusters
ggplot(spatial_clusters) +
  geom_sf(aes(fill = as.factor(cluster)), color = NA) + # Fill color by cluster, remove borders
  scale_fill_viridis_d(option = "C", begin = 0.3, end = 0.7, direction = 1, name = "Cluster") + # A diverging color palette
  theme_void() + # Remove extra elements for a clean map
  labs(title = "Spatial Variation Across Clusters") +
  theme(legend.position = "bottom")

2. Discussion

D1

Geospatial Datasets [3 points]

Note: Parts (a - e) will be graded as discussion and (f) as code.

D1a

NYC Building Footprints(https://uwmadison.box.com/shared/static/qfmrp9srsoq0a7oj0e7xmgu5spojr33e.geojson)

D1A:

library(sf)
NYC <- st_read("https://uwmadison.box.com/shared/static/qfmrp9srsoq0a7oj0e7xmgu5spojr33e.geojson")

## Reading layer `nyc-buildings' from data source 
##   `https://uwmadison.box.com/shared/static/qfmrp9srsoq0a7oj0e7xmgu5spojr33e.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 2726 features and 0 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.0048 ymin: 40.7032 xmax: -73.9648 ymax: 40.7432
## Geodetic CRS:  WGS 84

# Get the geometry types of the dataset
geometry_types <- st_geometry_type(NYC)

# Count occurrences of each geometry type
geometry_summary <- table(geometry_types)
print(geometry_summary)

## geometry_types
##           GEOMETRY              POINT         LINESTRING            POLYGON 
##                  0                  0                  0                  0 
##         MULTIPOINT    MULTILINESTRING       MULTIPOLYGON GEOMETRYCOLLECTION 
##                  0                  0               2726                  0 
##     CIRCULARSTRING      COMPOUNDCURVE       CURVEPOLYGON         MULTICURVE 
##                  0                  0                  0                  0 
##       MULTISURFACE              CURVE            SURFACE  POLYHEDRALSURFACE 
##                  0                  0                  0                  0 
##                TIN           TRIANGLE 
##                  0                  0

library(sf)
# Assume NYC is already read and contains MultiPolygon geometries
# Convert MultiPolygons to single Polygons if necessary
if (any(st_geometry_type(NYC) == "MULTIPOLYGON")) {
  NYC_single_polygons <- st_cast(NYC, "POLYGON")
}

# Now NYC_single_polygons should contain only Polygon geometries
geometry_types_single <- st_geometry_type(NYC_single_polygons)
geometry_summary_single <- table(geometry_types_single)
print(geometry_summary_single)

## geometry_types_single
##           GEOMETRY              POINT         LINESTRING            POLYGON 
##                  0                  0                  0               2752 
##         MULTIPOINT    MULTILINESTRING       MULTIPOLYGON GEOMETRYCOLLECTION 
##                  0                  0                  0                  0 
##     CIRCULARSTRING      COMPOUNDCURVE       CURVEPOLYGON         MULTICURVE 
##                  0                  0                  0                  0 
##       MULTISURFACE              CURVE            SURFACE  POLYHEDRALSURFACE 
##                  0                  0                  0                  0 
##                TIN           TRIANGLE 
##                  0                  0

(1) Vector data
The provided GeoJSON data for NYC Building Footprints represents vector data. 
(Since vector data is characterized by geometries such as points, lines, polygons (including multi-polygons), which represent discrete features on a map.)

It is evident that the geometries contained within this dataset are primarily Polygons, with at least one MultiPolygon. This is discerned from the "type": "Polygon" and "type": "MultiPolygon" properties within the geometry object of each Feature. These geometries represent the two-dimensional outlines of building footprints.

Moreover:

i. Each Feature represents a building, and the coordinates array within each Feature provides the vertices of the polygon that outlines the building's footprint. 

ii. A Polygon is typically used to represent areas, and in this case, it is used to delineate the spatial extent of individual buildings on a map. 

ii. The MultiPolygon type allows for representing a feature that consists of multiple polygons; this could be used for buildings that have disjointed parts, such as a complex of separate but related structures.

(2) Which types of geometries it contains

i. Polygon: A closed shape with linear edges that represent areas. In this NYC Building Footprints dataset, each polygon corresponds to the area covered by a building's footprint.

ii.MultiPolygon: A collection of polygons that represent a single feature. In mapping buildings, this could represent a group of buildings that are spaced closely together and function as a unit or complex.

D1b

Africa Population 2020(https://github.com/krisrs1128/stat436_s24/raw/main/data/afripop2020.tif)

D1B:

(1) Raster Data
The file extension .tif indicates that the dataset "Africa Population 2020" is a TIFF (Tagged Image File Format), which is commonly associated with raster data.

Raster data is a grid of pixels (or cells), where each pixel holds a value representing information, such as temperature, elevation, or in this case, population density. Raster files are often used for satellite images, photographs, and digital elevation models.

(2) To elaborate:
Raster Data: It is composed of pixels, with each pixel storing a numerical value. In the context of the "Africa Population 2020" dataset, each pixel represents the estimated population count or density in that area.

Since the file in question is a .tif and is named to suggest it contains population data, it is reasonable to deduce that the dataset is a raster representing population statistics across Africa, with the value in each pixel corresponding to population data for the geographic area that pixel covers.

D1c

Himalayan Glacial Lakes(https://raw.githubusercontent.com/krisrs1128/stat436_s23/main/data/GL_3basins_2015.topojson)

D1C:

library(sf)
HGL <- st_read("https://raw.githubusercontent.com/krisrs1128/stat436_s23/main/data/GL_3basins_2015.topojson")

## Reading layer `GL_3basins_2015' from data source 
##   `https://raw.githubusercontent.com/krisrs1128/stat436_s23/main/data/GL_3basins_2015.topojson' 
##   using driver `TopoJSON'
## Simple feature collection with 3609 features and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 80.07269 ymin: 27.43573 xmax: 88.7197 ymax: 30.62196
## CRS:           NA

# Get the geometry types of the dataset
geometry_types <- st_geometry_type(HGL)

# Count occurrences of each geometry type
geometry_summary <- table(geometry_types)
print(geometry_summary)

## geometry_types
##           GEOMETRY              POINT         LINESTRING            POLYGON 
##                  0                  0                  0               3609 
##         MULTIPOINT    MULTILINESTRING       MULTIPOLYGON GEOMETRYCOLLECTION 
##                  0                  0                  0                  0 
##     CIRCULARSTRING      COMPOUNDCURVE       CURVEPOLYGON         MULTICURVE 
##                  0                  0                  0                  0 
##       MULTISURFACE              CURVE            SURFACE  POLYHEDRALSURFACE 
##                  0                  0                  0                  0 
##                TIN           TRIANGLE 
##                  0                  0

(1) Vector Data
TopoJSON is an extension of GeoJSON, and it encodes topology in addition to the geometry and properties of GeoJSON. This means that it records lines and points shared between borders or edges to efficiently represent geometries and minimize redundancy. Similar to GeoJSON, TopoJSON files are used to encode vector data.

Based on the "type": "Topology" attribute and the presence of "arcs" and "features", we can conclude that this dataset contains vector data. 

(2) Which types of geometries it contains
Since the dataset is related to glacial lakes, which are specific areas of water collected in the glaciers' hollows, we can reasonably infer the geometries are likely to be either:
i. Polygons: To represent the area covered by each glacial lake. Polygons are used to represent two-dimensional shapes.
ii. MultiPolygons: If a single glacial lake has multiple disconnected bodies of water or if the dataset groups several close lakes into a single feature.

But according to the geometry_summary, the type of geometries is "Polygons".

D1d

Wisconsin EV Charging(https://github.com/krisrs1128/stat436_s24/raw/main/data/ev.geojson)

D1D:

# URL for the GeoJSON data
data_url <- "https://github.com/krisrs1128/stat436_s24/raw/main/data/ev.geojson"

# Read the GeoJSON file
ev_charging_data <- st_read(data_url, quiet = TRUE)

# Check geometry types in the dataset
geometry_types <- st_geometry_type(ev_charging_data)

# Summarize and print the geometry types
geometry_summary <- table(geometry_types)
print(geometry_summary)

## geometry_types
##           GEOMETRY              POINT         LINESTRING            POLYGON 
##                  0                846                  0                  0 
##         MULTIPOINT    MULTILINESTRING       MULTIPOLYGON GEOMETRYCOLLECTION 
##                  0                  0                  0                  0 
##     CIRCULARSTRING      COMPOUNDCURVE       CURVEPOLYGON         MULTICURVE 
##                  0                  0                  0                  0 
##       MULTISURFACE              CURVE            SURFACE  POLYHEDRALSURFACE 
##                  0                  0                  0                  0 
##                TIN           TRIANGLE 
##                  0                  0

(1) Vector Data
The "Wisconsin EV Charging" dataset is in vector data format. This is evident from the GeoJSON file format, indicated by the .geojson extension and the "type": "FeatureCollection" attribute within the data.

GeoJSON is a format for encoding various geographic data structures using JavaScript Object Notation (JSON) and is commonly used to represent vector geometries.

(2) Which types of geometries it contains
We can see that the "geometry" type for each feature is "Point." This suggests that the dataset consists of points, each representing a location of an electric vehicle (EV) charging station in Wisconsin. Here's a breakdown of the geometry contained:

i. Point: This is the simplest type of geometry, representing a single coordinate in space. In the context of this dataset, each point corresponds to the precise location of an EV charging station. The coordinates are given as a pair of longitude and latitude values.

The properties associated with each feature contain various details about the charging stations, such as the name, address, type of charging available, and contact information. The presence of latitude and longitude coordinates within the properties, and also within the geometry section, reinforces that this data is point-based vector data, each point corresponding to one EV charging station.

D1e

Zion Elevation(https://github.com/krisrs1128/stat436_s24/raw/main/data/landsat.tif)

D1E:

(1)Raster Data
The dataset "Zion Elevation" has the file extension .tif, which indicates that it is a TIFF file. TIFF files are commonly used for storing raster data, particularly in the context of geographic information systems (GIS) when referring to satellite imagery, aerial photography, and elevation data.

Raster data is made up of pixels (also called cells), each with an associated value that represents information for that area. In the case of an elevation dataset, each pixel would contain a value corresponding to the elevation at that point.

(2) To elaborate:
The file is likely a Digital Elevation Model (DEM), which is a raster representation of a terrain's surface created from terrain elevation data.

"landsat.tif" for Zion Elevation is a raster dataset where each pixel's value represents the elevation above sea level for that specific geographic location. The naming "Landsat" also implies it could be derived from satellite data, Landsat being a program that collects satellite imagery of the Earth.

PART (f) is in the "code" part at about the middle of the assignment.

D2

Olympics Interactive App [3 points] [Olympics Interactive App] The code below sets up a Shiny app for interactively visualizing athlete weight and heights in the 2012 London Olympics. We would like to have an interactive scatterplot of Weight vs. Height, cm that updates which points (athletes) are highlighted depending on which sports have been selected by a dropdown menu. Code for generating the scatterplot is provided in the function scatterplot.

D2a

Provide server code which would allow the scatterplot to update the highlighted athletes depending on the currently selected sports.

D2A:

library(shiny)
library(tidyverse)
library(readr)

# Load the data
olympics_url <- "https://uwmadison.box.com/shared/static/rzw8h2x6dp5693gdbpgxaf2koqijo12l.csv"
olympics <- read_csv(olympics_url)

#' Scatterplot with highlighted points
scatterplot <- function(df) {
  ggplot(df, aes(x = Weight, y = `Height, cm`, alpha = as.numeric(selected), size = as.numeric(selected))) +
    geom_point() +
    scale_alpha(range = c(0.05, .8)) +
    scale_size(range = c(0.1, 1)) +
    
    labs(title = "Height vs. Weight of Athletes", x = "Weight (kg)", y = "Height (cm)") +
    theme_minimal()
}

ui <- fluidPage(
  titlePanel("Olympic Athletes - 2012 London Olympics"),
  sidebarLayout(
    sidebarPanel(
      selectInput("dropdown", "Select a Sport", choices = unique(olympics$Sport), multiple = TRUE)
    ),
    mainPanel(
      plotOutput("scatterplot")
    )
  )
)

server <- function(input, output) {
  # fill this in...
  filtered_data <- reactive({
    if (length(input$dropdown) == 0) {
      olympics %>% mutate(selected = FALSE)
    } else {
      olympics %>% 
        mutate(selected = Sport %in% input$dropdown)
    }
  })

  output$scatterplot <- renderPlot({
    scatterplot(filtered_data())
  })
}

shinyApp(ui = ui, server = server)

Shiny applications not supported in static R Markdown documents

D2b

We have been asked to also print a table of the selected athletes. Assume the UI has the form,

ui <- fluidPage(
  selectInput("dropdown", "Select a Sport", choices = unique(olympics$Sport), multiple = TRUE),
  plotOutput("scatterplot"),
  dataTableOutput("table")
)

# Updated UI:
library(shiny)
library(tidyverse)
library(DT)  # Include DT for dataTableOutput

ui <- fluidPage(
  titlePanel("Olympic Athletes - 2012 London Olympics"),
  sidebarLayout(
    sidebarPanel(
      selectInput("dropdown", "Select a Sport", choices = unique(olympics$Sport), multiple = TRUE)
    ),
    mainPanel(
      plotOutput("scatterplot"),
      DTOutput("table")  # Data table output for displaying selected athletes
    )
  )
)

server <- function(input, output) {
  filtered_data <- reactive({
    if (length(input$dropdown) == 0) {
      olympics %>% mutate(selected = FALSE)
    } else {
      olympics %>% 
        mutate(selected = Sport %in% input$dropdown)
    }
  })

  
  #Updated Server:
  # Render the scatterplot
  output$scatterplot <- renderPlot({
    scatterplot(filtered_data())
  })

  # Render the data table
  output$table <- renderDT({
    # Select only the relevant columns to display in the table
    filtered_data() %>%
      filter(selected) %>%
      select(Name, Country, Sport, `Height, cm`, Weight, Age)
  }, options = list(pageLength = 10))
}

shinyApp(ui = ui, server = server)

Shiny applications not supported in static R Markdown documents

Describe changes to your solution to (a) to meet the new requirements. How would you minimize code duplication? Be as specific as possible.

D2B:

To minimize code duplication, I use the reactive expression I defined earlier (filtered_data) for both the scatterplot and the data table. This ensures that both visual components use the same filtered dataset, maintaining consistency across the app and reducing the need to duplicate data filtering logic.

Description of Changes:
(1) UI Components: The UI now includes an additional component (DTOutput("table")) to display the data table. I use the DT library, which needs to be loaded with library(DT).

(2) Server Logic: The server function now also handles rendering a data table. The renderDT function is used to render the table of filtered athletes who are selected based on the dropdown. I filter the filtered_data reactive data to include only selected athletes and then specify the columns that we want to display to make the table more readable and relevant.

(3) Reusing filtered_data: By using the same filtered_data() reactive expression for both the plot and the table, we avoid duplicating the logic for filtering data based on selected sports. This ensures that both outputs are always in sync and minimizes the amount of code.

(4) Data Table Configuration: The renderDT function is configured to display only 10 rows at a time by setting options = list(pageLength = 10).


Conclusion:
This structured approach maintains a clean separation of concerns within the app, ensuring that data manipulation is centralized and visual components are solely responsible for presentation. This design reduces redundancy and increases maintainability.

Homework 3

Linlin Zhang

1. Coding

Q1

Q1a

Q1b

Q1c

Q1d

Q2

Q2f

Q3

Q3a

Q3b

Q3c

Q4

Q4a

Q4b

Q4c

Q4d[Bonus]

2. Discussion

D1

D1a

D1b

D1c

D1d

D1e

D2

D2a

D2b